Research Article

Bayesian inference for biophysical neuron models enables stimulus optimization for retinal neuroprosthetics

Institute for Ophthalmic Research, University of Tübingen, Germany
Naturwissenschaftliches und Medizinisches Institut an der Universität Tübingen, Germany
Center for Integrative Neuroscience, University of Tübingen, Germany
Bernstein Center for Computational Neuroscience, University of Tübingen, Germany
Department of Neuroscience, University of Pennsylvania, United States
Institute for Bioinformatics and Medical Informatics, University of Tübingen, Germany

Oct 27, 2020

Open access
Copyright information

Abstract
Introduction
Materials and methods
Results
Discussion
Appendix 1
Data availability
References
Article and author information
Metrics

Abstract

While multicompartment models have long been used to study the biophysics of neurons, it is still challenging to infer the parameters of such models from data including uncertainty estimates. Here, we performed Bayesian inference for the parameters of detailed neuron models of a photoreceptor and an OFF- and an ON-cone bipolar cell from the mouse retina based on two-photon imaging data. We obtained multivariate posterior distributions specifying plausible parameter ranges consistent with the data and allowing to identify parameters poorly constrained by the data. To demonstrate the potential of such mechanistic data-driven neuron models, we created a simulation environment for external electrical stimulation of the retina and optimized stimulus waveforms to target OFF- and ON-cone bipolar cells, a current major problem of retinal neuroprosthetics.

Introduction

Mechanistic models have been extensively used to study the biophysics underlying information processing in single neurons and small networks in great detail (Gerstner and Kistler, 2002; Koch, 2004). In contrast to phenomenological models used for neural system identification, such models try to preserve certain physical properties of the studied system to facilitate interpretation and a causal understanding. For example, biophysical models can incorporate the detailed anatomy of a neuron (Golding et al., 2001; Poirazi et al., 2003; Hay et al., 2011), its ion channel types (Hodgkin and Huxley, 1952; Fohlmeister and Miller, 1997) and the distributions of these channels (Rattay et al., 2017) as well as synaptic connections to other cells (O'Leary et al., 2014). For all these properties, the degree of realism can be adjusted as needed. While increased realism may enable models to capture the highly non-linear dynamics of neural activity more effectively, it usually also increases the number of model parameters. While the classical Hodgkin-Huxley model with one compartment has already 10 free parameters (Hodgkin and Huxley, 1952), detailed multicompartment models of neurons can have dozens or even hundreds of parameters (Taylor et al., 2009; Hay et al., 2011).

Constraining many of these model parameters such as channel densities requires highly specialized and technically challenging experiments, and hence it is usually not viable to measure every single parameter for a neuron model of a specific neuron type. Rather, parameters for mechanistic simulations are often aggregated over different neuron types and even across species. Even though this may be justified in specific cases, it likely limits our ability to identify mechanistic models of individual cell types. Alternatively, parameter search methods have been proposed to identify the parameters of mechanistic neuron models from standardized patch-clamp protocols based on exhaustive grid-searches (Goldman et al., 2001; Prinz et al., 2003; Stringer et al., 2016) or evolutionary algorithms (Gerken et al., 2006; Keren et al., 2005; Achard and De Schutter, 2006; Rossant et al., 2011). Such methods are often inefficient and identify only a single point estimate consistent with the data (for discussion, see Gonçalves et al., 2020).

Here, we built on recent advances in Bayesian simulation-based inference to fit multicompartment models of neurons with realistic anatomy in the mouse retina. We used a framework called Sequential Neural Posterior Estimation (SNPE) (Lueckmann et al., 2017; Gonçalves et al., 2020) to identify model parameters based on high-throughput two-photon measurements of these neurons’ responses to light stimuli. SNPE is a Bayesian simulation-based inference algorithm that allows parameter estimation for simulator models for which the likelihood cannot be evaluated easily. The algorithm estimates the distribution of model parameters consistent with specified target data by evaluating the model for different sets of parameters and comparing the model output to the target data. To this end, parameters are drawn from a prior distribution, which is an initial guess about which parameters are likely to produce the desired model output. For example, the choice of prior distribution can be informed by the literature, without constraining the model to specific values. The model output for the sampled parameter sets can then be used to refine the distribution over plausible parameters given the data. This updated distribution, containing information from both the prior and the observed simulations, is known as the posterior. For high-dimensional parameter spaces, many samples are necessary to obtain an informative posterior estimate. Therefore, to make efficient use of simulation time, SNPE iteratively updates its sampling distribution, such that only in the first round samples are drawn from the prior, while in subsequent rounds samples are drawn from intermediate posteriors. This procedure increases the fraction of samples leading to simulations close to the target data. Since this approach for parameter estimation not only returns a point-estimate but also a posterior distribution over parameters consistent with the data, it allows one to straightforwardly determine how well the parameters are constrained. While the method has been used previously to fit simple neuron models (Lueckmann et al., 2017; Gonçalves et al., 2020), it has so far not been applied to models as complex and realistic as the ones presented here.

We estimated the posterior parameter distribution of multicompartment models of three retinal neurons, a cone photoreceptor (cone), an OFF- and an ON-bipolar cell (BC). The structure of the BC models was based on high-resolution electron microscopy reconstructions (Helmstaedter et al., 2013) and in six independently parameterized regions. We performed parameter inference based on the responses of these neurons to standard light stimuli measured with two-photon imaging of glutamate release using iGluSnFR as an indicator (Franke et al., 2017). Our analysis shows that many of the model parameters can be constrained well, yielding simulation results consistent with the observed data. After validating our model, we show that the inferred models and the inference algorithm can be used to efficiently guide the design of electrical stimuli for retinal neuroprosthetics to selectively activate OFF- or ON-BCs. This is an important step toward solving a long-standing question in the quest to provide efficient neuroprosthetic devices for the blind.

Materials and methods

Key resources table

Reagent type (species) or resource	Designation	Source or reference	Identifiers	Additional information
Genetic reagent (mouse)	B6;129S6-Chat^{tm2(cre)LowlJ}	Jackson laboratory JAX 006410	RRID:IMSR_JAX:006410
Genetic reagent (mouse)	Gt(ROSA)26Sor ^{tm9(CAG-tdTomato)Hze}	Jackson laboratory JAX 007905	RRID:IMSR_JAX:007905
Strain (mouse, female)	B6.CXB1-Pde6b^rd10	Jackson laboratory JAX 004297	RRID:IMSR_JAX:004297
Strain (Adeno-associated virus)	AAV2.7m8.hSyn.iGluSnFR	Virus facility, Institute de la Vision, Paris
Software, algorithm	NeuronC	https://retina.anatomy.upenn. edu/rob/neuronc.html	RRID:SCR_014148	Version 6.3.14
Software, algorithm	SNPE	https://github.com/mackelab/delfi		See Inference algorithm
Software, algorithm	COMSOL Multiphysics	COMSOL Multiphysics	RRID:SCR_014767

Biophysical neuron models

We created detailed models of three retinal cell types: a cone, an ON- (Figure 1A,Bi) and an OFF-BC (Figure 1Bii). From the different OFF- and ON-BC types, we chose to model the types 3a and 5o, respectively, because those were the retinal cone bipolar cell (CBC) types in mice for which we could gather most information. To model the light response, the OFF-BC model received input from five and the ON-BC from three cones (Behrens et al., 2016). Every cone made two synaptic connection with each BC. The postsynaptic conductances were set to 0.25 nS per connection.

Figure 1

Download asset Open asset

Multicompartment models

Request a detailed protocol

We used NeuronC (Smith, 1992) to implement multicompartment models of these neurons. A multicompartment model subdivides a neuron into a finite number of compartments. Every compartment is modeled as an electrical circuit, has a position in space, a spatial shape and is connected to at least one neighboring compartment (Figure 1C). The voltage in a compartment n, connected to the compartments $n - 1$ and $n + 1$ is described by:

\frac{δ}{δ t} V_{m}^{n} = \frac{1}{c_{m}^{n}} (\frac{V_{m}^{n + 1} - V_{m}^{n}}{r_{i}^{n + 1} + r_{i}^{n}} + \frac{V_{m}^{n - 1} - V_{m}^{n}}{r_{i}^{n - 1} + r_{i}^{n}} + \frac{V_{r}^{n} + V_{e x}^{n} - V_{m}^{n}}{r_{m}^{n}} + \sum_{e} \frac{V_{e}^{n} + V_{e x}^{n} - V_{m}^{n}}{r_{e}^{n} (. . .)}) + \frac{δ}{δ t} V_{e x}^{n} .

Here, compartments are either modeled as cylinders or spheres. The membrane capacitance $c_{m}^{n}$ , membrane resistance $r_{m}^{n}$ and axial resistance $r_{i}^{n}$ of a compartment n are assumed to be dependent on the compartment surface area $A_{m}^{n}$ and/or the compartment length $l_{c}^{n}$ :

r_{m}^{n} = R_{m} / A_{m}^{n}, c_{m}^{n} = C_{m} A_{m}^{n}, r_{i}^{n} = R_{i} l_{c}^{n} / A_{m}^{n} .

We assumed the specific membrane resistance $R_{m}$ , the specific membrane capacitance $C_{m}$ , and the axial resistivity $R_{i}$ to be constant over all compartments within a cell model. We used $R_{i} = 132 Ω c m$ for all cell types and informed our priors for $C_{m}$ and $R_{m}$ , which we estimated for every cell type individually, based on estimates for rod bipolar cells of rats (Oltedal et al., 2009). Parameters of NeuronC and the used values are summarized in Appendix 1—table 3.

Anatomy

Request a detailed protocol

We used a simplified cone morphology consisting of four compartments: one cone-shaped compartment for the outer segment, one spherical compartment for the combination of inner segment and soma, one cylindrical compartment for the axon, and another spherical one for the axonal terminals (Figure 1). The light collecting area in the outer segment was set to 0.2 µm² (Nikonov et al., 2006). The diameter of the soma $d_{S}^{c}$ , the axon $d_{A}^{c}$ and axonal terminals $d_{A T}^{c}$ , the length of the axon $l_{A}^{c}$ and the length of the outer segment $l_{O S}^{c}$ were based on electron microscopy data (Carter-Dawson and LaVail, 1979):

d_{S}^{c} = 5.13 μ m, d_{A}^{c} = 1.3 μ m, d_{A T}^{c} = 6 μ m, l_{A}^{c} = 15 μ m, l_{O S}^{c} = 14.4 μ m .

The BC morphologies in this study were based on serial block-face electron microscopy data of mouse bipolar cells (Helmstaedter et al., 2013). We extracted the raw voxel-based morphologies from the segmentation of the EM dataset and transformed them into a skeleton plus diameter representation using Vaa3D-Neuron2 auto tracing (Xiao and Peng, 2013). These where then manually refined using Neuromantic (Myatt et al., 2012) to correct errors originating from small segmentation errors (Figure 1).

The ON-BC morphology we chose was classified as type 5o, equal to the functional type of the model. For the OFF-BC, we decided for a morphology classified as type 3b, although we functionally modeled a type 3a cell, because the chosen reconstructed morphology was of higher quality than all available type 3a reconstructions and because type 3a and 3b BCs have very similar morphologies. Additionally, type 3a and 3b mostly differ in the average axonal field size (Wässle et al., 2009; Helmstaedter et al., 2013), with that of type 3a being larger than that of type 3b. The selected morphology had the largest axonal field among all cells classified as 3b in the dataset, well within the range of type 3a cells.

Because the computational time scales approximately linear with the number of BC compartments, using the full number of compartments of the EM reconstructions (>1000) during parameter inference was computationally infeasible. Therefore, we utilized the compartment condensation algorithm of NeuronC, which iteratively reduces the number of compartments while preserving biophysical properties (Smith, 1992). To be able to draw a sufficient number of samples, we reduced the number of compartments during parameter inference to 22 and 19 for the OFF- and ON-BC respectively (requiring approximately 3 min per simulation for a 31 s light stimulus). To simulate the electrical stimulation, more compartments are necessary to capture the effect of the electrical field on the neurites of the BC models. Therefore, we increased the number of compartments to 139 and 152 for the OFF- and an ON-BC, respectively, which is sufficient to accurately represent all major neurites without becoming computationally too expensive (requiring approximately 20 min and 13 min per simulation for a 31 s light stimulus for the OFF- and ON-BC, respectively).

Ion channels and synapses - literature review

Request a detailed protocol

The complement and distribution of voltage- or ligand-gated ion channels shapes the response of neurons. Here, ion channels are modeled as additional electrical elements in the compartments’ membrane with conductances dependent on time varying parameters, such as the membrane potential and the calcium concentration within the cell. In addition to the equations that govern a channel’s kinetics, their location in the cell has to be defined. After a literature review of retinal cone bipolar cell types in mice, we decided to model the OFF- and ON-type for which we could gather most information, namely BC3a and BC5. Currently, there are three accepted subtypes of BC5: 5o, 5i and 5t (Greene et al., 2016). Here, we modeled the BC5 subtype that expresses voltage-gated sodium channels (Hellmer et al., 2016) which probably also corresponds to the more transient BC5 subtype reported in Ichinose et al., 2014. The TTX sensitivity observed in Matsumoto et al., 2019 suggests that both, 5o and 5i express voltage-gated sodium channels. To make our model consistent, we used data from the same BC5 subtype (5o) for the morphology, the target data and the number of cone contacts. A summary of all used channels, their location within the models and the respective references can be found in Table 1. The following paragraphs describe which channels were included in the models and why. Note, however, that for all channels (except the L-type calcium channel in the axon terminals, as calcium channels are necessary in the model for neurotransmitter release) channel densities of zero were included in the prior distributions, thereby allowing the parameter inference to effectively remove ion channels from the model.

Table 1

Ion channels of biophysical models.

Channel	Cone	OFF-BC (type 3a)	ON-BC (type 5)	Cone references	BC references
${C a}_{L}$	AT	S, AT	S, AT	Morgans et al., 2005; Mansergh et al., 2005	Van Hook et al., 2019
${C a}_{T}$		S, AT			Van Hook et al., 2019
${C a}_{P}$	AT	S, AT	S, AT	Morgans et al., 1998	Morgans et al., 1998
${H C N}_{1}$	All		D, S, AT	Knop et al., 2008; Van Hook et al., 2019	Knop et al., 2008; Hellmer et al., 2016
${H C N}_{4}$		D, S, AT			Hellmer et al., 2016; Knop et al., 2008
$K_{V}$	IS/S	DD, PA, DA	DD, PA, DA	Knop et al., 2008; Van Hook et al., 2019	Ma et al., 2005
$K_{i r}$		S	S		Cui and Pan, 2008; Knop et al., 2008
${C l}_{C a}$	AT			Yang et al., 2008; Caputo et al., 2015
${N a}_{V}$		DA	DA		Hellmer et al., 2016

Regions of ion channels and the respective abbreviations as in Figure 1.

D refers to the combination of DD and PD. All refers to the combination of IS/S, A, and AT.
If multiple regions are stated for a neuron, the ion channel density differs between them.

In their axon terminals, cones express L-type calcium ( ${C a}_{L}$ ) channels that mediate release of the transmitter glutamate (Morgans et al., 2005; Mansergh et al., 2005; Ingram et al., 2020). We modeled calcium extrusion purely with calcium pumps ( ${Ca}_{P}$ ) since other mechanisms such as sodium-calcium-exchangers probably only play a minor functional role in cones (Morgans et al., 1998). Additionally, there is evidence that cones express hyperpolarization-activated cyclic nucleotide-gated cation (HCN) channels of the type 1, mostly in the inner segment but also in the axon (Knop et al., 2008; Van Hook et al., 2019). The presence of ${H C N}_{3}$ channels in mouse cones is more controversial. These channels have been observed in rat cones (Müller et al., 2003), and a more recent study also found evidence for ${H C N}_{3}$ channels at the synaptic terminals of mouse cones, but could not observe any functional differences between wild-type and ${H C N}_{3}$ -knockout mice. To restrict the number of model parameters, we did not include ${H C N}_{3}$ in our cone model. However, we added calcium-activated chloride ( ${C l}_{C a}$ ) channels to the axon terminals (Yang et al., 2008; Caputo et al., 2015) and voltage-gated potassium channels $K_{V}$ at the inner segment (Van Hook et al., 2019).

Our BC5 type expresses voltage-gated sodium ( ${N a}_{V}$ ) channels at the axon shaft (Hellmer et al., 2016). Another study found inward-rectifier potassium ( $K_{i r}$ ) channels at the soma of BC5 (Knop et al., 2008), which were also found in the homologous type in rat (Cui and Pan, 2008). Additionally, BC5 express HCN channels at the axon terminal, the soma and the dendrites (Knop et al., 2008; Hellmer et al., 2016). From the four subtypes of HCN, BC5 seem to almost exclusively express ${H C N}_{1}$ . In the rat, there is also evidence for the expression of ${H C N}_{4}$ channels in BC5 (Müller et al., 2003; Ivanova and Müller, 2006), but this could not be verified for mice. Data from rat suggests that BCs with ${N a}_{V}$ channels also express $K_{V}$ channels (Ma et al., 2005). We therefore added $K_{V}$ channels at the dendrites and the axon.

Similar to BC5, BC3a express HCN channels at the axon terminals, the soma and the dendrites. However, instead of ${H C N}_{4}$ they express ${H C N}_{4}$ (Hellmer et al., 2016; Knop et al., 2008). There is also evidence that BC3a express ${N a}_{V}$ channels at the axon shaft (Hellmer et al., 2016), which were also found in the homologous type in rat (Cui and Pan, 2008). Just like for BC5, we added also $K_{V}$ in BC3a were only reported for rat so far (Cui and Pan, 2008). As we could not find any evidence for the lack of $K_{i r}$ channels in mouse BC3a and the channel repertoires of BC3a in mouse and rat are overall very consistent, we included them in our model.

The distribution of calcium channels in mouse CBCs is largely unknown (Van Hook et al., 2019). In the rat retina, there is evidence for T-type calcium ( ${C a}_{T}$ ) channels in BC3a (Ivanova and Müller, 2006). Calcium currents of unspecified type were observed in BC5 (Cui and Pan, 2008). Generally, L-type calcium ( ${C a}_{L}$ ) channels are believed to mediate neurotransmitter release in almost all BCs across types and species (Van Hook et al., 2019). Therefore, we included them in both BC models. The literature review in Van Hook et al., 2019 suggests that T-type calcium channels might be exclusively expressed in BC3. In mouse BC3b, the simultaneous expression of both ${C a}_{T}$ and ${C a}_{L}$ has been described (Cui et al., 2012). Furthermore, the latter and other studies (Hu et al., 2009; Satoh et al., 1998) suggest that voltage-gated calcium channels might not be located in the axon terminals only, but also in the soma and might play a role in signal transmission within the cell. Based on the studies mentioned, we assumed that BC3a and BC5 express ${C a}_{L}$ in the axon terminals and potentially also at the soma. The BC3a model may additionally use ${C a}_{T}$ channels, both at the soma and at the axon terminals. For calcium extrusion, we added calcium pumps (Morgans et al., 1998).

BC5 receive input from cones via the metabotropic glutamate receptor 6 (mGluR6) (Van Hook et al., 2019). BC3a receive input from cones via kainate receptors (Ichinose and Hellmer, 2016). We modeled the kainate receptors by modifying the inactivation time constant $τ_{γ}$ of the AMPA receptors included with NeuronC.

Ion channels and synapses - implementation

Request a detailed protocol

All ion channels in this study were based on the models available in NeuronC. We used both Hodgkin-Huxley (HH) and Markov-Sequential-State (MS) channel implementations. Since we did not add channel noise to our model, every HH channel could have also been described as an equivalent MS channel. However, since HH channels are computationally less expensive, we used HH implementations wherever possible. Implementation details and references are listed in Table 2. The L-type calcium channel, for example, was based on the HH model defined by the following equations:

\frac{1}{r_{e}} := g_{e} = c^{3} \cdot g_{max}, \frac{δ}{δ t} c = (1 - c) \cdot α (V) - c \cdot β (V),

α (V) = η_{T} \cdot \frac{- 0.04 \cdot (V + 15)}{\exp (- 0.04 \cdot (V + 15)) - 1} \cdot \frac{1}{m s}, β (V) = η_{T} \cdot 5 \cdot \exp (\frac{V + 38}{- 18}) \cdot \frac{1}{m s} .

Here, $η_{T}$ corrects for differences between the temperature of the simulated cell $T_{sim}$ and the temperature for which the channel equations were defined $T_{eq}$ based on a temperature sensitivity $Q_{10}$ which can vary between ion channels and state transitions:

η_{T} = \exp (\log (Q_{10}) \cdot (T_{sim} - T_{eq}) / 10 K) .

There are several sources for model uncertainty about the exact channel kinetics. First, not all channel models used here were developed based on mouse data resulting in species-dependent differences. Second, we do not always know the exact subtypes of ion channels, for example in the case of the T-type calcium channel. Third, the exact temperature sensitivities $Q_{10}$ are not known. Therefore, we estimated transition rates and thresholds for state transitions during the parameter inference. For this, we allowed for offsets $Δ V$ relative to $V$ in the rate equations and additionally, we estimated relative time constants τ for the rates. For example Equation 4 was changed to:

\frac{δ}{δ t} c = (1 - c) \cdot \frac{1}{τ_{α}} \cdot α (V - Δ V_{α}) - c \cdot \frac{1}{τ_{β}} \cdot β (V - Δ V_{β}) .

Table 2

Ion channel implementation details and optimized channel parameters.

Channel	NeuronC	Type	States	Parameters	Channel remarks and references
Kainate rec.	AMPA1	MS	7	$S T C$ , $τ_{γ}$	Based on Jonas et al., 1993.
mGluR6	mGluR			$S T C$	See NeuronC documenation.
${C a}_{L}$	CA0	HH	(4)	$Δ V_{α}$ , $τ_{α}$	Based on Karschin and Lipton, 1989.
${C a}_{T}$	CA7	MS	12	$Δ V_{α}$ , $τ_{α}$	Modification of Lee et al., 2003.
${C a}_{P}$				${C a}_{P K}$	See NeuronC documenation.
${H C N}_{1 / 2 / 4}$	K4	MS	10		Based on Altomare et al., 2001.
$K_{V}$	K0	HH	(5)	$Δ V_{α}$ , $τ_{α}$	Based on Hodgkin and Huxley, 1952.
$K_{i r}$	K5	MS	3	$Δ V_{α}$	Modification of Dong and Werblin, 1995.
${C l}_{C a}$	CLCA1	MS	12		Modification of Hirschberg et al., 1998.
${N a}_{V}$	NA5	MS	9	$Δ V_{α}$ , $Δ V_{γ}$ , $τ_{all}$	Based on Clancy and Kass, 2004.

To keep the parameter space as small as possible, we only optimized the kinetics of ion channels with high uncertainty (e.g. $K_{V}$ ) or with high relevance for the exact timing of the neurotransmitter release (e.g. ${C a}_{L}$ and ${C a}_{T}$ ). Additionally, we constrained the channel parameters to physiologically plausible ranges. Table 2 summarizes which channel parameters were estimated during parameter optimization. Time constants τ and voltage offsets $Δ V$ not optimized were set to one and zero, respectively. For the ${N a}_{V}$ , a single time constant $τ_{all}$ was used to modify all time constants proportionally. The calcium pump dynamics were modified by changing the calcium concentration ${C a}_{P K}$ that causes half of the maximum calcium extrusion velocity. The BC glutamate receptors were optimized by allowing for a change in the synaptic transmitter concentration at the receptors by a factor of $S T C$ , which might be smaller for the OFF-BC than for the ON-BC given the greater distance between the release sites of the cones and the dendritic tips of the BCs (Behrens et al., 2016). The simulated cell temperature $T_{sim}$ was set to 37°C if not stated otherwise. For further information we refer to the NeuronC documentation (Smith, 1992).

Neurotransmitter release

Request a detailed protocol

The glutamate release of cones and BCs release is mediated through ribbon synapses that release vesicles in response to calcium influx in a nonlinear way (Matthews and Fuchs, 2010; tom Dieck and Brandstätter, 2006; Baden et al., 2013). We modeled the ribbon synapses with a standard model (Smith, 1992) including a readily releasable pool (RRP) from which vesicles can be released (Lagnado and Schmitz, 2015). The presence of multiple release pools shapes the dynamic of release at the ribbon synapse and make it state dependent, allowing for rapid adaptational processes at the synaptic site (Baden et al., 2013). In the model, the current release rate is dependent on the number of vesicles currently available $v_{RRP}$ in the RRP, the maximum number of vesicles $v_{RRP}^{max}$ in the RRP and the intracellular calcium concentration $[C a]$ . In NeuronC, calcium is modeled in radial shells through which calcium can diffuse deeper into the neuron. For the release of neurotransmitter, only the calcium concentration in the first shell ${[C a]}_{0}$ (equivalent to the concentration at the membrane) is considered. The release rate $r$ is computed as:

r (t) = {([C a]_{0} (t) \frac{1 e 6}{m o l})}^{2} \cdot \frac{v_{RRP} (t)}{v_{RRP}^{max}} \cdot g_{l} \cdot \frac{v e s i c l e s}{s},

where $g_{l}$ is a linear gain factor. $g_{l}$ and $v_{RRP}^{max}$ were optimized for every cell type individually. The RRP is constantly replenished with a constant rate that is equivalent to the maximum sustainable release rate $r_{m s r}$ . At a time t, for a simulation time step $Δ t$ , the vesicles in the pool are updated as follows:

v_{RRP} (t + Δ t) = v_{RRP} (t) - r (t) \cdot Δ t + r_{m s r} \cdot Δ t \cdot (1 - \frac{v_{RRP} (t)}{v_{RRP}^{max}}) .

For the cone model, $r_{m s r}$ was set to 100 vesicles per second based on Berntson and Taylor, 2003. The prior for $v_{RRP}^{max}$ was based on RRP sizes reported for salamander (Thoreson et al., 2016; Bartoletti et al., 2010). For the BCs, $r_{m s r}$ was set to eight vesicles per second based on the reported value for rat rod bipolar cells in Singer and Diamond, 2006. The prior for $v_{RRP}^{max}$ was based on Wan and Heidelberger, 2011.

Bayesian inference for model parameters

To estimate the free parameters of the multicompartment models, we used a Bayesian likelihood-free inference framework called Sequential Neural Posterior Estimation (SNPE) (Lueckmann et al., 2017; Gonçalves et al., 2020). The goal of the parameter estimation was to find parameter regions for which the model outputs match the experimentally observed glutamate release in response to a light stimulus. Details of the target data, the stimulus, the comparison between experimental and simulated data and the inference algorithm are described below. To be able to simulate the light response of the BC models, we inferred the parameters of the cone model first.

Target data of neuron models

Request a detailed protocol

As target data, we used two-photon imaging data recorded with an intensity-based glutamate-sensing fluorescent reporter (iGluSnFR) (Marvin et al., 2013). All animal procedures were approved by the governmental review board (Regierungspräsidium Tübingen, Baden-Württemberg, Konrad-Adenauer-Str. 20, 72072 Tübingen, Germany) and performed according to the laws governing animal experimentation issued by the German Government.

To constrain the cone models, we used glutamate traces of two cone axon terminals (Figure 5—figure supplement 1A) in response to a full-field chirp light stimulus (Figure 5A). The traces were recorded in one transgenic mouse (B6;129S6-Chat^tm2(cre)LowlJ, JAX 006410, crossbred with Gt(ROSA)26Sor^{tm9(CAG-tdTomato)Hze}, JAX 007905) that expressed the glutamate biosensor iGluSnFR ubiquitously across all retinal layers after intravitreal injection of the viral vector AAV2.7m8.hSyn.iGluSnFR (provided by D. Dalkara, Institut de la Vision, Paris). The cone glutamate release in the outer plexiform layer was recorded in x-z scans (64 × 56 pixels at 11.16 Hz; Zhao et al., 2019). Region-of-interest (ROIs) were drawn manually and traces of single ROIs were then normalized and upsampled to 500 Hz as described previously (Franke et al., 2017; Szatko et al., 2019). For each axon terminal, we computed the mean over five traces. Both means were then aligned by minimizing the mean squared error between them, and the mean of the two aligned means was used as target data for the cone model (Figure 5—figure supplement 1B).

For the BC models, we used mean glutamate traces of BC3a (n = 19 ROIs) and BC5o (n = 13 ROIs) (Figure 5—figure supplement 1C–F) in response to a chirp light stimulus (Figure 6A) from a recently published dataset (Franke et al., 2017). In that study, glutamate responses were recorded from BC terminals at different depths of the inner plexiform layer (x-y scans, 64 × 16 pixels at 31.25 Hz). ROIs were drawn automatically based on local image correlation and traces of single ROIs were normalized and upsampled to 500 Hz (see above). Since we simulated isolated BCs (except for the cone input), we used the responses to a local ‘chirp’ light stimulus recorded with the glycine receptor blocker strychnine, which means that the target data is less affected by inhibition from small-field amacrine cells. We did not consider input from GABAergic, wide-field amacrine cells, because these are not strongly activated by the local chirp stimulus (Franke et al., 2017). The shape of the BC stimulus differed from the cone stimulus as contrast was not linearized for the BC recordings and therefore intensity modulations below 20% brightness were weakly rectified.

Light stimulus and cell response

Request a detailed protocol

We first matched the experimental with the simulated stimulus. For this, we used the digital stimuli and corrected both timing and amplitude (using a sigmoid function) to minimize the mean squared error with respect to the experimentally recorded stimuli, correcting for delays and non-linearities in the displaying process. Then we linearly transformed the light stimulus such that the simulated photon absorption rates were $10 \times 10^{3} P^{*} / (s \cdot c o n e)$ for the lowest and $31 \times 10^{3} P^{*} / (s \cdot c o n e)$ for the highest stimulus intensity including the background illumination, approximating the values reported in Franke et al., 2017. In NeuronC, the photon absorption rate acts as input to a phototransduction model (Nikonov et al., 1998), which provides the hyperpolarizing current entering the inner segment. The membrane potential in the axon terminal compartment regulates the calcium influx into the cell which in turn influences the glutamate release rate. This glutamate release from the simulated cones modifies the opening probability (the fraction of open channels in the deterministic case) of postsynaptic receptors, which drive the BC models.

Discrepancy function

Request a detailed protocol

To compare model outputs to the experimentally observed target data, we defined a discrepancy function $𝜹$ . Since the target traces were relative fluorescence intensities, the absolute number of released glutamate vesicles could not be directly inferred from the target data, and the data only constrained relative variations in the release rate during simulation. Because we also wanted to constrain our models to plausible membrane potentials and release rates, we combined the following seven discrepancy measures:

$δ_{iGluSnFR}$ : The mean squared error between the experimental and simulated iGluSnFR trace.
$δ_{Rate}^{Rest}$ : A penalty for implausibly high resting release rates.
$δ_{V}^{Rest}$ : A penalty for implausibly low or implausible high resting membrane potentials.
$δ_{Rate}^{Δ}$ : A penalty for implausibly low release rate changes.
$δ_{V}^{Δ}$ : A penalty for implausibly low membrane potential changes.
$δ_{V}^{m i n}$ : A penalty for implausibly low membrane potentials.
$δ_{V}^{m a x}$ : A penalty for implausibly high membrane potentials.

The discrepancy between a model output m and the target data $𝝂_{t}$ was computed as:

𝜹 (𝝂_{t}, m) = {[δ_{iGluSnFR} (𝝂_{t}, m), δ_{Rate}^{Rest} (m), δ_{V}^{Rest} (m), δ_{Rate}^{Δ} (m), δ_{V}^{Δ} (m), δ_{V}^{m i n} (m), δ_{V}^{m a x} (m)]}^{⊤} .

To identify the overall ‘best’ samples, we computed the total discrepancy as the absolute-value norm of the discrepancy vector: $δ_{tot} (𝝂_{t}, m) = ∥ δ (𝝂_{t}, m) ∥$ .

The discrepancy function $δ_{iGluSnFR}$ (Equation 10) computes the distance between a simulated iGluSnFR trace $𝝂_{m}$ and an iGluSnFR target $𝝂_{t}$ . To estimate the simulated iGluSnFR trace, we convolved the glutamate release rate $𝐫_{m}$ with an iGluSnFR kernel $κ$ . Here, the time-dependent kernel function $κ$ was approximated with an exponential decay function, based on iGluSnFR intensity changes to spontaneous vesicle release reported in Marvin et al., 2013:

ν_{m} = r_{m} * κ, κ (t) = \exp (- t / 60 m s) .

The discrepancy was then computed as the euclidean distance between the simulated and the target iGluSnFR trace with respect to a distance minimizing linear transformation of the simulated trace. This linear transformation was necessary because the target traces only reflect relative fluorescence changes. The discrepancy was normalized to be between zero and one by dividing by the variance ${∥ 𝝂_{t} - μ_{t} ∥}^{2}$ , where $μ_{t}$ is the mean of the target data.

δ_{iGluSnFR} (𝝂_{t}, 𝝂_{m}) = \min_{a, b} \frac{{∥ 𝝂_{t} - (a + b \cdot 𝝂_{m}) ∥}^{2}}{{∥ 𝝂_{t} - μ_{t} ∥}^{2}}, b \geq 0 .

For all other discrepancies, specific values of the glutamate release rate (in the case of the BCs, the mean release rate over all synapses) or the somatic membrane potential were compared to a lower and an upper bound of target values $t_{l}$ and $t_{u}$ , such that values within these bounds were assigned a discrepancy of 0.0. Outside this range, the discrepancy was defined by additional bounds $p_{l}$ and $p_{u}$ . Given a specific value of the simulation $y_{m}$ , the respective discrepancy $δ_{∙}^{∙} (y_{m})$ was computed as:

δ_{∙}^{∙} (y_{m}) = {\begin{cases} - 1 + \exp (- 2 \frac{(y_{m} - t_{l})^{2}}{(t_{l} - p_{l})^{2}}), & i f y_{m} \in (- \infty, t_{l}), \\ 0, & i f y_{m} \in [t_{l}, t_{u}], \\ 1 - \exp (- 2 \frac{(y_{m} - t_{u})^{2}}{(t_{u} - p_{u})^{2}}), & i f y_{m} \in (t_{u}, \infty) . \end{cases}

To compute $δ_{Rate}^{Rest}$ and $δ_{V}^{Rest}$ , the resting release rate $r_{m}^{0}$ and resting membrane potential $v_{m}^{0}$ for the background light adapted state were extracted. For the BC models, the resting membrane potential was not penalized for values between $t_{l} = - 65 m V$ and $t_{u} = - 45 m V$ based on reported values for mice (Ichinose et al., 2014; Ichinose and Hellmer, 2016) and rat retina (Ma et al., 2005). For the cone model, the expected resting membrane potential was more depolarized between $t_{l} = - 55 m V$ and $t_{u} = - 40 m V$ (Cangiano et al., 2012).

The discrepancy of the resting release rate $δ_{Rate}^{Rest}$ was computed similarly. For the BC models, the lower bound t_l was set to zero. As mentioned earlier, we limited our BC models to have a maximum sustainable release rate of 8 vesicles per second based on Singer and Diamond, 2006. We allowed non-zero resting release rates due to the background light and spontaneous vesicle fusion but constrained it to values lower than the maximum sustainable release rate (Kavalali, 2015; Baden et al., 2014). For the OFF-BC we chose an upper bound of 4 vesicles per second (half the maximum sustainable release rate). For the ON-BC, we chose a slightly smaller value of 3 vesicles per second. This difference was based on the observation that the ON-BC target never falls significantly below the value of the resting state, indicating that the resting release rate is probably close to zero and can therefore not become smaller. In contrast, the OFF-BC target falls below the resting value right after stimulus onset, indicating a small but non-zero resting release rate. For the cone model, we assumed a comparably high resting release rate between $t_{l} = 50$ and $t_{u} = 80$ vesicles per second based on the assumed higher maximum sustainable release rate and the fact that cones show steady release in darkness (Choi et al., 2005; Sheng et al., 2007).

For the penalty on implausible release changes $δ_{Rate}^{Δ}$ , we computed the largest absolute difference $Δ r$ between the resting release rate $r_{m}^{0}$ and release rates $𝐫_{m}$ after stimulus onset. $δ_{V}^{Δ}$ was computed analogously but for the membrane potential $𝐯_{m}$ and the resting membrane potential $v_{m}^{0}$ :

Δ r = max | r_{m} - r_{m}^{0} | a n d Δ v = max | v_{m} - v_{m}^{0} | .

$δ_{Rate}^{Δ} (y_{m})$ and $δ_{V}^{Δ} (y_{m})$ were then computed by using the differences $y_{m} = Δ r$ and $y_{m} = Δ v$ , respectively, in Equation 13. For the BC release rate, we did not penalize differences larger than $t_{l} = 5$ vesicles per second. For the cone, we expected much larger differences between $t_{l} = 50$ to $t_{u} = 65$ vesicles per second due to their larger maximum sustainable release rate. For the membrane potential, we expected a difference of at least $t_{l} = 5 m V$ based on light step responses recorded with patch clamp in mouse BCs (Ichinose et al., 2014; Ichinose and Hellmer, 2016). Since here, the stimulus contrast was higher, we only used the reported values as lower bounds but allowed the model to have larger variation, namely up to $t_{u} = 25 m V$ for the OFF- and $t_{u} = 15 m V$ for the ON-BC. We allowed greater membrane potential variation in the OFF-BC, because it receives input from more cones.

For the discrepancy measures $δ_{V}^{m i n}$ and $δ_{V}^{m a x}$ , we computed the minimum and maximum of the membrane potential $𝐯_{𝐦}$ after stimulus onset and used again Equation 13. For $δ_{V}^{m i n}$ , we chose $t_{l} = - 80 m V$ for the BCs and $t_{l} = - 60 m V$ for the cone model, and in both cases $t_{u} = \infty$ . For $δ_{V}^{m a x}$ , we chose $t_{u} = - 10 m V$ for the BCs and $t_{u} = - 35 m V$ for the cone, and in both cases we set $t_{l} = - \infty$ . The BC values are based on data from rat (Ma et al., 2005) and ground squirrel (Saszik and DeVries, 2012); the cone values are based on Cangiano et al., 2012.

All values for $p_{l}$ and $p_{u}$ were based on pilot simulations with the goal to distribute the penalties where they most mattered. All discrepancies (except for $δ_{iGluSnFR}$ ) and their respective values $p_{l}$ , $p_{u}$ , $t_{l}$ and $t_{u}$ are illustrated in Figure 2 for clarity and summarized in Table 3.

Figure 2

Download asset Open asset

Discrepancy measures based on Equation 13 for the cone (red dashed line), the OFF- (blue solid line) and ON-BC (green dotted line).

The parameters defining the discrepancy measures are listed in Table 3. All discrepancy measures are between zero and one per definition.

Table 3

Parameters of discrepancy measures.

$δ_{∙}^{∙}$	Cone				BC (3a \| 5o)				References
	$p_{l}$	$t_{l}$	$t_{u}$	$p_{u}$	$p_{l}$	$t_{l}$	$t_{u}$	$p_{u}$
$δ_{Rate}^{Rest}$	0	50	80	100	0	0	3 \| 4	7	Choi et al., 2005; Sheng et al., 2007; Berntson and Taylor, 2003; Singer and Diamond, 2006
$δ_{V}^{Rest}$	−80	−55	−40	−20	−80	−65	−45	−30	Cangiano et al., 2012; Ichinose et al., 2014; Ichinose and Hellmer, 2016; Ma et al., 2005
$δ_{Rate}^{Δ}$	0	50	65	100	0	5	∞	∞
$δ_{V}^{Δ}$	0	5	10	20	0	5	15 \| 25	40	Ichinose et al., 2014; Ichinose and Hellmer, 2016
$δ_{V}^{m i n}$	−75	60	∞	∞	−100	−80	∞	∞	Cangiano et al., 2012; Ma et al., 2005; Saszik and DeVries, 2012
$δ_{V}^{m a x}$	$- \infty$	$- \infty$	−35	−20	$- \infty$	$- \infty$	−10	0	Cangiano et al., 2012; Ma et al., 2005; Saszik and DeVries, 2012

Priors

The inference method SNPE is a Bayesian method and therefore it needs a prior distribution $p (𝜽)$ for the parameters $𝜽$ to estimate the posterior. We chose truncated normal distributions for all priors because they allow for weighting of more plausible parameters (in contrast to e.g. uniform distributions), while they enable restrictions to plausible ranges (in contrast to e.g. normal distributions). A $d$ -dimensional truncated normal distribution $𝒩_{T}$ is defined by a mean $𝝁 = {(μ_{1}, \dots, μ_{d})}^{T}$ , a $d \times d$ covariance matrix $Σ$ and a $d$ -dimensional space $W = [a_{1}, b_{1}] \times \dots \times [a_{d}, b_{d}]$ :

𝒩_{T} (θ | μ, Σ, W) = {\begin{cases} \frac{\exp (- 0.5 (θ - μ)^{T} Σ^{- 1} (θ - μ))}{\int_{W} \exp (- 0.5 (ω - μ)^{T} Σ^{- 1} (ω - μ)) d ω} & if θ \in W, \\ 0 & otherwise . \end{cases}

The prior means $μ_{i}$ and truncation bounds $[a_{i}, b_{i}]$ were based on experimental data wherever possible (see Appendix 1—table 1 and 2), including data from rat and different cell types such as rod bipolar cells, as well as pilot simulations. For parameter inference, we normalized the parameter space such that the truncation bounds were [0, 1] in all dimensions. The diagonal entries of the prior covariance matrix $Σ$ were set to 0.3². Because it is difficult to find prior knowledge about the dependencies of parameters, we set all non-diagonal entries to zero. To sample from $𝒩_{T}$ , we implemented a rejection sampler, that samples from a normal distribution with the same mean $𝝁$ and covariance matrix $Σ$ and resamples all $𝜽$ not in $W$ .

Inference algorithm

Request a detailed protocol

SNPE estimates a posterior parameter distribution represented by a mixture-density network, based on sampling, that is, model evaluations for randomly drawn parameters. Inference is performed in several rounds. In every round $j$ , the algorithm draws $N$ parameters from a sampling distribution ${\tilde{p}}_{j} (𝜽)$ to estimate the posterior distribution $p (𝜽 | 𝐱_{target})$ , where $𝐱_{target}$ is a summary statistic of the target data.

In the first round, parameter samples $𝜽_{n}$ are drawn from the prior, that is, ${\tilde{p}}_{1} (𝜽) = p (𝜽)$ , and the multicompartment model is evaluated for all $𝜽_{n}$ . From each simulated response, a summary statistic $𝐱_{n}$ is computed, resulting in $N$ pairs of parameters and summary statistics $(𝜽_{n}, 𝐱_{n})$ . At the end of the round, a mixture-density network is trained with summary statistics $𝐱$ as input, and the parameters $ϕ$ of a mixture of Gaussian distribution $q_{ϕ} (𝜽 | 𝐱)$ as output. The network is trained by minimizing the loss function $ℒ$ :

ℒ (ϕ) = - \frac{1}{N} \sum_{n}^{N} \frac{p (θ_{n})}{\tilde{p} (θ_{n})} K (x_{n}) \log (q_{ϕ} (θ_{n} | x_{n})),

where $K$ is a kernel function with values between zero and one that weights the influence of samples on the network training. $K (𝐱_{n})$ is close to one for samples with summary statistics $𝐱_{n}$ close to the the target summary statistic $𝐱_{target}$ and becomes smaller with increasing discrepancy between $𝐱_{n}$ and $𝐱_{target}$ . This means, the network tries to find parameter distributions $q_{ϕ} (𝜽 | 𝐱)$ that describe the distribution of samples for any given summary statistic $𝐱$ . Or, in other words, the network is trained to find a mapping from summary statistics to parameter distributions. $K$ ensures that the network focuses its capacity on summary statistics close to the target summary statistic. After training the network, it can be evaluated at a summary statistic $𝐱^{*}$ to obtain the posterior parameter distributions for the given summary statistic. Evaluating at $𝐱^{*} = 𝐱_{target}$ yields an approximation of the true posterior distribution $p (𝜽 | 𝐱_{target}) \approx q_{ϕ} (𝜽 | 𝐱_{target})$ . This posterior can either be used as the sampling distribution of the next round ${\tilde{p}}_{j + 1} (𝜽)$ , or—if the algorithm is stopped—as the final posterior distribution. The relative probability $p (𝜽_{n}) / \tilde{p} (𝜽_{n})$ in Equation 16 weights samples not drawn from the prior, which ensures that Bayes’s rule is not violated. A detailed proof that this actually yields an approximation of the true posterior in the Bayesian sense can be found in Lueckmann et al., 2017.

We based our algorithm on the Python code available at https://github.com/mackelab/delfi version 0.5.1 (Macke Lab, 2020) with the following settings and modifications. We modeled $q_{ϕ}$ as a single Gaussian, because we noticed that mixture of Gaussians almost always collapsed to a single component after a few rounds. Both, intermediate and final posteriors were truncated using the truncation bounds of the prior. The truncation was performed after network training. For every neuron model, we drew $N = 2, 000$ samples per round and stopped the algorithm after the fourth round. Two hundred additional samples were drawn from the posterior for further analysis. Since we wanted to use the posterior samples to simulate the effects of electrical stimulation on the BCs, the number of compartments was increased in this last step to 139 and 152 for the OFF- and ON-BC, respectively.

As summary statistics of model outputs $m_{n}$ , we used the discrepancy function $δ (𝝂_{t}, m_{n}) = 𝐱_{n}$ (see Equation 10), which describes the discrepancy between model outputs and the target data. The target summary statistic was set to be a zero-vector $𝐱_{target} = {[0, \dots, 0]}^{⊤}$ , since the target should have a discrepancy of zero with respect to itself. The first dimension of $𝜹$ , $δ_{iGluSnFR}$ , computes the distance between the simulated and experimentally observed iGluSnFR trace. Considering the noise in the target data, observing a discrepancy of zero is virtually impossible. Therefore, evaluating the network at $𝐱_{target} = {[0, \dots, 0]}^{⊤}$ is based on extrapolation, that is, the mixture-density network is evaluated for an input where it was not trained on. This, as we observed during pilot experiments, often led to posterior estimates of poor quality or endless loops of resampling. So instead of evaluating the network at $𝐱_{target} = {[0, \dots, 0]}^{⊤}$ to obtain the posterior estimate, the network was instead evaluated at $𝐱_{target} = {[x_{min}^{iGluSnFR}, 0, \dots, 0]}^{⊤}$ , where $x_{min}^{iGluSnFR}$ is the the smallest $δ_{iGluSnFR}$ sampled during this round. This is roughly equivalent to assuming that the best strategy for extrapolation is to simply use the estimate at the boundary. For the weighting function $K$ , we used zero-centered Gaussian kernels $k$ with a bandwidth of $σ = 0.25$ in all dimensions but the first one. In the first dimension, that is the weighting kernel for $δ_{iGluSnFR}$ , we also used an adaptive strategy and both, the mean $μ_{iGluSnFR}$ and the bandwidth $σ_{iGluSnFR}$ of the kernel, were updated in every round:

μ_{iGluSnFR} = x_{min}^{iGluSnFR}, σ_{iGluSnFR} = q_{20}^{iGluSnFR} - x_{min}^{iGluSnFR},

where $q_{20}^{iGluSnFR}$ is the 20th percentile of all sampled iGluSnFR discrepancies of the same round. $K$ was computed as the product of all scalar kernels $k$ .

Some parameter combinations caused the neuron simulation to become numerically unstable. If a simulation could not successfully terminate for this reason, the sample was ignored during training of the mixture-density network by setting the kernel weight to zero. In other cases, the BC models had a second, strongly depolarized and therefore biologically implausible equilibrium state. To test for this, we simulated a somatic voltage clamp to 30 mV for 100 ms and checked whether the membrane potential would recover to a value of −30 mV or lower within additional 300 ms. Samples not recovering to ≤ -30 mV were also ignored during training.

Data analysis of simulated traces

Request a detailed protocol

The distance function $δ_{iGluSnFR} (𝝂_{t}, 𝝂_{m})$ (see Equation 12) was used not only to compute the discrepancy between simulations and the respective targets but also more generally to compare different experimental and simulated iGluSnFR traces. The distance between two iGluSnFR traces $𝝂_{1}$ and $𝝂_{2}$ was computed as $δ_{iGluSnFR} (𝝂_{1}, 𝝂_{2})$ .

To quantify the timing precision of our neuron models, we estimated peak times in simulated and target iGluSnFR traces to compute pairwise peak time differences. For every peak in the simulated trace, we computed the time difference to the closest peak of the same polarity (positive or negative) in the target. We did not consider peaks between 16 s and 23 s of the stimulation for the cone and between 16 s and 21 s for the BC models, because the targets were to noisy for precise peak detection in these time windows. This resulted in approximately 35 positive and negative peak time differences per trace.

Simulation of electrical stimulation

To simulate external electrical stimulation of our BC models, we implemented a two-step procedure. In the first step, the electrical field is estimated as a function of space and time across the whole retina for a given stimulation current. By setting a position of the BC multicompartment models within the retina, the extracellular voltage for every compartment can be extracted. In the second step, the extracellular voltages are applied to the respective compartments (Figure 1C) to simulate the neural response in NeuronC. To be able to perform the first step, we estimated the electrical properties of retinal tissue first. For this, we utilized the same algorithm that was used for parameter inference of the neuron models. To validate the framework, we simulated the electrical stimulation in Corna et al., 2018 and compared experimental and simulated neural responses. Finally, we utilized the framework to find electrical stimuli for selective stimulation of OFF- and ON-bipolar cells. Details of the implementation and the experimental data are described in the following.

Computing the extracellular voltage

Request a detailed protocol

We estimated the electrical field in the retina for a given electrical stimulus with the finite-element method using the software COMSOL Multiphysics (Comsol, 2019). We modeled the photoreceptor degenerated retina as a cylinder with a radius of 2 mm and a height of 105 µm (Pennesi et al., 2012). The stimulation electrodes were modeled as flat disks on the bottom of the retina. Above the retina, an additional cylinder with the same radius and a height of 2 mm was placed to model the electrolyte. The top of this cylinder was assumed to be the return electrode. The implementation of such a model with the subdivision into finite elements is shown in Figure 3. For a single circular stimulation electrode, the model was radially symmetric and could therefore be reduced to a half cross-section as shown in Figure 3 to increase the simulation speed without altering the results. The following initial and boundary conditions were applied to the model. The initial voltage was set to zero at every point $V (x, y, z, t = 0) = 0$ . The surface normal current density $j_{stim}^{⟂}$ of stimulation electrodes was always spatially homogeneous and dependent on the total stimulation current $i_{stim}$ and the total surface area of all electrodes $A_{electrode}$ :

j_{stim}^{⟂} = \frac{i_{stim}}{A_{electrode}} .

Figure 3

Download asset Open asset

Model for the external electrical stimulation of the retina.

(A) Schematic figure of the experimental setup for subretinal stimulation of ex vivo retina combined with epiretinal recording of retinal ganglion cells. Schematic modified from Corna et al., 2018. (**B,C**) Model for simulating the electrical field potential in the retina in 3D and 2D, respectively. The retina (darker blue) and the electrolyte above (lighter blue) are modeled as cylinders. The shown 3D model is radially symmetric with respect to the central axis (red dashed line). Therefore, the 3D and 2D implementations are equivalent, except that the computational costs for the 2D model are much lower. The 2D implementation is annotated with parameters that were either taken from the literature or inferred from experimental data. (D) Electrical field potential in the retina for a constant stimulation current of 0.5 µA for a single stimulation electrode with a diameter of 30 µm. Additionally, the compartments (black circles with white filling) of the ON-BC model are shown. The stimulation is subretinal meaning that the dendrites are facing the electrode (horizontal black line on bottom).

The potential of the return electrode was kept constant $V_{return} (t) = 0$ . At all other boundaries, the model was assumed to be a perfect insulator $j_{other}^{⟂} = 0$ . We assumed a spatially and temporally homogeneous conductivity and permittivity in both the retina and the electrolyte. The conductivity of the electrolyte was set to $σ_{ames} = 1.54 S/m$ based on Eickenscheidt and Zeck, 2014 and its relative permittivity was assumed to be $ε_{ames} = 78$ , based on the value for water. The conductivity $σ_{retina}$ and relative permittivity $ε_{retina}$ of the retina were optimized with respect to experimental target data as described below.

Target data to infer the electrical parameters of the retina

Request a detailed protocol

To estimate the electrical properties of the retina, we first recorded target data. All procedures were approved by the governmental review board (Regierungspräsidium Tübingen, Baden-Württemberg, Konrad-Adenauer-Str. 20, 72072 Tübingen, Germany, AZ 35/9185.82–7) and performed according to the laws governing animal experimentation issued by the German Government. We applied different sinusoidal stimulation voltages $v_{stim}$ and recorded the evoked currents. Currents were recorded with ( $i_{retina}^{rec}$ ) and without ( $i_{ames}^{rec}$ ) retinal tissue placed on the micro-electrode array. In both cases, the recording chamber was filled with an electrolyte (Ames’ medium, A 1420, Sigma, Germany). A single Ag/AgCl pellet (E201ML, Science Products) was used as a reference electrode and located approximately 1 cm above a customized micro-electrode array. The electrodes, made of sputtered iridium oxide had diameters of 30 µm and center-to-center distance of 70 µm. The stimulation current was calculated from the voltage drop across a serial 10 resistor in series with the Ag/AgCl electrode (Corna et al., 2018). The voltage drop was amplified using a commercial voltage amplifier (DLPVA, Femto Messtechnik GmbH, Berlin, Germany) and recorded using the analog input (ME 2100, Multi Channel Systems MCS GmnH, Germany). Stimulation currents were measured across an ex vivo retina of a rd10 mouse (female; post-natal day 114; strain: Pde6b^rd10 JAX Stock No: 004297).

We applied sinusoidal voltages of 25 and 40 Hz. For 25 Hz, we applied amplitudes from 100 to 600 mV with steps of 100 mV. For 40 Hz all amplitudes were halved.

Procedure to infer the electrical parameters of the retina

Request a detailed protocol

We estimated the conductivity $σ_{retina}$ and relative permittivity $ε_{retina}$ of the retina in three steps based on the experimental voltages $v_{stim}$ and the respective recorded currents $i_{retina}^{rec}$ and $i_{ames}^{rec}$ . To facilitate the following steps, we fitted sinusoids $i_{retina}$ and $i_{ames}$ to the slightly skewed recorded currents and used them in the following (Figure 8C). To fit the sinusoids, we minimized the mean squared error between recorded currents and idealized sinusoidal currents of the same frequency f, resulting in estimates of the phase $ϕ (i_{ames})$ and the amplitude $A (i_{ames})$ of the currents:

ϕ (x), A (x) = \underset{ϕ, A}{a r g m i n} \int_{t} (x - A \cdot \sin (2 π f t + ϕ))^{2} d t .

During parameter inference, we only used two voltage amplitudes per frequency, resulting in four voltage and eight current traces. The other amplitudes were used for model validation. First, we estimated the electrical properties of the electrode. Here, ‘electrode’ is meant to include the electrical double layer and all parasitic resistances and capacitances in the electrical circuit. We simulated the voltage $v_{ames}$ across the electrolyte without retinal tissue by applying the currents $i_{ames}$ as stimuli (Figure 8Ai). Since this setup does not contain anything besides the electrolyte and the electrode, the difference between the experimental stimulus $v_{stim}$ , which was applied to record $i_{ames}$ , and the simulated voltage $v_{ames}$ was assumed to have dropped over the electrode:

v_{electrode} = v_{stim} - v_{ames} .

Based on that assumption, we could estimate the electrical properties of the electrode. We modeled the electrode as a RC parallel circuit (Figure 8Aii). Having both, sinusoidal voltages ( $v_{electrode}$ ) over and the respective sinusoidal currents ( $i_{ames}$ ) through the electrode, we analytically computed the values for $R_{e}$ and $C_{e}$ as follows. We assumed $R_{e}$ and $C_{e}$ to be dependent on $v_{stim}$ and therefore to be dependent on the stimulus frequency and amplitude. From the data we derived the phase $ϕ_{Z}$ and amplitude $| Z |$ of the impedance formed by the RC circuit. For every $v_{electrode}$ , we estimated $ϕ (v_{electrode})$ and $A (v_{electrode})$ using Equation 19. $ϕ_{Z}$ and $| Z |$ were then computed as:

ϕ_{Z} = ϕ (v_{electrode}) - ϕ (i_{ames}), | Z | = A (v_{electrode}) / A (i_{ames}) .

Then, knowing the frequency f, $ϕ_{Z}$ and $| Z |$ are sufficient to compute $R_{e}$ and $C_{e}$ :

R_{e} = | Z | \sqrt{1 + t a n {(ϕ_{Z})}^{2}}, C_{e} = - t a n (ϕ_{Z}) / (2 π f R_{e}) .

With the estimated values of the RC circuit, we created a model with only two unknowns, the conductivity $σ_{retina}$ and the relative permittivity $ε_{retina}$ of the retina (Figure 8Aiii). To estimate the unknown parameters of this model, we used the same inference algorithm as for the neuron models but with a different discrepancy function. Here, the discrepancy $δ_{R} (v_{stim})$ for a stimulus $v_{stim}$ was computed as the mean squared error between the respective experimental current (now with retinal tissue) $i_{retina}$ and the simulated current $i_{retina}^{sim}$ :

δ_{R} (v_{stim}) = \sum_{v_{stim}} \int_{t} {(i_{retina} - i_{retina}^{sim})}^{2} 𝑑 t .

The total discrepancy was computed as the sum of all discrepancies $δ_{R} (v_{stim})$ for the four different $v_{stim}$ stimuli that were used. To cover a wider range of possible parameters, we first estimated the parameters in a logarithmic space by sampling the exponents $p_{σ}$ and $p_{ε}$ of the parameters:

σ_{retina} = 2^{p_{σ}} \cdot 0.1 S / m, ϵ_{retina} = 2^{p_{ε}} \cdot 10^{6} .

We used normal distributions (without truncation) as priors for $p_{σ}$ and $p_{ε}$ and set the means to 1.0 and the standard deviations to 2.0. After three rounds with 50 samples each, we computed the minimum ( $a_{σ}$ , $a_{ε}$ ), maximum ( $b_{σ}$ , $b_{ε}$ ) and mean ( $μ_{σ}$ , $μ_{ε}$ ) for both parameters $σ_{retina}$ and $ε_{retina}$ from the 10% best samples. Then, we then ran the parameter inference algorithm again, but now in a linear parameter space around the best samples observed in the logarithmic space. For the priors of $σ_{retina}$ and $ε_{retina}$ , we used truncated normal priors bound to $[a_{σ}, b_{σ}]$ and $[a_{ε}, b_{ε}]$ with means $μ_{σ}$ and $μ_{ε}$ , respectively. As for the cell parameter inference, we normalized the parameter space to values between in [0, 1]. The diagonal entries of the prior covariance matrix were set to 0.3², with non-diagonal entries of zero. The parameters resulting in the lowest sampled discrepancy during optimization are referred to as the optimized parameters and were used to simulate the neural responses to electrical stimulation.

Simulation of the neural response to electrical stimulation

Request a detailed protocol

With the optimized parameters for the electrical properties of the retina, we were able to compute the BC responses for any given stimulation current. Note that for this, we used the model illustrated in Figure 3 as described earlier but with the optimized parameters for $σ_{retina}$ and $ε_{retina}$ . To simulate the neural response, we first used the stimulation current to simulate the extracellular voltage over time within the retina. After defining the relative position of the multicompartment model with respect to the retinal cylinder, we extracted the extracellular voltage for each compartment at its the central position (Figure 3C). Finally, these extracellular voltages were applied to the compartment models in NeuronC to simulate their response (Figure 1C). To estimate the uncertainty of the BC responses to electrical stimulation, we simulated different cell parametrizations in every stimulation setting. For this, we used the five best posterior samples, that is, the five (out of 200) samples with the smallest $δ_{tot}$ , for both BC models. In all simulations, we modeled subretinal stimulation of photoreceptor degenerated retina (Zrenner, 2002). For this, we removed all cone input from the BCs and virtually placed the multicompartment models in the retinal cylinder such that the dendrites were facing towards the electrode. The z-position of BC somata, that is, the distance to the bottom of the retinal cylinder, was set to 30 µm.

Model validation

Request a detailed protocol

To validate the model for electrical stimulation, we compared simulated BC responses to experimentally recorded retinal ganglion cell (RGC) thresholds to 4 ms biphasic current pulses reported in Corna et al., 2018. In this study, the RGC thresholds were recorded epiretinally under subretinal stimulation of photoreceptor degenerated (rd10) mouse retina using a micro-electrode array (Figure 3A). The stimulation threshold was defined as the charge delivered during the anodic stimulation phase evoking 50% of the firing rate of a specific RGC. On the micro-electrode array. The 30 µm diameter electrodes were arranged on a regular grid with a center-center spacing of 70 µm. The RGC thresholds were measured for different numbers $N$ of $N$ × $N$ active electrodes.

We simulated the electrical field in the retina for the configurations with 1×1, 2×2, 4×4 and 10×10 active electrodes using the respective currents from the experimental data. The electrodes were centered with respect to the retina. For every stimulation current, we simulated the response of the OFF- and ON-BC at six xy-positions with distances from 0 to 500 µm relative to the center. Simulation temperature $T_{sim}$ was set to 33.5°C to match experimental conditions. For every 40 ms simulation, we computed the mean number of vesicles released per synapse.

Optimizing electrical stimulation to separately activate ON- and OFF-BCs

Request a detailed protocol

To find stimuli for selective stimulation of ON- and OFF-BCs, we simulated the response of the BC models to different electrical stimuli. For this, we used a single 30 µm diameter electrode and centered the dendrites of the simulated BCs above this electrode. To find stimuli that stimulate the OFF-BC without stimulating the ON-BC or vice versa, we utilized the same algorithm used for estimating the BC parameters. Here, the inference algorithm was used to estimate parameters of a 40 ms stimulation current $i_{stim}$ parametrized by four free parameters $p_{1}, \dots, p_{4}$ . The current was defined as a cubic spline fit through the knot vector $𝐚 = (0, p_{1}, \dots, p_{4}, p^{*}, 0)$ spaced equidistantly in time between zero and 40 ms, where $p^{*}$ is chosen such that the stimulus is charge neutral (i.e. the integral over the current is zero). For all stimuli, the maximum stimulus amplitude was normalized to 0.5 µA. An illustration is shown in Figure 10.

Here, the priors over $p_{1}, \dots, p_{4}$ were Gaussian with zero means and standard deviations of 0.3. For every sampled stimulus $i_{stim}^{n}$ , we simulated the response of the BCs for $Δ t = 60 m s$ starting with the stimulus onset. For parameter estimation, we defined the discrepancy measure $δ_{stim}^{t}$ as the ratio between the relative release $R_{n}$ of the OFF- and ON-BC which was defined as:

R_{n} = \frac{(μ (𝐫_{n}) - μ (𝐫_{base})) Δ t}{v_{RRP}^{max} + (r_{m s r} - μ (𝐫_{base})) Δ t},

where $μ (𝐫_{n})$ is the evoked mean release and $μ (𝐫_{base})$ is the base release rate in the absence of electrical stimulation; that is, the numerator is equal to the number of released vesicles (as a mean over all synapses) caused by the stimulation. The denominator is equal to the theoretical maximum of releasable vesicles per synapse (see Equation 9). $δ_{stim}^{t}$ was computed as:

δ_{stim}^{t} (i_{stim}^{n}) = R_{n} (o t h e r B C) / R_{n} (t a r g e t B C) .

We ran the parameter inference twice (each with one round only), once with the ON- and once with the OFF-BC as target. We drew 400 samples from the prior that were reused for the second run of inference, and 100 more samples from the posterior. Here, the posteriors were two-component mixture of Gaussians without truncation.

Code and data availability

Request a detailed protocol

Models and simulation code is available at https://github.com/berenslab/CBC_inference (Oesterle, 2020; copy archived at swh:1:rev:2b8ec4ac0ca916d42cba0404229298f8ff79c3a3). Experimental and inference data is available at https://zenodo.org/record/4185955.

Results

We used a high-resolution electronmicroscopy data set (Helmstaedter et al., 2013) to create biophysically realistic multicompartment models of three neuron types from the mouse retina including cones, an OFF- and an ON-bipolar cell (BC) type. These neurons form the very beginning of the visual pathway, with cones converting light into electrochemical signals and providing input via sign-preserving and -reversing synapses to OFF- and ON-BCs, respectively. The parameters of these models include the basic electrical properties of the cells as well as the density of different ion channel types in different compartments. Given a set of parameters, simulations from the model can easily be generated; however, it is not possible to evaluate the likelihood for a given set of parameters, which would be required for standard Bayesian inference procedures for example through MCMC.

To overcome the challenge of choosing the resulting 20 to 32 parameters of these models, we adapted a recently developed technique called Sequential Neural Posterior Estimation (SNPE) (Lueckmann et al., 2017) (for details, see Materials and methods). Starting from prior assumptions about the parameters, the algorithm compared the result of simulations from the model to data obtained by two-photon imaging of the glutamate release from the neurons (Franke et al., 2017) and measured a discrepancy value between the simulation and the data. Based on this information, the algorithm used a neural network to iteratively find a distribution over parameters consistent with the measured data. This yielded optimized biophysically realistic models for the considered neuron types.

Inference of cone parameters

We first estimated the posterior distribution over the parameters of a cone based on the glutamate release of a cone stimulated with a full-field chirp light stimulus, consisting of changes in brightness, temporal frequency, and temporal contrast (Figure 4A and Figure 5). The cone model had a simplified morphology and consisted of four compartments (Figure 1, see Materials and methods). We included a number of ion channels in the model reported to exist in the cones of mice or closely related species (see Table 1). Prior distributions were chosen based on the literature. For inference, we drew 2000 samples of different parameter settings per round and stopped the algorithm after the fourth round. Then, 200 more parameter samples were drawn from the respective posteriors for further analysis. The chosen discrepancy functions penalized discrepancies between the target and simulated iGluSnFR trace $δ_{iGluSnFR}$ , implausible membrane potentials, and implausible release rates. To compare different model fits, the discrepancy measures were added to yield a total discrepancy $δ_{tot}$ . We found that the total discrepancy $δ_{tot}$ of the cone model was relatively high for most samples drawn from the prior but decreased over four rounds of sampling (Figure 4A). The discrepancy measuring the fit quality to the glutamate recording $δ_{iGluSnFR}$ was already relatively small in the first round for most, but not all samples. In the following rounds, the number of samples with large $δ_{iGluSnFR}$ was strongly reduced (Figure 4A).

Figure 4

Download asset Open asset

Discrepancies of samples from the cone and the BC models during and after parameter estimation.

(**A, B, C**) Sampled discrepancies for the cone (A), the OFF- (B), and ON-BC (C) respectively. For every model, the total discrepancy $δ_{tot}$ (left) and the discrepancy between the simulated and target iGluSnFR trace $δ_{iGluSnFR}$ (right) are shown. For every model, four optimization rounds with 2000 samples each were drawn (indicated by gray vertical lines). After the last round (indicated by dashed vertical lines, ‘p.’), 200 more samples were drawn from the posteriors. For the BCs, the number of compartments was increased in this last step to 139 and 152 for the OFF- and ON-BC, respectively. Additionally, 200 samples were drawn from assuming independent posterior marginals for comparison (indicated by blue vertical lines, ‘m.’). For every round, the discrepancy distribution (horizontal histograms), the median discrepancies (red vertical lines), the 25th to 75th percentile (blue shaded area) and the 5th to 95th percentile (gray-shaded area) are shown. (D) Discrepancies between different iGluSnFR traces of BCs to demonstrate the high precision of the model fit. The pairwise discrepancy computed with equation Equation 12 between eight iGluSnFR traces is depicted in a heat map. The column and row labels indicate which $𝝂_{t}$ and $𝝂_{m}$ were used in equation Equation 12 respectively. The traces consists of the optimized BC models (‘Model’), the targets used during optimization (‘Target’), experimental data from the same cell type without the application of any drug (‘No drug’) and experimental data from another retinal CBC type with the application of strychnine (‘BC4’ and ‘BC7’). Note that strychnine was also applied to record the targets.

Figure 4—source data 1 Sample discrepancies of all samples shown in (A-C).: https://cdn.elifesciences.org/articles/54997/elife-54997-fig4-data1-v2.zip
Download elife-54997-fig4-data1-v2.zip
Figure 4—source data 2 Discrepancies shown in (D), and the respective (mean) iGluSnFR traces.: https://cdn.elifesciences.org/articles/54997/elife-54997-fig4-data2-v2.zip
Download elife-54997-fig4-data2-v2.zip

Figure 5 with 2 supplements see all

Download asset Open asset

Optimized cone model.

(A) Normalized light stimulus. (B) Somatic membrane potential relative to the resting potential for the best parameters (blue line) and for 200 samples from the posterior shown as the median (gray dashed line) and 10th to 90th percentile (gray shaded area). A histogram over all resting potentials is shown on the right. (C) Release rate relative to the resting rate. Otherwise as in (B). (D) Simulated iGluSnFR trace (as in (B)) compared to target trace (orange). Three regions (indicated by gray dashed boxes) are shown in more detail below without samples from the posterior. Estimates of positive and negative peaks are highlighted (up- and downward facing triangles respectively) in the target (brown) and in the simulated trace (blue and cyan respectively). Pairwise time differences between target and simulated peaks (indicated by triangle pairs connected by a black line) are shown as histograms for positive (blue) and negative (cyan) peaks on the right. The median over all peak time differences is shown as a black vertical line.

Figure 5—source data 1 iGluSnFR traces used for constraining the cone and BC models.: https://cdn.elifesciences.org/articles/54997/elife-54997-fig5-data1-v2.zip
Download elife-54997-fig5-data1-v2.zip
Figure 5—source data 2 Stimulus, target, and cell responses, including responses with removed ion channels.: https://cdn.elifesciences.org/articles/54997/elife-54997-fig5-data2-v2.zip
Download elife-54997-fig5-data2-v2.zip

The parameter setting with lowest discrepancy ( $δ_{tot} = 0.10$ ) modeled accurately the response of the cone to full-field stimulation with the chirp light stimulus (Figure 5). The simulated iGluSnFR signal nicely matched the data both on a coarse timescale and in the millisecond regime (Figure 5D). Indeed, for this sample, all discrepancies besides $δ_{iGluSnFR}$ were zero or almost zero ( $δ_{tot} - δ_{iGluSnFR} < 0.0001$ ) and most of the remaining discrepancy could be attributed to the noisy target data.

As our inference algorithm returned not only a single best set of parameters, but also a posterior distribution, we could obtain additional parameter samples from the model which should produce simulations consistent with the data. Almost all samples from the posterior yielded simulations that matched the target data well (median $δ_{iGluSnFR}$ : 0.12) and the overall total discrepancy was small (median $δ_{tot}$ : 0.21). Besides the discrepancy between the experimental and simulated glutamate trace $δ_{iGluSnFR}$ , most of the remaining discrepancy in the posterior samples was caused by rate variation (mean $∥ δ_{Rate}^{Δ} ∥$ : 0.18) and resting rates (mean $∥ δ_{Rate}^{Rest} ∥$ : 0.14) that were too low in some of the models. While in principle we could propagate the remaining uncertainty about the model parameters provided by the posterior to the inference about BC models, we used only the parameter set with the smallest total discrepancy $δ_{tot}$ for efficiency and refer to this as the optimized cone model. To analyze the role of active ion channels, we removed ion channels individually (except for the ${C a}_{L}$ channel with is necessary to simulate the vesicle release) from the optimized cone and simulated the light response (Figure 5—figure supplement 2). We found that the HCN channel contributed most, while the contribution of ${C l}_{C a}$ was negligible. Since ${C l}_{C a}$ did not alter the light response for both, the cone and BC light stimulus, we removed it in the following steps for computational efficiency.

Inference of bipolar cell parameters

We next turned to anatomically detailed multicompartment models of two BC types. We chose to model type 3a and type 5o because these were the OFF- and ON-CBC types for which we could gather most prior information from the literature. The anatomy of the cells was extracted from a 3D reconstruction of examples of these cell types based on electron microscopy data (Helmstaedter et al., 2013) and divided into regions sharing ion channel parameters (Figure 1). As for the cone model, the channels included in the model and the prior distributions were chosen based on the literature (see Table 1). This yielded 32 and 27 free parameters for the OFF- and ON-BC, respectively.

We fitted the BC type models to published two-photon glutamate imaging data reflecting the glutamate release from the BC axon terminals (Franke et al., 2017). In this case, we used responses to a local chirp light stimulus activating largely the excitatory center of the cells. In addition, the responses were measured in the presence of the drug strychnine to block locally acting inhibitory feedback through small-field amacrine cells (Franke et al., 2017) (see Materials and methods for details). Similar to what we observed in cones, the total discrepancy $δ_{tot}$ for parameter sets sampled for the OFF- and ON-BC model decreased over the four rounds of optimization (Figure 4B and C). In contrast to the the cone model, the discrepancy measure penalizing deviations from the glutamate trace $δ_{iGluSnFR}$ was relatively large for prior samples and declined approximately equally fast as the total discrepancy $δ_{tot}$ .

We found that simulations generated with the parameter set with minimal total discrepancy or parameters sampled from the posterior matched the target traces very well for both OFF- and ON-BC models (Figure 6). For these parameters, the cells were relatively ispotential units throughout the light stimulus (Figure 6—figure supplement 1 and Figure 6—figure supplement 2) with a larger voltage gradients from dendrites to the axon in the ON-BC. The optimized BC models, that is the BC samples with the lowest total discrepancies $δ_{tot}$ , had discrepancies of zero except for the iGluSnFR discrepancy $δ_{iGluSnFR}$ . To get a more quantitative impression of the quality of the model fits, we compared the pairwise iGluSnFR discrepancies $δ_{iGluSnFR}$ between the optimized BC models, the experimentally measured response traces as used during optimization, traces recorded from the same cell type without application of strychnine and responses of another OFF- and ON-BC. For both optimized cell model outputs, the discrepancy was smallest for the targets used during optimization. This shows that the optimized models were able to reproduce cell-type specific details in light response properties that go beyond the simple distinction of ON and OFF responses. While the discrepancies between traces of different ON-BC types were overall relatively small for local stimulation (Franke et al., 2017), the discrepancies between traces from OFF cells were larger likely due to network modulation of the target cell type by amacrine cells (indicated by the difference between the target and the no-drug condition) and larger response differences between the two compared OFF-BC types. The posterior samples of both BC models had a low discrepancy, except for a few samples (median $δ_{tot}$ : 0.29 and 0.26 of the OFF- and ON-BC, respectively). The only discrepancy measure with a non-zero median of the absolute values was $δ_{iGluSnFR}$ , which also accounts for 88% and 82% of the mean total discrepancy for the OFF- and ON-BC respectively.

Figure 6 with 6 supplements see all

Download asset Open asset

Optimized BC models.

(A) Normalized light stimulus. Responses of the OFF- and ON-BC are shown in (**B–D**) and (**E–G**), respectively. (**B, E**) Somatic membrane potential relative to the resting potential for the best parameters (blue line) and for 200 samples from the posterior shown as the median (gray dashed line) and 10th to 90th percentile (gray-shaded area). A histogram over all resting potentials is shown on the right. (**C, F**) Mean release rate over all synapses relative to the mean resting rate. Otherwise as in (B). (**D, G**) Simulated iGluSnFR trace (as in (B)) compared to respective target trace (orange). Three regions (indicated by gray dashed boxes) are shown in more detail below without samples from the posterior. Estimates of positive and negative peaks are highlighted (up- and downward facing triangles, respectively) in the target (brown) and in the simulated trace (blue and cyan, respectively). Pairwise time differences between target and simulated peaks (indicated by triangle pairs connected by a black line) are shown as histograms for positive (blue) and negative (cyan) peaks on the right. The median over all peak time differences is shown as a black vertical line.

Figure 6—source data 1 Stimulus, target, and cell responses, including responses with removed ion channels.: https://cdn.elifesciences.org/articles/54997/elife-54997-fig6-data1-v2.zip
Download elife-54997-fig6-data1-v2.zip

Despite the overall high resemblance between optimized model outputs and targets, there were some visible systematic differences. For the ON-BC, the target showed a skewed sinusoidal response with faster rise than fall times during the frequency and amplitude increasing sinusoidal light stimulation between 10 s and 18 s and between 20 s and 27 s respectively. In contrast, the optimized model output showed approximately equal rise and fall times, resulting in a systematic delay of positive and negative peaks (median delay of all peaks: 15.6 ms) in the simulated iGluSnFR trace relative to the target (Figure 6G). Additionally, some of the positive peaks of the optimized ON-BC model during sinusoidal light stimulation were too small (e.g. at 11.5 s). This effect might have been a side-effect of the peak timing difference between target and model: Amplitude increases were inefficient in reducing the discrepancy as long as the peaks were not precisely aligned. In contrast, the peak time precision of the OFF-BC model (Figure 6D) was much higher (median delay of all peaks: 0.0 ms). In this case, the main difficulty for the model appeared to be its inability to reproduce the non-linearity in the cell response to the increasing amplitude sinusoidal light stimulation between 20 s and 27 s.

After having verified that the posterior over parameters provided a good fit to the experimental data, we inspected the one-dimensional marginal distributions to learn more about the resulting cellular models (Figure 7). For most parameters, the marginal posteriors had smaller variances than the priors, indicating that the parameter bounds were not chosen too narrow. For some parameters, the posterior mean differed substantially from the prior mean (e.g. the ${C a}_{T}$ channel density at the axon terminal of OFF-BC) while it was largely unchanged for others (e.g. the ${C a}_{L}$ channel density at the soma for the OFF-BC). The algorithm also inferred the dependencies of some parameters, visible in the two-dimensional marginals (Figure 7—figure supplement 1 and Figure 7—figure supplement 2). Because of these correlations, the full posterior in the high-dimensional parameter space led to simulations which were on average better (median: 0.29 vs. 0.31 and 0.26 vs. 0.33 for the OFF- and ON-BC, respectively) and less variable in their quality (95%-CIs: 0.53 vs. 1.01 and 0.64 vs. 1.42 for the OFF- and ON-BC, respectively) than parameters drawn from a posterior obtained by assuming independent marginal distributions. In most cases, the parameters resulting in the lowest total discrepancy were close to the means of the respective posteriors. For some parameters there was a strong difference between the marginal posteriors of the OFF- and ON-BC. For example, the two parameters controlling the leak conductance, $V_{r}$ and $R_{m}$ , were much lower for the OFF-BC consistent with the strong variation of membrane resistances reported in Oltedal et al., 2009. The membrane conductance was lower for the ON-BC, which could increase signal transduction speed in the longer axon. Even though the posteriors were narrower than the priors, they still covered a wide range of different parameters. To some extent, this may reflect the fact that we fit the model parameters solely on the cells output, and for example dendritic parameters may be underconstrained by such data; in addition, it may also reflect variability between cells of the same type seen in the experimental data that has also been reported in other studies (Franke et al., 2017).

Figure 7 with 2 supplements see all

Download asset Open asset

Parameter distributions of the BC models.

1D-marginal prior (dashed gray line) and posterior distributions (solid lines) are shown for the OFF- (blue) and ON-BC (green). The parameters of the posterior samples with the lowest total discrepancy are shown as dashed vertical lines in the respective color. $X_{Y} @ Z$ refers to the channel density of channel $X_{Y}$ at location $Z$ . Locations are abbreviated; S: soma, A: axon, D: dendrite and AT: axon terminals (see Figure 1 and main text for details). Note that although these 1D-marginal distributions seem relatively wide in some cases, the full high-dimensional posterior has much more structure than a posterior distribution obtained from assuming independent marginals (see Figure 4). Not all parameter distributions are shown.

Figure 7—source data 1 Prior and posterior parameters for the OFF- and ON-BC.: https://cdn.elifesciences.org/articles/54997/elife-54997-fig7-data1-v2.zip
Download elife-54997-fig7-data1-v2.zip

After the fourth optimization round, 200 samples were drawn from the posterior distribution with an increased number of compartments to find model parameters to simulate electrical stimulation (see Methods). For comparison, we also ran simulations with the same parameters but the original number of compartments (data not shown). Interestingly, more than 85% of the samples had a lower discrepancy if the models were simulated with more compartments for both BCs. For the best 20% (i.e. 40 samples) of the posterior samples (sorted with respect to samples with fewer compartments), the samples with more compartments had lower discrepancies with only one exception per cell. This indicates that, given enough computational power, the same parameter inference approach but with more compartments could further improve the model outputs. From these samples, we used the five samples with the smallest total discrepancies $δ_{tot}$ for the simulation of electrical stimulation.

Additionally, we used these five samples to analyze the effect of active ion conductances on the light response by removing individual ion channels types from the BCs (Figure 6—figure supplement 3 and Figure 6—figure supplement 4). Similar to the optimized cone model, the HCN channels played the most important role in shaping the light response. For both cells, the ${N a}_{V}$ and somatic calcium channels barely had any influence on the membrane voltage or the vesicle release rate.

Simulating electrical stimulation of the retina

To provide an exemplary use-case for our data-driven biophysical models of retinal neurons, we asked whether we could use our simulation framework to optimize the stimulation strategy for retinal neuroprosthetic implants. These implants have been developed for patients who lose sight due to degeneration of their photoreceptors (Zrenner, 2002). While existing implants have been reported to provide basic vision (Zrenner, 2002; Edwards et al., 2018; Luo and da Cruz, 2016), they are far from perfect. For example, most current stimulation strategies likely activate OFF- and ON-pathways at the same time (Barriga-Rivera et al., 2017). To this end, we created a simulation framework for subretinal electrical stimulation of retinal tissue with micro-electrode arrays. We estimated the conductivity and relative permittivity of the retina based on experimentally measured currents evoked by sinusoidal voltages and then validated simulations of the electrical stimulation of our fitted BC models with standard pulse like stimuli against responses measured in RGCs (Corna et al., 2018). Finally, we used the simulation framework to find stimuli that can separately stimulate OFF- and ON-BCs.

Our framework for simulating the effect of external electrical stimulation using the inferred BC models consisted of two steps: we first estimated the electrical field resulting from a stimulation protocol as a function of space and time across the whole retina (Figure 3). The corresponding extracellular voltages were then applied to the respective compartments to simulate the neural response. To do so, we needed a model of the electrical properties of the electrode and the retinal tissue. We assumed disk electrodes and a simplified model assuming homogeneous electrical properties of the retina and the surrounding medium (see Materials and methods). This model contained two free parameters that needed to be estimated from data: the conductivity and relative permittivity of the retinal tissue.

To estimate these parameters we recorded electrical currents resulting from sinusoidal voltage stimulation with different frequencies in a dish with and without a retina (Figure 8B,C). We used the data without a retina to estimate the properties of the stimulation electrode (Figure 8A,D,E and Materials and methods). Based on the estimates of the electrode properties and the data recorded with a retina, we estimated the conductivity and relative permittivity of the retina with the same parameter inference method as for the neuron models.

Figure 8

Download asset Open asset

Estimation of the conductivity $σ_{retina}$ and the relative permittivity $ε_{retina}$ of the retina for simulating electrical stimulation.

(A) Electrical circuits used during parameter estimation. (B) Stimulation voltages $v_{stim}$ at 25 (left) and 40 Hz (right). From the six experimentally applied amplitudes, only the amplitudes used for parameter inference are shown. (C) Experimentally measured currents $i_{ames}$ through electrolyte (Ames’ solution) without retinal tissue for the stimulus voltages $v_{stim}$ in (B). The mean traces over all (but the first two) repetitions are shown (black dashed lines). Sine waves were fitted to the mean traces (solid lines) with colors referring to the voltages in (B). (D) Simulated voltages over the electrolyte $v_{ames}$ using the fitted currents in (C) as stimuli applied to the circuit in (Ai). (Aii) Electrical circuit used to model the electrode plus interface. (E) Stimulus frequency and amplitude dependent estimates of $R_{e}$ (i) and $C_{e}$ (ii) based on the electrical circuit shown in (Aii) for 25 (black) and 40 Hz (gray). Note that the values were derived analytically (see main text). The values corresponding to the stimulus voltages shown in (B) are highlighted with respective color. (Aiii) Electrical circuit used to estimate $σ_{retina}$ and $ε_{retina}$ . The respective values for $R_{e}$ and $C_{e}$ are shown in (E) and are dependent on $v_{stim}$ . The current $i_{retina}$ through the model is measured for a given stimulus voltage $v_{stim}$ . (F) Sampled parameters of $ε_{retina}$ and $σ_{retina}$ and the respective sample losses. First, samples were drawn in a wide logarithmic space (gray dots) and then in a narrower linear parameter space. The best sample (lowest discrepancy) is highlighted in red. (G) Simulated currents $i_{retina}$ (solid lines) through the circuit in (Aiii) with optimized parameters (red dot in (F)) and respective experimentally measured currents (broken lines). Here, results for all six stimulus amplitudes are shown for both frequencies.

Figure 8—source data 1 Currents, voltages, and samples from optimization.: https://cdn.elifesciences.org/articles/54997/elife-54997-fig8-data1-v2.zip
Download elife-54997-fig8-data1-v2.zip

We found that both parameters are very well constrained by the measured data (Figure 8F). The parameters resulting in the lowest discrepancy were $σ_{r e t i n a} = 0.076 S / m$ and $ε_{r e t i n a} = 1.1 \times 10^{7}$ in accordance with the conductivity of $0.077 S / m$ reported for rabbit Eickenscheidt et al., 2012 and the relative permittivity of gray matter estimated in Gabriel et al., 1996. With these parameters, we simulated currents for all stimulus amplitudes we recorded experimentally. The simulated and experimental currents matched for the amplitudes used during parameter inference but also for amplitudes reserved for model validation (Figure 8G). Therefore, we used them in all following experiments.

To validate our simulation framework, we compared simulated BC responses to experimentally measured RGC thresholds (Corna et al., 2018). We simulated BCs at different positions for four different electrode configurations (Figure 9A) and nine stimulation current wave forms (Figure 9B). For small stimulus charge densities, the BCs barely responded, whereas for very high charge densities the cells released all glutamate vesicles available in the readily release pool (Figure 9C and D). In between, the response increased from no response to saturation, dependent on the distances of the simulated cell to the active electrodes. Across stimulation conditions, these threshold regions coincide with the measured RGC thresholds to the same stimulation, when assuming that the stimulated RGCs were not too far away from the stimulation electrode. Since the reported RGC thresholds likely result from indirect stimulation via BCs, the consistency between the RGC and simulated BC thresholds can be interpreted as evidence that our model was well calibrated to simulate electrical stimulation.

Figure 9

Download asset Open asset

Threshold of electrical stimulation for experimentally measured RGCs and simulated BCs of photoreceptor degenerated mouse retina.

(A) Stimulation currents measured experimentally and used as stimuli in the simulations. (B) xy-positions of BCs (crosses) and electrodes (red dots) for 1×1, 2×2, 4×4, 10×10 stimulation electrodes, respectively. Every electrode is modeled as a disc with a 30 µm diameter. Except for the electrode configuration, the models were as in Figure 8. (**C, D**) Mean synaptic glutamate release relative to the size of the readily releasable pool $v_{RRP}^{max}$ of simulated OFF- and ON-BCs, respectively. Values can be greater than 1, because the pool is replenished during simulation. Errorbars indicate the standard deviation between simulations of BCs with different parametrizations; for both BCs, the five best posterior samples were simulated. Glutamate release is shown for different charge densities (x-axis) and cell positions (colors correspond to xy-positions in (A); for example the darkest blue corresponds to the central BC). Experimentally measured RGC thresholds (gray dashed lines) plus-minus one standard deviation (gray-shaded ares) are shown in the same plots.

Figure 9—source data 1 Current traces and vesicle release for all simulations.: https://cdn.elifesciences.org/articles/54997/elife-54997-fig9-data1-v2.zip
Download elife-54997-fig9-data1-v2.zip

Optimized electrical stimulation for selective OFF- and ON-BC stimulation

We finally used our framework for electrical stimulation to find stimuli that excite OFF- or ON-BCs selectively. To this end, we performed Bayesian inference over an electrical charge-neutral stimulus (Figure 10A) with the SNPE algorithm, using the response ratio between the two BC types (Figure 10B) as a discrepancy function to minimize. Using this procedure, we found that triphasic, anodic first stimuli with a cathodic middle phase of high amplitude (Figure 10C) evoked substantial neurotransmitter release in the OFF-BC (Figure 10Di) while evoking almost no response in the ON-BC (Figure 10Dii). The stimuli optimized to target the ON-BC were biphasic, with no clear preference for anodic or cathodic first currents (Figure 10E). In contrast to the stimuli optimized for the OFF-BC, these stimuli did not exclusively stimulate the targeted ON-BC (Figure 10Fii) but also the OFF-BC (Figure 10Fi). We did not find stimuli evoking stronger release (defined as in Equation 24) in the ON-BC than the OFF-BC. This lower threshold of the OFF-BC, which we also observed for the biphasic current pulses (Figure 9), was partially caused by the shorter axon of the OFF-BC resulting in slightly larger changes of the extracellular voltage at the axon terminals during stimulation (Figure 10—figure supplement 1A,B). While the morphologies and ion channel distributions contributed to differences in the membrane voltage at the axon terminals (Figure 10—figure supplement 1C), the decisive factor for differences in the OFF and ON response presumably lies in the presence or absence of ${C a}_{T}$ channels at the axon terminals. Removing these channels from the OFF-BC had almost no effect on the membrane voltage, but resulted in qualitatively similar responses for both cell types (Figure 10—figure supplement 1D,E). Removing the ${C a}_{L}$ channels at the axon terminals on the other hand, left the OFF response for the triphasic and biphasic cathodic first stimuli largely unchanged (Figure 10—figure supplement 1F).

Figure 10 with 1 supplement see all

Download asset Open asset

Electrical currents optimized for selective ON- and OFF-BC stimulation.

(A) Illustration of the stimulus parametrization. A stimulation current was computed by fitting a cubic spline through points defined by $p_{1} \dots p_{4}$ and $p^{*}$ (black dots) that were placed equidistantly in time. $p_{1} \dots p_{4}$ were free parameters and $p^{*}$ was chosen such that the stimulus was charge neutral, that is the integral over the current (gray-shaded area) was zero. Currents were normalized such that the maximum amplitude was always 0.5 µA. (B) Illustration of the electrical stimulation of an OFF- (i) and ON-BC (ii) multicompartment model. (**C, E**) Stimuli optimized for selective OFF- and ON-BC stimulation, respectively. The two best stimuli observed during optimization are shown. (**D, F**) Cumulative vesicle release (as a mean over all synapses) relative to the size of the readily releasable pool $v_{RRP}^{max}$ in response to the electrical stimuli shown in (C) and (E), respectively. Stimuli and responses are shown in matching colors. Release is shown for the five best posterior cell parametrizations for the OFF- (i) and ON-BC (ii) as a mean (line) plus-minus one standard deviation (shaded area).

Figure 10—source data 1 Optimized stimuli and BC responses, including responses with removed ion channels.: https://cdn.elifesciences.org/articles/54997/elife-54997-fig10-data1-v2.zip
Download elife-54997-fig10-data1-v2.zip

Discussion

In this study, we showed how the recently developed Bayesian likelihood-free parameter inference method called Sequential Neural Posterior Estimation (SNPE) (Lueckmann et al., 2017) can be used to estimate the parameters of multicompartment models of retinal neurons based on light-response measurements. In addition, we built a model for electrical stimulation of the retina, and optimized electrical stimulation protocols for retinal neuroprosthetic devices designed to support vision in the blind.

Performing Bayesian inference for mechanistic models is difficult, as they typically do not come with an analytical solution of the likelihood function. The SNPE algorithm — like many approximate Bayesian computation (ABC) methods (Sisson et al., 2018) — does not require such an analytical formulation, as it builds on simulations of the model. In contrast to standard ABC methods, the SNPE algorithm makes efficient use of simulation time by using all simulated data to train a mixture density network to update the parameter distributions (Lueckmann et al., 2017; Gonçalves et al., 2020). Moreover, SNPE makes minimal assumptions about the simulator, giving full flexibility to use it with different simulation environments and software. As all Bayesian methods, SNPE allows the incorporation of literature knowledge in the form of priors which can be used to constrain parameters and to put more weight on more plausible parameters. Finally, it does not only yield a point-estimate of the best parameters — like exhaustive grid-search techniques (Goldman et al., 2001; Prinz et al., 2003; Stringer et al., 2016) or evolutionary algorithms (Gerken et al., 2006; Keren et al., 2005; Achard and De Schutter, 2006; Rossant et al., 2011) — but also returns a posterior distribution that reflects remaining parameter uncertainties and allows one to detect dependencies between parameters.

Recently, there has been a surge of interest in Bayesian simulator-based inference with many recently published algorithms (Gutmann and Corander, 2016; Papamakarios and Murray, 2016; Lueckmann et al., 2017; Lintusaari et al., 2017; Papamakarios et al., 2018; Wood, 2010; Durkan et al., 2018; Sisson et al., 2018; Gonçalves et al., 2020; Bittner et al., 2019). While we initially evaluated different algorithms, we did not perform a systematic comparison or benchmarking effort, which is beyond the scope of this project. Much of the literature on simulator-based inference evaluates the proposed algorithms on fairly simple models. In contrast, we used SNPE here to perform parameter inference of comparatively complex multicompartment models of neurons. Importantly, we did not need targeted experiments to constrain these models, but based our framework on two-photon imaging data of glutamate release in response to simple light stimuli using a genetically encoded indicator called iGluSnFR (Marvin et al., 2013; Franke et al., 2017). This methods allows direct experimental access to the neurotransmitter release of all excitatory neurons of the retina (Euler et al., 2019). Using this data, we inferred the distributions of model parameters relevant for all the intermediate steps between light stimulation of cones to the glutamate release from synaptic ribbons. While we optimized many parameters in the models using SNPE, we chose to keep some parameters on sensible default values to avoid issues with computational complexity. Of course, it is possible that optimization of the full parameter space would have yielded slightly better results or that some parameters would have been set to slightly different values, as a mechanism whose parameter was allowed to vary compensated for the one whose parameter was held fixed. As an alternative to our approach, one can combine classical systems identification approaches with inference for only some of the biophysical mechanisms such as the ribbon synapse (Schröder et al., 2019). Our approach, however, allows the exploration of mechanisms within neurons which are not or only barely experimentally accessible. For example, in BCs, it is currently difficult to obtain direct voltage measurements from any part of the cell but the soma. If one is interested in how the electrical signal propagates through the axon or the axon terminal, our simulations may help to obtain mechanistic insights and develop causal hypotheses.

Because our inference framework works with experimental measurements which can be performed in a high-throughput manner, it allows for a comparably easy scaling to infer model parameters of a larger number of multicompartment models e.g. of different neuron types. In principle it could even be possible to infer the parameters of a neuron by imaging another neuron. For example, one could attempt to infer parameters of an amacrine cell by observing the neurotransmitter release of a connected BC — although such an indirect inference would likely result in larger uncertainties. Ideally, such a large-scale approach would also include realistic morphologies, for example from electron microscopy as shown here. In fact, BCs are anatomically relatively ‘simple’ neurons, and it will be interesting to test our inference methods on neurons with more complicated structure such as some amacrine cells (Masland, 2012). While we did not aim at an exhaustive analysis of the effect of morphology on the neuron responses, one could easily explore how details of the morphology influence the distribution of optimal biophysical parameters. Further, we extended our model to simulate and optimize external electrical stimulation of the retina. For blind patients suffering from photoreceptor degeneration, for example caused by diseases like Retinitis Pigmentosa, neuroprosthetic devices electrically stimulating remaining retinal neurons can restore basic visual abilities (Edwards et al., 2018; Luo and da Cruz, 2016). The spatial and temporal resolution of such retinal implants is, however, still very limited (Weitz et al., 2015) with the highest reported restored visual acuity of 20/546 (Chuang et al., 2014). While many experimental studies have explored different strategies of stimulation (Jensen et al., 2005; Jensen and Rizzo, 2008; Tsai et al., 2009; Eickenscheidt et al., 2012), most of them are restricted to very specific stimulus types such as current or voltages pulses. As a consequence, retinal implants are not able to specifically target cell types such as the independent stimulation of the ON and OFF pathways of the retina (Lee and Im, 2019; Barriga-Rivera et al., 2017; Twyford et al., 2014). To facilitate a systematic stimulus optimization in silico, we developed a simulation framework for electrical stimulation of the retina. While the idea to simulate responses of BCs to electrical stimuli is not new, previous studies restricted their models to point/ball source electrodes (Resatz and Rattay, 2004; Rattay et al., 2017), simplified BCs to passive cables (Gerhardt et al., 2010) or used simplified BC models that only express L- or T-type channels (Werginz et al., 2015). Our framework combines the simulation of micro-electrode arrays used in neuroprosthetic devices (Edwards et al., 2018; Luo and da Cruz, 2016) with detailed models of an OFF- and ON-BC. This allowed us to test a large number of stimulus waveforms, optimizing for stimuli selectively targeting either OFF- or ON-BCs, which could help to better encode visual scenes into electrical signals of retinal implants. We found stimuli selectively targeting the OFF-BC without stimulating the ON-BC, but not vice versa. Likely, the main reason for the differential response of the two BC types was that only the OFF-BC had T-type calcium channels at the axon terminals. These channels were more sensitive to transient changes in the membrane voltage which were evoked by the stimuli optimized to selectively target the OFF-BC. The ON-BC, having no T-type calcium channels and an overall higher threshold, did not respond to these stimuli. However, it could be stimulated with longer anodic stimulus phases activating the L-type calcium channels. Since we modeled the cells in isolation, network effects through synaptic activation of amacrine cells might further modulate the activity of the BCs. However, the neurites and somata of amacrine cells are substantially farther away from the stimulation electrode than those of the BCs, and these effects might be comparably small. That notwithstanding, simulations including network effects and also more diverse BC types will be required in the future. Ideally, in silico optimized stimulation strategies would be first verified in ex vivo experiments before implementing them in neuroprosthetic devices to improve the quality of visual impressions.

To be able to simulate the electrical stimulation of the retina, we first inferred the conductivity and relative permittivity of the rd10 retina based on recorded currents evoked by sinusoidal stimulation voltages. While the estimated conductivity ( $σ_{r e t i n a} = 0.076 S / m$ ) is consistent with the value ( $σ_{r e t i n a} = 0.077 S / m$ ) reported in Eickenscheidt et al., 2012, also smaller ( $0.025 S / m$ , Karwoski and Xu, 1999) and larger (≈ $0.75 S / m$ , Wang and Weiland, 2015) conductivities have been reported for the retina. These differences may be due to different ways in tissue handling and preparation. Comparing the estimated values of the relative permittivity ( $ε_{r e t i n a} = 1.1 \times 10^{7}$ ) to the literature is more difficult, and most simulation studies choose to ignore its effects. The relative permittivity of retinal tissue has been reported for very high frequencies (10 MHz to 10 GHz), but the strong frequency dependence makes a direct comparison to frequencies several orders of magnitude smaller (e.g. 40 Hz) not meaningful. Additionally, data from gray matter suggest a relative permittivity of 1.5 × 10⁷ at 40 Hz very close to our estimate (Gabriel et al., 1996). This opens the question weather the common assumption to ignore the effects of the relative permittivity in other simulations (Gerhardt et al., 2010; Werginz et al., 2015; Rattay et al., 2017) is valid.

In summary, mechanistic models in neuroscience such as biophysically realistic multicompartment models have long been regarded as requiring many manual parameter choices or highly specific experiments to constrain them. We showed here how relatively standard, high-throughput imaging data in combination with likelihood-free inference techniques can be used to perform Bayesian inference on such models, allowing unprecedented possibilities for efficient optimization and analysis of such models. Importantly, this allow us to understand which parameters in such models are well constrained, and which are not, and determine which parameter combinations lead to similar simulation outcomes (Gonçalves et al., 2020; Alonso and Marder, 2019) — questions that have hindered progress in the field for years.

Appendix 1

Cone prior parameters

Appendix 1—table 1

Cone model parameter priors.

Parameter	Unit	a_i	µ_i	b_i	References: direct and additional
$C_{m}$	$μ F {c m}^{- 2}$	0.9	1	1.3	Oltedal et al., 2009
$R_{m}$	$k Ω {c m}^{2}$	1	5	100	Oltedal et al., 2009
$V_{r}$	$mV$	-90	-70	-50
${C a}_{L}$	$m S {c m}^{- 2}$	0.1	2	10	Morgans et al., 2005; Mansergh et al., 2005; Cui et al., 2012
$K_{V}$	$m S {c m}^{- 2}$	0	0.1	3	Van Hook et al., 2019
$H C N @ S$	$m S {c m}^{- 2}$	0	3	10	Knop et al., 2008; Hellmer et al., 2016; Van Hook et al., 2019
$H C N @ A$	$m S {c m}^{- 2}$	0	1	10	Knop et al., 2008; Hellmer et al., 2016; Van Hook et al., 2019
$H C N @ A T$	$m S {c m}^{- 2}$	0	0.1	10	Knop et al., 2008; Hellmer et al., 2016; Van Hook et al., 2019
${C l}_{C a}$	$m S {c m}^{- 2}$	0	0.1	10	Yang et al., 2008; Caputo et al., 2015
${C a}_{P}$	$μ S {c m}^{- 2}$	0.1	10	100	Morgans et al., 1998
$τ_{α} ({C a}_{L})$		0.75	1	1.5
$τ_{α} (K_{V})$		0.1	1	10
$Δ V_{α} ({C a}_{L})$	$mV$	−5	0	5
$Δ V_{α} (K_{V})$	$mV$	−10	0	10
${C a}_{P K}$	$μ m o l$	0.01	5	100
$v_{RRP}^{\max}$	$v e s i c l e s$	10	20	30	Thoreson et al., 2016; Bartoletti et al., 2010
$g_{l}$		0.3	1	3

BC prior parameters

Appendix 1—table 2

BC model parameter priors.

Parameter	Unit	a_i (3a\|5)	µ_i (3a\|5)	b_i (3a \| 5)	References: direct and additional
$C_{m}$	$μ F {c m}^{- 2}$	0.9	1.18	1.3	Oltedal et al., 2009
$R_{m}$	$k Ω {c m}^{2}$	1	26	100	Oltedal et al., 2009
$V_{r}$	$mV$	−90	−70	−50
${C a}_{L} @ S$	$m S {c m}^{- 2}$	0	0.5	3	Cui et al., 2012, Van Hook et al., 2019
${C a}_{L} @ A T$	$m S {c m}^{- 2}$	0.1	0.5	3	Cui et al., 2012, Van Hook et al., 2019
${C a}_{T} @ S$	$m S {c m}^{- 2}$	0\|n/a	0.5\|n/a	3\|n/a	Cui et al., 2012, Van Hook et al., 2019; Hu et al., 2009; Satoh et al., 1998
${C a}_{T} @ A T$	$m S {c m}^{- 2}$	0\|n/a	0.5\|n/a	3\|n/a	Cui et al., 2012, Van Hook et al., 2019; Hu et al., 2009; Satoh et al., 1998
$K_{V} @ D$	$m S {c m}^{- 2}$	0	0.4	2	Ma et al., 2005
$K_{V} @ P A$	$m S {c m}^{- 2}$	0	0.4	2	Ma et al., 2005
$K_{V} @ A$	$m S {c m}^{- 2}$	0	0.4	2	Ma et al., 2005
$K_{i r} @ S$	$m S {c m}^{- 2}$	0	1	2	Cui and Pan, 2008, Knop et al., 2008
$H C N @ D$	$m S {c m}^{- 2}$	0	0.2	2	Hellmer et al., 2016, Knop et al., 2008
$H C N @ S$	$m S {c m}^{- 2}$	0	0.2	2	Hellmer et al., 2016, Knop et al., 2008
$H C N @ A T$	$m S {c m}^{- 2}$	0	0.2	2	Hellmer et al., 2016, Knop et al., 2008
${N a}_{V} @ D A$	$m S {c m}^{- 2}$	0	20	100	Hellmer et al., 2016
${C a}_{P} @ S$	$μ S {c m}^{- 2}$	0.1	10	100	Morgans et al., 1998
${C a}_{P} @ A T$	$μ S {c m}^{- 2}$	0.1	10	100	Morgans et al., 1998
$τ_{γ} (K a i n a t e)$		1\|n/a	5\|n/a	20\|n/a
$τ_{α} ({C a}_{L})$		0.5	1	2
$τ_{α} ({C a}_{T})$		0.5\|n/a	1\|n/a	2\|n/a
$τ_{α} (K_{V})$		0.1	1	10
$τ_{a l l} ({N a}_{V})$		0.5	1	2
$Δ V_{α} ({C a}_{L})$	$mV$	−10	0	10
$Δ V_{α} ({C a}_{T})$	$mV$	−10\|n/a	0\|n/a	10\|n/a
$Δ V_{α} (K_{V})$	$mV$	−10	0	10
$Δ V_{α} (K_{i r})$	$mV$	−5	0	5
$Δ V_{α} ({N a}_{V})$	$mV$	−5	0	5
$Δ V_{γ} ({N a}_{V})$	$mV$	−5	0	5
${Ca}_{PK}$	$μ m o l$	0.1\|0.01	20	100
$S T C$	$mmol$	0.05\|1	0.5\|1.5	1\|3
$v_{RRP}^{max}$	$v e s i c l e s$	4	8	15	Singer and Diamond, 2006
$g_{l}$		0.5	1	3

NeuronC parameters

Appendix 1—table 3

Other NeuronC parameters.

Parameter	Unit	Value	Remarks
vna	$mV$	65	Reversal potential sodium
vk	$mV$	−89	Reversal potential potassium
vcl	$mV$	−70	Reversal potential chloride
dnao	$mmol$	151.5	Extracellular sodium concentration
dko	$mmol$	2.5	Extracellular potassium concentration
dclo	$mmol$	133.5	Extracellular chloride concentration
dcao	$mmol$	2	Extracellular calcium concentration
dicafrac		1	Fraction of calcium pump current that is added to total current
use_ghki		1	Use Goldman-Hodgkin-Katz equation
cone_timec		0.2	Time constant of cone phototransduction
cone_loopg		0.0	Gain of calcium feedback loop in cones
cone_maxcond	$n S$	0.2	Max. conductance of of outer segment
timinc	$m s$	0.001 \| 0.01	Simulation time step (Figure 9)
ploti	$m s$	0.01 \| 1	Recording time step (Figure 9 and Figure 10 \| otherwise)
stiminc	$m s$	0.01 \| 0.1	Synaptic time step (Figure 9 and Figure 10 \| otherwise)
srtimestep	$m s$	0.001 \| 0.01 \| 0.1	Stimulus update time step (Figure 9 \| Figure 10 \| otherwise)
predur	$s$	5 \| ≥10	Simulation time to reach equilibrium potential (Cone \| BCs)

Data availability

Models and simulation code is available at https://github.com/berenslab/CBC_inference (copy archived at https://archive.softwareheritage.org/swh:1:rev:2b8ec4ac0ca916d42cba0404229298f8ff79c3a3/). Experimental and inference data is available at https://zenodo.org/record/4185955.

The following data sets were generated

1. Oesterle J
2. Behrens C
3. Schröder C
4. Herrmann T
5. Euler T
6. Franke K
7. Smith RG
8. Zeck G
9. Berens P
(2020) Zenodo
Data for "Bayesian inference for biophysical neuron models enables stimulus optimization for retinal neuroprosthetics".

https://doi.org/10.5281/zenodo.4185955

References

1. Achard P
2. De Schutter E
(2006) Complex parameter landscape for a complex neuron model
PLOS Computational Biology 2:e94.

https://doi.org/10.1371/journal.pcbi.0020094
- PubMed
- Google Scholar
1. Alonso LM
2. Marder E
(2019) Visualization of currents in neural models with similar behavior and different conductance densities
eLife 8:e42722.

https://doi.org/10.7554/eLife.42722
- PubMed
- Google Scholar
(2001) Integrated allosteric model of voltage gating of HCN channels
Journal of General Physiology 117:519–532.

https://doi.org/10.1085/jgp.117.6.519
- PubMed
- Google Scholar
(2013) Spikes and ribbon synapses in early vision
Trends in Neurosciences 36:480–488.

https://doi.org/10.1016/j.tins.2013.04.006
- PubMed
- Google Scholar
1. Baden T
2. Nikolaev A
3. Esposti F
4. Dreosti E
5. Odermatt B
6. Lagnado L
(2014) A synaptic mechanism for temporal filtering of visual signals
PLOS Biology 12:e1001972.

https://doi.org/10.1371/journal.pbio.1001972
- PubMed
- Google Scholar
(2017) Visual prosthesis: interfacing stimulating electrodes with retinal neurons to restore vision
Frontiers in Neuroscience 11:620.

https://doi.org/10.3389/fnins.2017.00620
- PubMed
- Google Scholar
(2010) Vesicle pool size at the salamander cone ribbon synapse
Journal of Neurophysiology 103:419–423.

https://doi.org/10.1152/jn.00718.2009
- PubMed
- Google Scholar
1. Behrens C
2. Schubert T
3. Haverkamp S
4. Euler T
5. Berens P
(2016) Connectivity map of bipolar cells and photoreceptors in the mouse retina
eLife 5:e20041.

https://doi.org/10.7554/eLife.20041
- PubMed
- Google Scholar
1. Berntson A
2. Taylor WR
(2003) The unitary event amplitude of mouse retinal on-cone bipolar cells
Visual Neuroscience 20:621–626.

https://doi.org/10.1017/S0952523803206040
- PubMed
- Google Scholar
Preprint
1. Bittner SR
2. Palmigiano A
3. Piet AT
4. Duan CA
5. Brody CD
6. Miller KD
7. Cunningham JP
(2019) Interrogating theoretical models of neural computation with deep inference
bioRxiv.

https://doi.org/10.1101/837567
- Google Scholar
(2012) The photovoltage of rods and cones in the dark-adapted mouse retina
The Journal of Physiology 590:3841–3855.

https://doi.org/10.1113/jphysiol.2011.226878
- PubMed
- Google Scholar
1. Caputo A
2. Piano I
3. Demontis GC
4. Bacchi N
5. Casarosa S
6. Della Santina L
7. Gargini C
(2015) TMEM16A is associated with voltage-gated calcium channels in mouse retina and its function is disrupted upon mutation of the auxiliary α2δ4 subunit
Frontiers in Cellular Neuroscience 9:422.

https://doi.org/10.3389/fncel.2015.00422
- PubMed
- Google Scholar
1. Carter-Dawson LD
2. LaVail MM
(1979) Rods and cones in the mouse retina. I. structural analysis using light and electron microscopy
The Journal of Comparative Neurology 188:245–262.

https://doi.org/10.1002/cne.901880204
- PubMed
- Google Scholar
1. Choi SY
2. Borghuis BG
3. Borghuis B
4. Rea R
5. Levitan ES
6. Sterling P
7. Kramer RH
(2005) Encoding light intensity by the cone photoreceptor synapse
Neuron 48:555–562.

https://doi.org/10.1016/j.neuron.2005.09.011
- PubMed
- Google Scholar
(2014) Retinal implants: a systematic review
British Journal of Ophthalmology 98:852–856.

https://doi.org/10.1136/bjophthalmol-2013-303708
- PubMed
- Google Scholar
1. Clancy CE
2. Kass RS
(2004) Theoretical investigation of the neuronal na+ channel SCN1A: abnormal gating and epilepsy
Biophysical Journal 86:2606–2614.

https://doi.org/10.1016/S0006-3495(04)74315-5
- PubMed
- Google Scholar
Software
1. Comsol AB
(2019) Comsol Multiphysics, version 5.4
Comsol Multiphysics, Stockholm, Sweden.

http://www.comsol.com
(2018) Electrode-size dependent thresholds in subretinal neuroprosthetic stimulation
Journal of Neural Engineering 15:045003.

https://doi.org/10.1088/1741-2552/aac1c8
- PubMed
- Google Scholar
1. Cui J
2. Ivanova E
3. Qi L
4. Pan ZH
(2012) Expression of CaV3.2 T-type Ca²⁺ channels in a subpopulation of retinal type-3 cone bipolar cells
Neuroscience 224:63–69.

https://doi.org/10.1016/j.neuroscience.2012.08.017
- PubMed
- Google Scholar
1. Cui J
2. Pan ZH
(2008) Two types of cone bipolar cells express voltage-gated Na⁺ channels in the rat retina
Visual Neuroscience 25:635–645.

https://doi.org/10.1017/S0952523808080851
- PubMed
- Google Scholar
1. Dong CJ
2. Werblin FS
(1995) Inwardly rectifying potassium conductance can accelerate the hyperpolarizing response in retinal horizontal cells
Journal of Neurophysiology 74:2258–2265.

https://doi.org/10.1152/jn.1995.74.6.2258
- PubMed
- Google Scholar
Preprint
(2018) Sequential neural methods for likelihood-free inference
arXiv.

https://arxiv.org/abs/1811.08723
- Google Scholar
(2018) Assessment of the electronic retinal implant alpha AMS in restoring vision to blind patients with End-Stage retinitis pigmentosa
Ophthalmology 125:432–443.

https://doi.org/10.1016/j.ophtha.2017.09.019
- PubMed
- Google Scholar
(2012) Electrical stimulation of retinal neurons in epiretinal and subretinal configuration using a multicapacitor array
Journal of Neurophysiology 107:2742–2755.

https://doi.org/10.1152/jn.00909.2011
- PubMed
- Google Scholar
1. Eickenscheidt M
2. Zeck G
(2014) Action potentials in retinal ganglion cells are initiated at the site of maximal curvature of the extracellular potential
Journal of Neural Engineering 11:036006.

https://doi.org/10.1088/1741-2560/11/3/036006
- PubMed
- Google Scholar
Book
1. Euler T
2. Franke K
3. Baden T
(2019) Studying a light sensor with light: multiphoton imaging in the retina
In: Hartveit E, editors. Multiphoton Microscopy. Springer. pp. 225–250.

https://doi.org/10.1007/978-1-4939-9702-2_10
- Google Scholar
1. Fohlmeister JF
2. Miller RF
(1997) Impulse encoding mechanisms of ganglion cells in the tiger salamander retina
Journal of Neurophysiology 78:1935–1947.

https://doi.org/10.1152/jn.1997.78.4.1935
- PubMed
- Google Scholar
1. Franke K
2. Berens P
3. Schubert T
4. Bethge M
5. Euler T
6. Baden T
(2017) Inhibition decorrelates visual feature representations in the inner retina
Nature 542:439–444.

https://doi.org/10.1038/nature21394
- PubMed
- Google Scholar
(1996) The dielectric properties of biological tissues: iii. Parametric models for the dielectric spectrum of tissues
Physics in Medicine and Biology 41:2271–2293.

https://doi.org/10.1088/0031-9155/41/11/003
- PubMed
- Google Scholar
Conference
(2010) Electric field stimulation of bipolar cells in a degenerated retina--a theoretical study
IEEE Transactions on Neural Systems and Rehabilitation Engineering. pp. 1–10.

https://doi.org/10.1109/TNSRE.2009.2037323
- Google Scholar
(2006) Genetic algorithm for optimization and specification of a neuron model
Neurocomputing 69:1039–1042.

https://doi.org/10.1016/j.neucom.2005.12.041
- Google Scholar
Book
1. Gerstner W
2. Kistler WM
(2002) Spiking Neuron Models: Single Neurons, Populations, Plasticity
Cambridge university press.

https://doi.org/10.1017/CBO9780511815706
- Google Scholar
(2001) Dichotomy of action-potential backpropagation in CA1 pyramidal neuron dendrites
Journal of Neurophysiology 86:2998–3010.

https://doi.org/10.1152/jn.2001.86.6.2998
- PubMed
- Google Scholar
(2001) Global structure, robustness, and modulation of neuronal models
The Journal of Neuroscience 21:5229–5238.

https://doi.org/10.1523/JNEUROSCI.21-14-05229.2001
- PubMed
- Google Scholar
(2020) Training deep neural density estimators to identify mechanistic models of neural dynamics
eLife 9:e56261.

https://doi.org/10.7554/eLife.56261
- PubMed
- Google Scholar
1. Greene MJ
2. Kim JS
3. Seung HS
4. EyeWirers
(2016) Analogous convergence of sustained and transient inputs in parallel on and off pathways for retinal motion computation
Cell Reports 14:1892–1900.

https://doi.org/10.1016/j.celrep.2016.02.001
- PubMed
- Google Scholar
1. Gutmann MU
2. Corander J
(2016)
Bayesian optimization for likelihood-free inference of simulator-based statistical models

The Journal of Machine Learning Research 17:4256–4302.
- Google Scholar
1. Hay E
2. Hill S
3. Schürmann F
4. Markram H
5. Segev I
(2011) Models of neocortical layer 5b pyramidal cells capturing a wide range of dendritic and perisomatic active properties
PLOS Computational Biology 7:e1002107.

https://doi.org/10.1371/journal.pcbi.1002107
- PubMed
- Google Scholar
1. Hellmer CB
2. Zhou Y
3. Fyk-Kolodziej B
4. Hu Z
5. Ichinose T
(2016) Morphological and physiological analysis of type-5 and other bipolar cells in the mouse retina
Neuroscience 315:246–258.

https://doi.org/10.1016/j.neuroscience.2015.12.016
- PubMed
- Google Scholar
1. Helmstaedter M
2. Briggman KL
3. Turaga SC
4. Jain V
5. Seung HS
6. Denk W
(2013) Connectomic reconstruction of the inner plexiform layer in the mouse retina
Nature 500:168–174.

https://doi.org/10.1038/nature12346
- PubMed
- Google Scholar
(1998) Gating of recombinant small-conductance Ca-activated K+ channels by calcium
Journal of General Physiology 111:565–581.

https://doi.org/10.1085/jgp.111.4.565
- PubMed
- Google Scholar
1. Hodgkin AL
2. Huxley AF
(1952) A quantitative description of membrane current and its application to conduction and excitation in nerve
The Journal of Physiology 117:500–544.

https://doi.org/10.1113/jphysiol.1952.sp004764
- PubMed
- Google Scholar
1. Hu C
2. Bi A
3. Pan ZH
(2009) Differential expression of three T-type calcium channels in retinal bipolar cells in rats
Visual Neuroscience 26:177–187.

https://doi.org/10.1017/S0952523809090026
- PubMed
- Google Scholar
(2014) Roles of ON cone bipolar cell subtypes in temporal coding in the mouse retina
Journal of Neuroscience 34:8761–8771.

https://doi.org/10.1523/JNEUROSCI.3965-13.2014
- PubMed
- Google Scholar
1. Ichinose T
2. Hellmer CB
(2016) Differential signalling and glutamate receptor compositions in the OFF bipolar cell types in the mouse retina
The Journal of Physiology 594:883–894.

https://doi.org/10.1113/JP271458
- PubMed
- Google Scholar
(2020) Membrane conductances of mouse cone photoreceptors
Journal of General Physiology 152:e201912520.

https://doi.org/10.1085/jgp.201912520
- PubMed
- Google Scholar
1. Ivanova E
2. Müller F
(2006) Retinal bipolar cell types differ in their inventory of ion channels
Visual Neuroscience 23:143–154.

https://doi.org/10.1017/S0952523806232048
- PubMed
- Google Scholar
(2005) Thresholds for activation of rabbit retinal ganglion cells with relatively large, extracellular microelectrodes
Investigative Opthalmology & Visual Science 46:1486–1496.

https://doi.org/10.1167/iovs.04-1018
- Google Scholar
1. Jensen RJ
2. Rizzo JF
(2008) Activation of retinal ganglion cells in wild-type and rd1 mice through electrical stimulation of the retinal neural network
Vision Research 48:1562–1568.

https://doi.org/10.1016/j.visres.2008.04.016
- PubMed
- Google Scholar
(1993) Quantal components of unitary EPSCs at the mossy fibre synapse on CA3 pyramidal cells of rat Hippocampus
The Journal of Physiology 472:615–663.

https://doi.org/10.1113/jphysiol.1993.sp019965
- PubMed
- Google Scholar
1. Karschin A
2. Lipton SA
(1989) Calcium channels in solitary retinal ganglion cells from post-natal rat
The Journal of Physiology 418:379–396.

https://doi.org/10.1113/jphysiol.1989.sp017847
- PubMed
- Google Scholar
1. Karwoski CJ
2. Xu X
(1999) Current source-density analysis of light-evoked field potentials in rabbit retina
Visual Neuroscience 16:369–377.

https://doi.org/10.1017/S0952523899162163
- PubMed
- Google Scholar
1. Kavalali ET
(2015) The mechanisms and functions of spontaneous neurotransmitter release
Nature Reviews Neuroscience 16:5–16.

https://doi.org/10.1038/nrn3875
- PubMed
- Google Scholar
(2005) Constraining compartmental models using multiple voltage recordings and genetic algorithms
Journal of Neurophysiology 94:3730–3742.

https://doi.org/10.1152/jn.00408.2005
- PubMed
- Google Scholar
1. Knop GC
2. Seeliger MW
3. Thiel F
4. Mataruga A
5. Kaupp UB
6. Friedburg C
7. Tanimoto N
8. Müller F
(2008) Light responses in the mouse retina are prolonged upon targeted deletion of the HCN1 channel gene
European Journal of Neuroscience 28:2221–2230.

https://doi.org/10.1111/j.1460-9568.2008.06512.x
- PubMed
- Google Scholar
Book
1. Koch C
(2004)
Biophysics of Computation: Information Processing in Single Neurons

Oxford University Press.
- Google Scholar
1. Lagnado L
2. Schmitz F
(2015) Ribbon synapses and visual processing in the retina
Annual Review of Vision Science 1:235–262.

https://doi.org/10.1146/annurev-vision-082114-035709
- PubMed
- Google Scholar
(2003) Availability of low-threshold Ca2+ current in retinal ganglion cells
Journal of Neurophysiology 90:3888–3901.

https://doi.org/10.1152/jn.00477.2003
- PubMed
- Google Scholar
1. Lee JI
2. Im M
(2019) Optimal electric stimulus amplitude improves the selectivity between responses of ON versus OFF types of retinal ganglion cells
IEEE Transactions on Neural Systems and Rehabilitation Engineering 27:2015–2024.

https://doi.org/10.1109/TNSRE.2019.2939012
- PubMed
- Google Scholar
(2017) Fundamentals and recent developments in approximate bayesian computation
Systematic Biology 66:e66–e82.

https://doi.org/10.1093/sysbio/syw077
- PubMed
- Google Scholar
Conference
(2017)
Flexible statistical inference for mechanistic models of neural dynamics

Advances in Neural Information Processing Systems. pp. 1289–1299.
- Google Scholar
1. Luo YH
2. da Cruz L
(2016) The argus II retinal prosthesis system
Progress in Retinal and Eye Research 50:89–107.

https://doi.org/10.1016/j.preteyeres.2015.09.003
- PubMed
- Google Scholar
1. Ma YP
2. Cui J
3. Pan ZH
(2005) Heterogeneous expression of voltage-dependent Na⁺ and K⁺ channels in mammalian retinal bipolar cells
Visual Neuroscience 22:119–133.

https://doi.org/10.1017/S0952523805222010
- PubMed
- Google Scholar
Software
1. Macke Lab
(2020) DELFI, version 0.5.1
GitHub.

https://github.com/mackelab/delfi
1. Mansergh F
2. Orton NC
3. Vessey JP
4. Lalonde MR
5. Stell WK
6. Tremblay F
7. Barnes S
8. Rancourt DE
9. Bech-Hansen NT
(2005) Mutation of the calcium channel gene Cacna1f disrupts calcium signaling, synaptic transmission and cellular organization in mouse retina
Human Molecular Genetics 14:3035–3046.

https://doi.org/10.1093/hmg/ddi336
- PubMed
- Google Scholar
1. Marvin JS
2. Borghuis BG
3. Tian L
4. Cichon J
5. Harnett MT
6. Akerboom J
7. Gordus A
8. Renninger SL
9. Chen TW
10. Bargmann CI
11. Orger MB
12. Schreiter ER
13. Demb JB
14. Gan WB
15. Hires SA
16. Looger LL
(2013) An optimized fluorescent probe for visualizing glutamate neurotransmission
Nature Methods 10:162–170.

https://doi.org/10.1038/nmeth.2333
- PubMed
- Google Scholar
1. Masland RH
(2012) The tasks of amacrine cells
Visual Neuroscience 29:3–9.

https://doi.org/10.1017/S0952523811000344
- PubMed
- Google Scholar
(2019) Spatiotemporally asymmetric excitation supports mammalian retinal motion sensitivity
Current Biology 29:3277–3288.

https://doi.org/10.1016/j.cub.2019.08.048
- PubMed
- Google Scholar
1. Matthews G
2. Fuchs P
(2010) The diverse roles of ribbon synapses in sensory neurotransmission
Nature Reviews Neuroscience 11:812–822.

https://doi.org/10.1038/nrn2924
- PubMed
- Google Scholar
(1998) Calcium extrusion from mammalian photoreceptor terminals
The Journal of Neuroscience 18:2467–2474.

https://doi.org/10.1523/JNEUROSCI.18-07-02467.1998
- PubMed
- Google Scholar
1. Morgans CW
2. Bayley PR
3. Oesch NW
4. Ren G
5. Akileswaran L
6. Taylor WR
(2005) Photoreceptor calcium channels: insight from night blindness
Visual Neuroscience 22:561–568.

https://doi.org/10.1017/S0952523805225038
- PubMed
- Google Scholar
1. Müller F
2. Scholten A
3. Ivanova E
4. Haverkamp S
5. Kremmer E
6. Kaupp UB
(2003) HCN channels are expressed differentially in retinal bipolar cells and concentrated at synaptic terminals
European Journal of Neuroscience 17:2084–2096.

https://doi.org/10.1046/j.1460-9568.2003.02634.x
- PubMed
- Google Scholar
(2012) Neuromantic - from semi-manual to semi-automatic reconstruction of neuron morphology
Frontiers in Neuroinformatics 6:4.

https://doi.org/10.3389/fninf.2012.00004
- PubMed
- Google Scholar
(1998) Kinetics of recovery of the dark-adapted salamander rod photoresponse
Journal of General Physiology 111:7–37.

https://doi.org/10.1085/jgp.111.1.7
- PubMed
- Google Scholar
1. Nikonov SS
2. Kholodenko R
3. Lem J
4. Pugh EN
(2006) Physiological features of the S- and M-cone photoreceptors of wild-type mice from single-cell recordings
Journal of General Physiology 127:359–374.

https://doi.org/10.1085/jgp.200609490
- PubMed
- Google Scholar
(2014) Cell types, network homeostasis, and pathological compensation from a biologically plausible ion channel expression model
Neuron 82:809–821.

https://doi.org/10.1016/j.neuron.2014.04.002
- PubMed
- Google Scholar
Software
1. Oesterle J
(2020) CBC_inference, version swh:1:rev:2b8ec4ac0ca916d42cba0404229298f8ff79c3a3
Software Heritage.

https://archive.softwareheritage.org/browse/directory/74cb51f26fd3961f853488aeb3b8ef959ce4e541/?branch=refs/heads/master&origin_url=https://github.com/berenslab/CBC_inference&snapshot=e584de7e33f4121abd7232331b53d84b5863674b
(2009) Passive membrane properties and electrotonic signal processing in retinal rod bipolar cells
The Journal of Physiology 587:829–849.

https://doi.org/10.1113/jphysiol.2008.165415
- PubMed
- Google Scholar
Preprint
(2018) Sequential neural likelihood: fast likelihood-free inference with autoregressive flows
arXiv.

https://arxiv.org/abs/1805.07226
- Google Scholar
Conference
1. Papamakarios G
2. Murray I
(2016)
Fast ε-free inference of simulation models with bayesian conditional density estimation

Advances in Neural Information Processing Systems. pp. 1028–1036.
- Google Scholar
1. Pennesi ME
2. Michaels KV
3. Magee SS
4. Maricle A
5. Davin SP
6. Garg AK
7. Gale MJ
8. Tu DC
9. Wen Y
10. Erker LR
11. Francis PJ
(2012) Long-term characterization of retinal degeneration in rd1 and rd10 mice using spectral domain optical coherence tomography
Investigative Ophthalmology & Visual Science 53:4644–4656.

https://doi.org/10.1167/iovs.12-9611
- PubMed
- Google Scholar
(2003) Pyramidal neuron as two-layer neural network
Neuron 37:989–999.

https://doi.org/10.1016/S0896-6273(03)00149-1
- PubMed
- Google Scholar
(2003) Alternative to hand-tuning conductance-based models: construction and analysis of databases of model neurons
Journal of Neurophysiology 90:3998–4015.

https://doi.org/10.1152/jn.00641.2003
- PubMed
- Google Scholar
(2017) Impact of electrode position on the elicitation of sodium spikes in retinal bipolar cells
Scientific Reports 7:17590.

https://doi.org/10.1038/s41598-017-17603-8
- PubMed
- Google Scholar
1. Resatz S
2. Rattay F
(2004) A model for the electrically stimulated retina
Mathematical and Computer Modelling of Dynamical Systems 10:93–106.

https://doi.org/10.1080/13873950412331318080
- Google Scholar
(2011) Fitting neuron models to spike trains
Frontiers in Neuroscience 5:9.

https://doi.org/10.3389/fnins.2011.00009
- PubMed
- Google Scholar
1. Saszik S
2. DeVries SH
(2012) A mammalian retinal bipolar cell uses both graded changes in membrane voltage and all-or-nothing na+ spikes to encode light
Journal of Neuroscience 32:297–307.

https://doi.org/10.1523/JNEUROSCI.2739-08.2012
- PubMed
- Google Scholar
1. Satoh H
2. Aoki K
3. Watanabe SI
4. Kaneko A
(1998) L-type calcium channels in the axon terminal of mouse bipolar cells
NeuroReport 9:2161–2165.

https://doi.org/10.1097/00001756-199807130-00002
- PubMed
- Google Scholar
Conference
(2019)
Approximate bayesian inference for a mechanistic model of vesicle release at a ribbon synapse

Advances in Neural Information Processing Systems.
- Google Scholar
1. Sheng Z
2. Choi SY
3. Dharia A
4. Li J
5. Sterling P
6. Kramer RH
(2007) Synaptic Ca2+ in darkness is lower in rods than cones, causing slower tonic release of vesicles
Journal of Neuroscience 27:5033–5042.

https://doi.org/10.1523/JNEUROSCI.5386-06.2007
- PubMed
- Google Scholar
1. Singer JH
2. Diamond JS
(2006) Vesicle depletion and synaptic depression at a mammalian ribbon synapse
Journal of Neurophysiology 95:3191–3198.

https://doi.org/10.1152/jn.01309.2005
- PubMed
- Google Scholar
Book
(2018)
Handbook of Approximate Bayesian Computation

Chapman and Hall/CRC.
- Google Scholar
1. Smith RG
(1992) NeuronC: a computational language for investigating functional architecture of neural circuits
Journal of Neuroscience Methods 43:83–108.

https://doi.org/10.1016/0165-0270(92)90019-A
- PubMed
- Google Scholar
1. Stringer C
2. Pachitariu M
3. Steinmetz NA
4. Okun M
5. Bartho P
6. Harris KD
7. Sahani M
8. Lesica NA
(2016) Inhibitory control of correlated intrinsic variability in cortical networks
eLife 5:e19695.

https://doi.org/10.7554/eLife.19695
- PubMed
- Google Scholar
Preprint
1. Szatko KP
2. Korympidou MM
3. Ran Y
4. Berens P
5. Dalkara D
6. Schubert T
7. Euler T
8. Franke K
(2019) Neural circuits in the mouse retina support color vision in the upper visual field
bioRxiv.

https://doi.org/10.1101/745539
- Google Scholar
(2009) How multiple conductances determine electrophysiological properties in a multicompartment model
Journal of Neuroscience 29:5573–5586.

https://doi.org/10.1523/JNEUROSCI.4438-08.2009
- PubMed
- Google Scholar
(2016) Modeling and measurement of vesicle pools at the cone ribbon synapse: changes in release probability are solely responsible for voltage-dependent changes in release
Synapse 70:1–14.

https://doi.org/10.1002/syn.21871
- PubMed
- Google Scholar
1. tom Dieck S
2. Brandstätter JH
(2006) Ribbon synapses of the retina
Cell and Tissue Research 326:339–346.

https://doi.org/10.1007/s00441-006-0234-0
- PubMed
- Google Scholar
1. Tsai D
2. Morley JW
3. Suaning GJ
4. Lovell NH
(2009) Direct activation and temporal response properties of rabbit retinal ganglion cells following subretinal stimulation
Journal of Neurophysiology 102:2982–2993.

https://doi.org/10.1152/jn.00545.2009
- PubMed
- Google Scholar
1. Twyford P
2. Cai C
3. Fried S
(2014) Differential responses to high-frequency electrical stimulation in ON and OFF retinal ganglion cells
Journal of Neural Engineering 11:025001.

https://doi.org/10.1088/1741-2560/11/2/025001
- PubMed
- Google Scholar
(2019) Voltage- and calcium-gated ion channels of neurons in the vertebrate retina
Progress in Retinal and Eye Research 72:100760.

https://doi.org/10.1016/j.preteyeres.2019.05.001
- PubMed
- Google Scholar
1. Wan QF
2. Heidelberger R
(2011) Synaptic release at mammalian bipolar cell terminals
Visual Neuroscience 28:109–119.

https://doi.org/10.1017/S0952523810000453
- PubMed
- Google Scholar
Conference
1. Wang B
2. Weiland JD
(2015) Resistivity profiles of wild-type, rd1, and rd10 mouse retina
37th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC) IEEE. pp. 1650–1653.

https://doi.org/10.1109/EMBC.2015.7318692
- Google Scholar
(2009) Cone contacts, mosaics, and territories of bipolar cells in the mouse retina
Journal of Neuroscience 29:106–117.

https://doi.org/10.1523/JNEUROSCI.4442-08.2009
- PubMed
- Google Scholar
(2015) Improving the spatial resolution of epiretinal implants by increasing stimulus pulse duration
Science Translational Medicine 7:318ra203.

https://doi.org/10.1126/scitranslmed.aac4877
- PubMed
- Google Scholar
1. Werginz P
2. Benav H
3. Zrenner E
4. Rattay F
(2015) Modeling the response of ON and OFF retinal bipolar cells during electric stimulation
Vision Research 111:170–181.

https://doi.org/10.1016/j.visres.2014.12.002
- PubMed
- Google Scholar
1. Wood SN
(2010) Statistical inference for noisy nonlinear ecological dynamic systems
Nature 466:1102–1104.

https://doi.org/10.1038/nature09319
- PubMed
- Google Scholar
1. Xiao H
2. Peng H
(2013) APP2: automatic tracing of 3D neuron morphology based on hierarchical pruning of a gray-weighted image distance-tree
Bioinformatics 29:1448–1454.

https://doi.org/10.1093/bioinformatics/btt170
- PubMed
- Google Scholar
1. Yang YD
2. Cho H
3. Koo JY
4. Tak MH
5. Cho Y
6. Shim WS
7. Park SP
8. Lee J
9. Lee B
10. Kim BM
11. Raouf R
12. Shin YK
13. Oh U
(2008) TMEM16A confers receptor-activated calcium-dependent chloride conductance
Nature 455:1210–1215.

https://doi.org/10.1038/nature07313
- PubMed
- Google Scholar
Preprint
1. Zhao Z
2. Klindt D
3. Chagas AM
4. Szatko KP
5. Rogerson L
6. Protti D
7. Behrens C
8. Dalkara D
9. Schubert T
10. Bethge M
(2019) The temporal structure of the inner retina at a single glance
bioRxiv.

https://doi.org/10.1101/743047
- Google Scholar
1. Zrenner E
(2002) Will retinal implants restore vision?
Science 295:1022–1025.

https://doi.org/10.1126/science.1067996
- PubMed
- Google Scholar

Article and author information

Author details

Jonathan Oesterle

Institute for Ophthalmic Research, University of Tübingen, Tübingen, Germany

Contribution
Software, Formal analysis, Investigation, Visualization, Methodology, Writing - original draft, Writing - review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-8919-1445
Christian Behrens

Institute for Ophthalmic Research, University of Tübingen, Tübingen, Germany

Contribution
Supervision, Investigation, Writing - review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-3623-352X
Cornelius Schröder

Institute for Ophthalmic Research, University of Tübingen, Tübingen, Germany

Contribution
Methodology, Writing - review and editing

Competing interests
No competing interests declared
Thoralf Hermann

Naturwissenschaftliches und Medizinisches Institut an der Universität Tübingen, Reutlingen, Germany

Contribution
Data curation, Investigation

Competing interests
No competing interests declared
Thomas Euler
1. Institute for Ophthalmic Research, University of Tübingen, Tübingen, Germany
2. Center for Integrative Neuroscience, University of Tübingen, Tübingen, Germany
3. Bernstein Center for Computational Neuroscience, University of Tübingen, Tübingen, Germany
Contribution
Supervision, Writing - review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-4567-6966
Katrin Franke
1. Institute for Ophthalmic Research, University of Tübingen, Tübingen, Germany
2. Bernstein Center for Computational Neuroscience, University of Tübingen, Tübingen, Germany
Contribution
Data curation, Investigation

Competing interests
No competing interests declared
Robert G Smith

Department of Neuroscience, University of Pennsylvania, Philadelphia, United States

Contribution
Conceptualization, Resources, Funding acquisition, Investigation, Writing - review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-5703-1324
Günther Zeck

Naturwissenschaftliches und Medizinisches Institut an der Universität Tübingen, Reutlingen, Germany

Contribution
Conceptualization, Supervision, Funding acquisition, Writing - review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-3998-9883
Philipp Berens
1. Institute for Ophthalmic Research, University of Tübingen, Tübingen, Germany
2. Center for Integrative Neuroscience, University of Tübingen, Tübingen, Germany
3. Bernstein Center for Computational Neuroscience, University of Tübingen, Tübingen, Germany
4. Institute for Bioinformatics and Medical Informatics, University of Tübingen, Tübingen, Germany
Contribution
Conceptualization, Supervision, Funding acquisition, Writing - original draft, Writing - review and editing

For correspondence
philipp.berens@uni-tuebingen.de

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-0199-4727

Funding

Bundesministerium für Bildung und Forschung (01GQ1601)

Philipp Berens

Bundesministerium für Bildung und Forschung (01IS18052)

Philipp Berens

Deutsche Forschungsgemeinschaft (EXC 2064 - 390727645)

Philipp Berens

Deutsche Forschungsgemeinschaft (BE5601/4-1)

Philipp Berens

Baden-Württemberg Stiftung (NEU013)

Günther Zeck
Philipp Berens

National Institutes of Health (EY022070)

Robert G Smith

National Institutes of Health (EY023766)

Robert G Smith

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank Yves Bernaerts for preliminary work on electrical stimulation of biophysical models, Timm Schubert for support with animal protocols and Gordon Eske and Adam vor Daranyi for technical support. Additionally, we thank Jakob Macke, Philipp Hennig and Lara Hoefling for discussion. This work was funded by the Federal Ministry of Education and Research (BMBF, FKZ 01GQ1601 and 01IS18052 to PB), the German Science Foundation through the Excellence Cluster 2064 ‘Machine Learning - New Perspectives for Science’ (390727645), a Heisenberg Professorship (BE5601/4-1 to PB), the Baden-Württemberg Stiftung (NEU013 to PB and GZ), and the National Institutes of Health (NIH) (EY022070 and EY023766 to RGS).

Ethics

Animal experimentation: All procedures were approved by the governmental review board (Regierungspräsidium Tübingen, Baden-Württemberg, Konrad-Adenauer-Str. 20, 72072 Tübingen, Germany, AZ 35/9185.82-7) and performed according to the laws governing animal experimentation issued by the German Government.

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.