Stochastic inference with spiking neurons in the high-conductance state∗
Mihai A. Petrovici1 *, Johannes Bill1,2 *, Ilja Bytschok1 , Johannes Schemmel1 , Karlheinz Meier1

arXiv:1610.07161v1 [q-bio.NC] 23 Oct 2016

2

1
Kirchhoff Institute for Physics, University of Heidelberg
Institute for Theoretical Computer Science, Graz University of Technology
(Dated: October 25, 2016)

The highly variable dynamics of neocortical circuits observed in vivo have been hypothesized
to represent a signature of ongoing stochastic inference, but stand in apparent contrast to the
deterministic response of neurons measured in vitro. Based on a propagation of the membrane
autocorrelation across spike bursts, we provide an analytical derivation of the neural activation
function that holds for a large parameter space, including the high-conductance state. On this
basis, we show how an ensemble of leaky integrate-and-fire neurons with conductance-based synapses
embedded in a spiking environment can attain the correct firing statistics for sampling from a welldefined target distribution. For recurrent networks, we examine convergence towards stationarity in
computer simulations and demonstrate sample-based Bayesian inference in a mixed graphical model.
This points to a new computational role of high-conductance states and establishes a rigorous link
between deterministic neuron models and functional stochastic dynamics on the network level.

Introduction

In responding to environmental stimuli, brains have
to make predictions based on incomplete, noisy and ambiguous data. The recent hypothesis that the brain copes
with this challenge by performing Bayesian, rather than
logical inference [1–3], has been strengthened by electrophysiological data which identified neural correlates of
the involved computations [4, 5] and theoretical work on
potential spiking network implementations [6–8].
In probabilistic inference, possible values of a quantity are described by a random variable (RV) zk and
all dependencies between RVs are stored in a joint
distribution p(z1 , . . . , zK ). The belief about a set of
unobserved RVs {z1 , . . . , zM } given an observed set
of RVs is represented by the posterior distribution
p(z1 , . . . , zM | zM +1 , . . . , zK ). In particular, the posterior
contains information on the most likely conclusion, as
well as on all potential alternatives.
With regard to the representation of probability distributions in the brain, theoretical work [2] has argued
in favor of sample-based codes, i.e., instead of providing the entire distribution at any point in time, samples
z (t) ∼ p(z1 , . . . , zK ) are used as a proxy. When modeling large systems, this offers three important advantages.
First, approximate solutions can be provided at any time,
with increasingly reliable results as the calculation progresses (“anytime computing”). Secondly, marginalization comes at no cost, as p(zk ) can be determined by simply neglecting the values of all other RVs. Thirdly, some
sampling algorithms support a high degree of parallelization with an algorithmic structure that is reminiscent
of neural networks [9]. Recently, a theory has been suggested which combines these advantages by implementing
Markov chain Monte Carlo (MCMC) sampling in spiking
networks of abstract model neurons [7]. In this framework, spike patterns are interpreted as samples of binary

FIG. 1. (A) Spike patterns as samples of a random vector z.
The RV zk is active for duration τon (gray bar) after each
spike. (B) In stochastic neuron models, internal state variables (red) modulate the firing probability (gray). (C) In
contrast, deterministic neurons elicit a spike when the membrane potential (blue) reaches a threshold voltage.

RVs as follows (see Fig. 1A):
(t)

zk = 1 ⇔ Neuron k fired in (t − τon , t] .

(1)

The duration τon of the active state following a spike is a
free parameter for the timescale on which a neuron affects
downstream cells. The neuron model in [7] is inherently
stochastic, with an instantaneous firing probability
 exp(v )
k
if zk = 0
p(spike in [t, t+∆t))
τ
rk (t) = lim
=
(2)
∆t
∆t→0
0
if zk = 1
where vk is an abstract membrane potential (Fig. 1B).
In contrast to this idealized stochastic neuron model,
in vitro experiments have demonstrated the largely deterministic nature of single neurons [10]. Similarly, microscopic models of neural circuits typically rely on deterministic dynamics of their constituents (see Fig. 1C).
The aim of this letter is to demonstrate how a network of
deterministic neurons in a biologically plausible spiking
environment can quantitatively reproduce the stochastic
dynamics required for sampling from a well-defined distribution and perform inference given observations. In
this letter, we extend our earlier discussion of this approach [11].
We start by calculating the dynamics of a single
leaky integrate-and-fire (LIF) neuron in a spiking noisy

2
environment and derive its activation function in the
high-conductance state. For this biologically relevant
regime, existing analytical descriptions of neuronal response functions (e.g., [16, 17]) are not applicable. Here,
we provide an approach based on the propagation of
the membrane autocorrelation throughout spike bursts.
This establishes an equivalence to the abstract, inherently stochastic units (Eq. 2). On the network level, we
show how conductance-based synapses approximate the
interaction for sampling from a well-defined target distribution. Furthermore, we show that the distribution
sampled by LIF networks remains a good approximation
of the target distribution even for large networks with
strong recurrent interaction. We complement our study
with a demonstration of probabilistic inference by implementing a small graphical model for pattern recognition
in a recurrent network of LIF neurons.

LIF dynamics in a high-conductance state

We consider deterministic LIF neurons defined
d
uk = gl (El − uk ) + Ik , with membrane potenby Cm dt
tial uk , capacitance Cm , leak potential El , leak conductance gl and input current Ik . When uk crosses a
threshold ϑ from below, a spike is emitted and uk is reset to % for a refractory period τref . We formally partition the total input current Ik into recurrent synaptic
input, diffuse synaptic noise and additional external currents: Ik = Ikrec + Iknoise + Ikext . The currents Ikrec and
noise
IP
are mediated through synapses and obey Iksyn =
k
rev
− uk ) with reversal potential Eirev of the ith
i gki (Ei
synapse. The temporal evolution of the conductance gki
is modelled as a low-pass
filter on pre-synaptic spikes:
P
d
g
=
−g
/τ
+
w
δ(t − tsi ), with synaptic time
ki
ki
syn
ki
s
dt
constant τsyn , weight wki and spike times tsi .
We start by considering a single neuron that receives
diffuse synaptic noise Iknoise in the form of Poisson spike
trains from its surrounding. The capacity of recurrent
networks to produce such noise has been shown in, e.g.,
[12]. For the following analysis of individual neurons in a
noisy environment, we omit the index k and set I rec = 0.
When a conductance-based LIF neuron receives
strong synaptic stimulation, it enters a so-called highconductance state (HCS, [13]), characterized by accelerated membrane dynamics. For an analytical treatment,
it is advantageous to rewrite the membrane dynamics as
d
τeff dt
u = ueff (t) − u, such that u decays towards an effective leak potential ueff with an effective time constant
τeff = Cm /gtot [14] (see Appendix AI). The total conductance gtot subsumes both leakage and synaptic conductances. In a high input rate regime, τeff → 0 and the
effective potential ueff simply becomes a linear transformation of the synaptic noise input (see AII). Using methods similar to [15], it can be shown that, in this regime,
ueff (t) can be described as an Ornstein-Uhlenbeck (OU)

process du =
AIII):

1
τsyn (ū

− u)dt + σdW with parameters (see

P
ū = {I ext + gl El + i νi wi Eirev τsyn } / h gtot i ,
nP
o
2
2
rev
σ2 =
− ū)] τsyn / h gtot i
.
i νi [wi (Ei

(3)
(4)

The activation function of LIF neurons in the HCS

Similarly to the abstract model [7], we define the refractory state of a neuron as z(t) = 1. The mapping of
spikes to RV states (Eq. 1) naturally leads to the concept
of an activation function p(z = 1|ū), where ū = h ueff it .
In the following, we derive a general expression for the
activation function of an LIF neuron under Poisson stimulus, which is of particular use in the cortically relevant
HCS regime. Fig. 2A shows an exemplary simulation
(membrane potential and spike train) of such a scenario.
The activation function in Fig. 2B is obtained by sweeping over I ext (see AVII and AVIII).
Related setups have already been examined in [16, 17].
However, the methods in [16, 17] are tailored to certain
parameter ranges, which, in particular, do not include the
HCS regime with refractoriness. This is because [16] requires τsyn  τeff , whereas [17] assumes τref  τeff , τsyn ,
which leads to discrepancies in the predicted activation
functions (Fig. 2B). The deeper reason for the observed
discrepancies is found in the lack of an appropriate propagation of the autocorrelation of ueff through τref . Here,
we propose a derivation that explicitly includes this propagation and thereby covers a large parameter space, including the cases studied in [16, 17], as well as the HCS.
In Fig. 2A, two modes of firing can be observed: a
“bursting” mode, where the effective membrane potential ueff after the refractory period is still above threshold, and a freely evolving mode, where the neuron does
not spike again immediately after the refractory period.
This is illustrated by the distributions in Fig. 2C. Our
approach relies on the calculation of burst lengths and
their associated occurrence probability Pn .
Denoting the relative occurrence of burst lengths n by
Pn , the average drift time from % to ϑ between the kth
and (k + 1)st spike in a burst by τkb and the average
duration of the freely evolving mode that follows an nspike-burst by Tn , we identify the following relation:
P
P ·n·τ
 n n P ref
 . (5)
p(z = 1|ū) = P
n−1 b
n Pn · nτref +
k=1 τk + Tn
Knowing the OU process governing ueff , recursive expressions for Pn and Tn can be derived. These terms
have considerable impact on the activation function if
τref ≈ τsyn . In the limit of strong noise stimuli, we can
calculate the average drift time τkb in a quasistatic approximation [17] of ueff and assume that u ≈ ueff for the

3

FIG. 2. (A) Membrane potential u(t) and spikes of an LIF
neuron in a spiking noisy environment. (B) Prediction of
the activation function (red) compared to simulation results
(green), as well as to other predictions from literature [16, 17].
(C) In a HCS, u (blue) and ueff (red) are nearly identical when
the neuron is not refractory. After each refractory period
(gray), the predicted distribution of ueff (pink) is used for the
propagation in (Eq. 6). (D) High-noise regime: theoretical
prediction (red) vs. simulation results (green), with a fitted
logistic function σ(ū) (blue).

freely evolving mode, thus obtaining (see AIV):

Pn−1  R ∞
Pn = 1 − i=1 Pi ϑ dun−1 p(un−1 |un−1 ≥ ϑ)
hR
i
ϑ
× −∞ dun p(un |un−1 ) ,
(6)
R∞
Tn = ϑ dun−1 p(un−1 |un−1 ≥ ϑ)
hR
i
ϑ
× −∞ dun p(un |un < ϑ, un−1 ) h T (ϑ, un ) i , (7)


R∞
%−uk
τkb = ϑ duk ln ϑ−u
k
R∞
× ϑ duk−1 p(uk |uk ≥ ϑ, uk−1 ) .
(8)
Fig. 2C displays an intuitive picture of the integrals in
(Eq. 6) and (Eq. 7). The transfer function p(un |un−1 ) is
the Green’s function of the OU process at time t−ts = τon
and h T (ϑ, un ) i denotes the average time needed for the
membrane to reach ϑ starting from un , which can be
given in closed form [18]. The dependency of Pn , Tn and
τkb on the moments of the OU process (Eq. 3,4) renders
Eq. (5) a function of ū.
So far, we have only considered τeff = 0. To further
improve the prediction, we take into
paccount finite values
of τeff by means of an expansion in τeff /τsyn . Due to the
symmetry of the PSP shape in τeff and τsyn , this can be
done analogously to [16, 17], where the opposite limit of
large τeff and small τsyn is used (see AIV). A comparison
between our prediction of p(z = 1|ū) and results from a
numerical simulation is shown in Fig. 2B.
Sampling with networks of LIF neurons

We can now reconcile the response of LIF neurons with
the inherently stochastic neuron model (Eq. 2) that requires a logistic activation function for constant poten−1
tial v: p(z = 1|v) = σ(v) := [1 + exp(−v)] . In the HCS

FIG. 3. (A) Spike pattern of a recurrent network of K = 5 LIF
neurons during sampling from a randomly generated Boltzmann machine. (B) Sampled distribution pN (z) of network
states (blue) and target distribution pB (z) (red). pN (z) was
estimated from ten 10 s simulation runs (errorbars: std. deviation between runs). (C) DKL (pN || pB ) as a function of
integration time T for 10 trials. The red dotted line shows
convergence for the theoretically optimal abstract model [7].
(D) DKL (pN || pB ) when sampling for T = 106 ms from 100
different randomly generated target distributions.

regime, this logistic activation can be approximated by
LIF neurons with high accuracy. Fig. 2D shows our theoretical prediction and simulation results for the activation function alongside a fitted logistic function. For the
translation from the LIF domain to the abstract model
(Eq. 2) we have employed a linear mapping
v = (ū − ū0 )/α

,

(9)

with scaling factor α and ū0 denoting the potential for
which p(z = 1) = 12 .
We next connect the neurons to form a recurrent network. In addition to noise stimuli, an LIF neuron in a
network receives synaptic currents Ikrec from other neurons. For certain connectivity structures, it is possible
to predict the target distribution [7, 19] of states z (t)
that arise from the stochastic dynamics of the recurrent
network. We use the emulation of Boltzmann machines
(BMs) as an example case. The joint distribution reads:

pB (z) = Z −1 exp z T W z/2 + z T b
,
(10)
where W is a symmetric zero-diagonal weight matrix, b
is a bias vector and Z is the normalizing partition function. This probabilistic model underlies state-of-the-art
machine learning algorithms for image [20] and speech
recognition [21]. It has been shown [7] that a network of
abstract neurons (Eq. 2) with linear membrane potentials
vk = bk +

PK

j=1

Wkj zj

(11)

4

FIG. 4. Sampling in large networks. Spike raster (top left), joint distribution (bottom) and KL divergence (top right) as
in Fig. 3A,B,C. Since the target distribution cannot be computed analytically for 500 RVs, we define the Gibbs-sampling
estimate after 106 steps as a reference distribution. For the joint, the sampled distribution over 5 RVs (out of 500) obtained
from the LIF network (after Tsim = 104 ms) is plotted alongside the Gibbs estimate. Error bars indicate standard deviation
between 10 simulation runs. In all panels, biases were drawn from a beta distribution: bk ∼ 1.2 · [B(0.5, 0.5) − 0.5]. (A)
Small weights: Wkj ∼ 0.6 · [B(0.5, 0.5) − 0.5]. (B) Moderate weights: Wkj ∼ 1.2 · [B(0.5, 0.5) − 0.5]. (C) Strong weights:
Wkj ∼ 2.4 · [B(0.5, 0.5) − 0.5]. Systematic deviations between the sampled distribution and the reference manifest mainly in
low-probability modes and are due to the difference in PSP shapes between the abstract model and LIF neurons.

will sample from the desired target distribution (Eq. 10).
This finding uses the fact that individual neurons sample from the conditionals p(zk = 1 | z \k ) = σ(vk ), with
z \k = {zj | j 6= k}, in an MCMC updating scheme.
As shown above, LIF neurons in a spiking noisy environment closely approximate this logistic activation function if the synaptic currents Ikrec shift the mean membrane potential ūk according to the linear interaction
(Eq. 11). Using the linear transformation (Eq. 9) between
vk and ūk , and estimating the effect of a conductancebased synapse of weight wkj , we arrive at the following
translation between the abstract and the LIF domain (see
AV):
bk = (ūbk − ū0k )/α


rev
E
−
µ
w
kj
kj
1
Wkj =
1
1
αCm
τsyn − τeff


1 − e τeff  − ττsyn
×
−
e eff − 1 ,
e
τsyn

(12)

(13)

where ūbk is the mean potential ūk that establishes p(zk =
rev
1|ūk = ūbk ) = σ(bk ) in Eq. (5), and Ekj
denotes the
reversal potential for synapse wkj . The idea behind
(Eq. 13) is to match the integrals of individual postsynaptic potentials (PSPs) on vk and ūk . Since the membrane
loses any memory following a reset, in contrast to [7], we
use the synaptic conductance as a memory carrier.
The remaining systematic difference to the abstract
model lies in the additive – instead of renewing – nature
of PSPs elicited by the same presynaptic neuron, which

has a noticeable effect in case of fast consecutive spikes
(bursts). Since the membrane potential closely follows
the effective potential in the HCS, renewing PSPs can be
achieved by using renewing postsynaptic conductances.
For this, we have used short-term synaptic depression
[22] with a recovery time constant equal to τsyn .
The sampling quality with networks of LIF neurons
was examined in computer simulations of BMs, with randomly drawn parameters bk and Wkj (Fig. 3A,B). The
target distribution pB (z) is approximated by the distribution pN (z) of network states with τon = 10 ms. The
chosen integration time T = 10 s displays a conservative
estimate of the maximum duration a neuronal ensemble
will experience stable stimulus conditions in a behaving
organism and can thus be expected to sample from a
stable target distribution. For this integration time, we
find that the recurrent network of LIF neurons accurately
encodes the target distribution, within the precision imposed by the sample-based representation. The network
distribution pN (z) becomes increasingly more reliable as
more samples are considered (Fig. 3C). After few samples, the network has generated a coarse approximation
of pB (z) that could serve as an “educated guess” in online computation tasks. Only for simulation times T well
beyond biologically relevant timescales do systematic errors in pN (z) become apparent. The sampling quality
holds for a variety of target distributions (Fig. 3D, see
AIX for simulation details).
These observations hold for larger-scale networks as
well (Fig. 4). For weak coupling, the networks are in
an asynchronous irregular state of firing, which enables

5

FIG. 5. (A) Graphical model used for the probabilistic inference task. (B) A Gaussian likelihood model provides input to the
sampling neurons. (C) Two dimensional projection of states z (t) when sampling from the prior p(z), which was trained to
store hand-written digits 0, 3 and 4. Solid line: network trajectory over 200 ms. Color maps: Marginals of zk averaged over
25 ms. The time arrow covers the duration of the red trajectory and consecutive snapshots are 25 ms apart. (D) As in (C)
when provided with incomplete input y that is incompatible with digit 0 and ambiguous with respect to digits 3 and 4.

highly accurate sampling of the target distribution. As
synaptic weights increase, the network activity expectedly becomes more synchronous and the sampled distribution overall less accurate, especially in low-probability
regions of the state space, but the high-probability modes
are always sampled from with high fidelity.

Demonstration of probabilistic inference

We conclude our investigation of sampling in recurrent
networks of LIF neurons with an example of Bayesian
inference based on incomplete observations. A fully connected BM of K = 144 neurons, aligned on a 12×12 grid,
was trained as an associative network [23] to store patterns of hand-written digits 0, 3 and 4 in the weights Wkj
and biases bk . Each pixel of the image grid was assigned
to one network neuron. The resulting joint distribution
p(z) displays “prior knowledge” stored by the network.
The probabilistic model was augmented by adding
real-valued input channels for each pixel, associated with
random variables yk ∈ R, 1 ≤ k ≤ K. The resulting generative model p(y, z) has the structure shown
in Fig. 5A and connects the latent network variables zk
to observable inputs yk by means of likelihood functions p(yk | zk ), which we have chosen to be Gaussian
with unit variance (Fig. 5B). The likelihoods p(yk | zk )
tend to align the network state with the observation, i.e.
zk = 1 for yk > 0, while the prior p(z) reconciles the observations with knowledge on consistent activation patterns z. The task for the network is to calculate and
represent the posterior distribution according to Bayes’
rule: p(z | y) ∝ p(z) · p(y | z). A short derivation (see
AVI) shows that the posterior p(z | y) is a BM for any
input y with the following abstract membrane potential:
P
vk = bk + yk + j Wkj zj .
(14)
In the LIF domain, the sum bk + yk is equivalent to an
effective bias (Eq. 74) and corresponds to an external current Ikext = Ikb + Iky that shifts ūk appropriately. Thus,
a neuron receives synaptic input from recurrent connec-

tions and noise sources, as well as an external current,
i.e., Ik = Ikrec + Iknoise + Ikext .
In case of Iky = 0 ∀k, the network samples from
the prior distribution p(z) = p(z | y = 0). A twodimensional projection of network states z (t) ∼ p(z) is
shown in Fig. 5C. The sampled distribution has three distinct modes that correspond to the three hand-written
digits stored in the recurrent weight matrix. A closer
look at the network trajectory reveals that the system
stays in one mode (“digit”) for some duration, traverses
the state space and then samples from a different mode
of the distribution.
A typical inference scenario with incomplete observations is shown in Fig. 5D. Four input channels at the center were picked to inject positive currents Iky > 0 to the
network while all other inputs remained uninformative.
Positive currents Iky were chosen such that the observation y appeared incompatible with digit 0, and remained
ambiguous with respect to digits 3 and 4 (see AX). In
accordance with Bayes’ rule, the resulting bimodal posterior distribution p(z | y) has a suppressed 0-mode, but
preserves the 3 and 4 modes.

Discussion

We have shown how recurrent networks of
conductance-based neurons in a spiking noisy environment can perform probabilistic inference through
sampling from a well-defined posterior distribution.
Our approach extends Bayesian spiking network implementations to deterministic neuron models widely
used in computational neuroscience. We have provided
an analytical derivation of the bursty firing response of
LIF neurons under Poisson bombardment, which holds
for a wide range of parameter regimes (both high and
low ratios of τsyn /τeff , as well as both high and low
τref ). Our approach is based on the propagation of the
membrane autocorrelation throughout bursts and can
thereby provide a prediction of the activation function in
a regime where existing approaches [16, 17] do not hold

6
(see also Fig. 6 in AXI). We have further shown how
high-frequency spiking inputs that could be provided
by the surrounding network can lead to fast membrane
dynamics, which enable individual LIF neurons to correctly encode conditional distributions given information
from their presynaptic partners. Thereby, our derivation
also identifies a potential functional role of biologically
observed high-conductance states and synaptic memory
within a Bayesian framework of brain computation.
For mathematical tractability, simplifying modeling
assumptions had to be made. The neuron model uses
an absolute refractory time τref , which matches the activation time constant τon , and neglects any gradual recovery effects. On the network level, we have assumed
statistically independent noise sources and instantaneous
axonal transmission. One important difference to cortical
structure is the requirement of a symmetric connectivity
matrix. The precise symmetry is a consequence of all
neurons sharing the same parameters and can be relaxed
as neurons become diverse (α 7→ αk in Eq. 13). Furthermore, while Dale’s principle is known to not hold universally [24], negative coupling in large networks between
otherwise excitatory neurons could be, in principle, introduced through populations of inhibitory interneurons.
LIF PSPs differ from the theoretically optimal rectangular shape, which could impair convergence to the target
distribution outside of the high noise regime [25]. However, computer simulations indicate that in many biologically relevant scenarios the above approximations are
not critical (see also Fig. 7 in AXII). In particular, the
sampling properties of LIF networks remain preserved as
network size and interaction strength increase (Fig. 4).
Embedding LIF sampling in cortical-size networks appears feasible in light of our results, but remains a matter
for future work.
For neuroscientific modeling, our analysis of LIF neurons can be readily transferred to other neuron [26] and
synapse models. Beyond neuroscience, the ability to perform probabilistic inference with deterministic neurons
displays a promising computing paradigm for neuromorphic hardware systems, which typically implement physical models of integrate-and-fire neurons [27, 28]. The distributed nature of the proposed LIF sampling networks
allows to exploit the inherent parallelism of neuromorphic architectures and fosters their application to online
data evaluation and robotics. In this context, our results
have already provided the basis for the implementation
of Bayesian networks [29] and learning [30].

APPENDIX
Appendix I:
Conductance-based LIF neuron model

We restate the set of equations that govern the
conductance-based LIF model. For the membrane potential, we have
Cm

du
= gl (El − u) + I
dt

,

(15)

with capacitance Cm , membrane potential uk , leak potential El , leak conductance gl and input current I. The
input current I can be formally partitioned into recurrent synaptic inputs, diffuse synaptic noise and additional “external” currents, i.e.,
I = I rec + I noise + I ext

.

(16)

I ext may represent external current stimuli, average
synaptic stimulus currents, or changes in the leak mechanism (i.e., changes in gl or El ). The total synaptic current
I syn = I rec + I noise obeys the equation
X syn
I syn =
gi (Eirev − u) ,
(17)
syn i

where gisyn
and Eirev

represents the conductance at the ith synapse
the corresponding reversal potential. The
synaptic conductance obeys the ODE
X
g syn
dgisyn
=− i +
wi δ(t − ts ) ,
dt
τsyn

(18)

spk s

with synaptic time constant τsyn and weight wi . The sum
runs over all presynaptic spikes s. The solution to this
equation is a superposition of exponentials:


X
t − ts
syn
gi =
wi Θ(t − ts ) exp −
.
(19)
τsyn
spk s

Putting all of the above together, we obtain the full ODE
for the membrane potential
Cm

du
= gl (El − u)
dt


XX
t − ts
(Eirev − u)
+
wi Θ(t − ts ) exp −
τ
syn
i
spk s

+ I ext

.

(20)

We can now divide the RHS of (Eq. 20) by gtot = gl +
P
syn
and rearrange the terms in order to obtain
i gi
τeff

du
= ueff − u
dt

,

(21)

with a new effective membrane time constant
τeff =

Cm
gtot

(22)

7
and

and effective leak potential

ueff (t) =

gl El +

syn
rev
i gi (t)Ei

P

+ I ext

gtot (t)

Var [gisyn ]
.

(23)

This transformation is routinely used in studies of
conductance-based neurons (see, e.g., [14]). Here, we
have made the time dependencies explicit, since the following HCS approximation will serve to eliminate t in
the denominator.




t − ts
=
Var wi Θ(t − ts ) exp −
τsyn
spk s
(* 
2 +

t − ts
= lim h N i
wi Θ(t − ts ) exp −
T →∞
τsyn
)


 2
t − ts
+
wi Θ(t − ts ) exp −
τsyn
(


Z T
t
1 2
exp −2
= lim νi T
w
dt
T →∞
T i 0
τsyn
"Z

 #)
T
t
1
exp −
− 2
dt
T
τsyn
0
X

=
Appendix II:
The high-conductance state (HCS)

In a first approximation, assuming a rapidly-firing
Poisson background (νsyn → ∞), the total average conductance can become arbitrarily large (h gtot i → ∞),
causing the membrane potential to follow the effective
potential nearly instantaneously (h τeff i → 0). Eq. (23)
can then be rewritten as

u ≈ ueff =

gl El +

syn
i h gi

P

P
iEirev + i ∆gisyn Eirev +I ext
P
,
h gtot i + i ∆gisyn
(24)

1 2
w νi τsyn
2 i

,

(27)

so the relative fluctuations of g syn are of the order
s
p
Var [gisyn ]
1
(28)
=
h gisyn i
2 νi τsyn
and vanish in the limit of large firing rates.
This warrants an expansion of Eq. (24) in ∆gisyn , ∀i.
Considering only the first-order term we obtain
P
I ext + gl El + i gisyn (t)Eirev
u(t) =
,
(29)
h gtot i
which renders u simply a linear
P transformation of the
synaptic noise current J syn = i gisyn Eirev .

where
Appendix III:
Derivation of the equivalence to an OU process

∆gisyn = gisyn − h gisyn i

(25)
From Eq. (18), we can find that the synaptic noise
J syn obeys the first-order inhomogenous ODE

denotes the fluctuations of the synaptic conductances.
For a single Poisson source with rate νi connected to
the neuron by a synapse with weight wi and time constant τsyn , the conductance course can be seen as a sum
of independent random variables, each of them representing the conductance change caused by a single spike. In
the limit of large νi , the central limit theorem guarantees the convergence of the conductance distribution to
a Gaussian, with moments given by

h gisyn


t − ts
i=
wi Θ(t − ts ) exp −
τsyn
spk s


Z T
hN i
t
exp −
= lim
wi
dt
T →∞ T
τsyn
0
X

= wi νi τsyn



J syn X X
dJ syn
=−
+
∆Jisyn δ(t − ts ) ,
dt
τsyn syn i

(30)

spk s

where ∆Jisyn = wi Eirev . This equation is highly reminiscent of the ODE that defines the OU process
dx(t) = θ [µ − x(t)] dt + σdW (t) .

(31)

It is well-known that the PDF of the OU process
s
θ
f (x, t|x0 ) =
πσ 2 (1 − e−2θt )



−θ (x − µ + (µ − x0 )e−θt )2
× exp
(32)
σ2
1 − e−2θt
is the unique solution of the Fokker-Planck equation

(26)

∂
σ2 ∂ 2 f
1 ∂f (x, t)
=
[(x − µ)f ] +
θ ∂t
∂x
2θ ∂x2

(33)

8
with starting condition x0 := x(t = 0). In the following,
we prove that, under certain assumptions, the distribution of the synaptic input J syn obeys the same FokkerPlanck equation. To this end, we follow an approach
similar to [15].
Consider the PDF of the synaptic input f (J syn , t). We
can use the Chapman-Kolmogorov equation to describe
its evolution after a short time interval ∆t as an integral
over all possible intermediate states J 0 :

f (J syn , t + ∆t) =

Z

∞

f (J syn , t + ∆t|J 0 , t)f (J 0 , t)dJ 0

−∞

(34)
For a small enough ∆t, the probability of the occurrence
of multiple spikes within ∆t can be neglected. As incoming spikes are assumed to be generated by Poisson
processes,Pthe probability of a single spike occurring in
∆t is ∆t i νi . By summing over the two possible histories of J syn within ∆t (either a single incoming spike or
no spike at all), we can use Eq. (30) to find

and rearrange the terms to obtain
f (J syn , t + ∆t) − f (J syn , t)
∆t
∂f (J syn , t + ∆t)
=
∂∆t
∆t=0
(

 

 
X
∆t
∆t
syn
= −
νi exp
f J exp
,t
τsyn
τsyn
i
!
X
1
+ 1 − ∆t
νi
τ
syn
i
(

 

 
∆t
∆t
syn
× exp
f J exp
,t
τsyn
τsyn
 i)
h

∆t


syn
∂f
J
exp
τsyn , t
∆t


+ exp 2
J syn
τsyn
∂J syn exp τ∆t
syn

 



X
∆t
∆t
+
νi exp
f J syn exp
− ∆Jisyn , t
τsyn
τsyn
i
)
+ (. . . )∆t

(38)
∆t=0

f (J syn , t + ∆t|J 0 )
"
# 


X
∆t
= 1 − ∆t
νi δ J syn − J 0 exp −
τsyn
i



X
∆t
, (35)
+ ∆t
νi δ J syn − (J 0 + ∆Jisyn ) exp −
τsyn
i

By taking the limit ∆t → 0, we obtain:
1
∂
∂f (J syn , t)
=
[J syn f (J syn , t)]
∂t
τsyn ∂J syn
X
+
νi [f (J syn − ∆Jisyn , t) − f (J syn , t)]

.

i

where νi represents the afferent firing frequency at the
ith synapse. Plugging this into Eq. (34) and integrating
over J 0 yields

f (J syn , t + ∆t)
!

 

 
∆t
∆t
f J syn exp
,t
τsyn
τsyn
i

 



X
∆t
∆t
+ ∆t
νi exp
f J syn exp
− ∆Jisyn , t .
τsyn
τsyn
i
=

1 − ∆t

X

νi

(39)
In the limit of small synaptic weights (i.e., ∆Jisyn → 0),
we can expand the second term on the RHS up to the
second order in ∆Jisyn . This yields, after some rearrangement



exp

(36)

∂f (J syn , t)
∂t
"
!
#
X
∂
1
syn
syn
syn
=
J
−
νi ∆Ji τsyn f (J , t)
τsyn ∂J syn
i
P
2
νi ∆Jisyn ∂ 2 f (J syn , t)
+ i
,
(40)
2
∂J syn 2

We can now expand f (x, t + ∆t) up to first order in ∆t

f (J syn , t + ∆t) ≈ f (J syn , t) +

∂f (J syn , t + ∆t)
∂∆t

∆t
∆t=0

(37)

which is the exact equivalent of the Fokker-Planck equation of the OU process (Eq. 33). Since u(t) is only a
linear transformation of J syn (t), it can also be approximated by an OU process in the limit of large input frequencies and small synaptic weights, with Eq. (29) and
Eq. (40) giving the specific time constant, mean value

9
and variance:
1
,
τsyn
P
I ext + gl El + i νi wi Eirev τsyn
µ=
,
h gtot i
P
2
νi [wi (Eirev − µ)] τsyn
σ2
= i
.
2
2
2 h gtot i
θ=

(41)
(42)
(43)

We conclude this section with two important notes.
Firstly, for the above methodology to be generally appliccable, we must be able to take the limit ∆Jisyn → 0 for
arbitrary first and second moments of f (J syn , t) without
modifying them. This is possible if at least one excitatory
and one inhibitory input is present, which then give us
two degrees of freedom with a proper choice of ν exc → ∞
and ν inh → ∞. Secondly, all higher moments (3 and
above) need to vanish in the abovementioned limit. This
has been shown to also be the case under the above conditions [31].
Appendix IV:
Derivation of the activation function

We can distinguish between two firing modes of the
neuron. The first mode can be classified as “burst spiking” and occurs when multiple spikes occur in rapid
succession with an expected ISI of h ∆tk i = τref +
ts+1
− tsk s = τref + τkb , where τkb represents the avk
erage drift time from the reset to the threshold potential
following the kth refractory period within a burst. In
this case, for each spike within a burst,
ueff (ts ) ≥ ϑ

,

(44)

and also, for all but the last spike,
ueff (ts + τref ) ≥ ϑ

.

(45)

The second mode appears between such bursts, where the
membrane potential evolves freely in the subthreshold
regime. If we define, just like in the abstract model, that
the kth neuron is in the state zk = 1 for a duration
τon = τref following a spike, we can write
P
P ·n·τ
 n n P ref
 , (46)
p(z = 1) = P
n−1 b
P
·
nτ
+
τ
+
T
n
ref
n
n
k=1 k
where Pn represents the distribution of burst lengths
(conditioned on the existence of the first spike) and Tn is
the mean time interval between the end of a burst (i.e.,
the endpoint of its last refractory period) and the next
spike. The variables Pn , Tn and τkb depend on all neuron
and noise parameters, but for calculating the activation
function (Fig. 2 in the main manuscript and Fig. 6 in the
Appendix), we only vary ū.

We can now calculate both Pn and Tn iteratively. The
idea behind this approach is to propagate the membrane
potential PDF from spike to spike within a burst and cut
off the irrelevant parts for a particular burst length n.
We denote the spike times within a burst of length n
by t0 , . . . , tn−1 and the endpoint of such a burst by tn :=
tn−1 +τref . For brevity, we also use ui := u(ti ). Assuming
a first spike at some time t0 (u0 := u(t0 ) = ϑ), a “burst”
of length n = 1 requires a subthreshold free membrane
potential after the first refractory period (u1 := u(t0 +
τref ) < ϑ). This occurs with probability
P1 : = p(u1 < ϑ|u0 = ϑ)
Z ϑ
=
du1 p(u1 |u0 = ϑ),
−∞
{z
}
|

(47)

I1

where p(ui+1 |ui ) := f (u, τref |ui ), which was defined in
Eq. (32). On average, the neuron then stays in the subthreshold regime for a period equal to the mean first passage time from u1 to ϑ, so the mean duration of the time
interval until the onset of the next burst can be expressed
as
Z ϑ
T1 =
du1 p(u1 |u0 = ϑ) h T (ϑ, u1 ) i .
(48)
−∞

The first-passage time problem of the OU process has
often been discussed [32]. While no closed-form expression for the distribution of first-passage times T (b, a) =
inf t ≥ 0 : x(t) = b|x(0) = a is known, its moments can be
computed analytically [18]. In particular, the mean first
passage time reads
h T (b, a) i
r Z




θ π b
(x − µ)2
x−µ
=
dx exp
. (49)
1+erf √
σ 2 a
2σ 2
2σ
A burst of n = 2 spikes can only occur when the effective membrane potential lies above the spiking threshold
(u1 ≥ ϑ) after the first refractory period and below after
the second (u2 < ϑ). This makes P2 and T2 recursive
functions of P1 :
P2 = p(u2 < ϑ, u1 ≥ ϑ|u0 = ϑ)
= p(u1 ≥ ϑ|u0 = ϑ) p(u2 < ϑ|u1 ≥ ϑ, u0 = ϑ)
"Z
#
Z ∞
ϑ
(Eq.47)
= (1 − I1 )
du1 p(u1 |u1 ≥ ϑ)
du2 p(u2 |u1 )
| {z } ϑ
−∞
=1−P1 |
{z
}
I2

(50)
Z

∞

du1 p(u1 |u1 ≥ ϑ)

T2 =
ϑ

"Z

ϑ

×

#
du2 p(u2 |u2 > ϑ, u1 ) h T (u2 , ϑ) i ,

−∞

(51)

10
where p(ui |ui ≥ ϑ) is a shorthand notation for p(ui |ui ≥
ϑ, ui−1 ≥ ϑ, . . . , u1 ≥ ϑ, u0 = ϑ). In particular, this
represents a renormalization of the PDF of the effective
membrane potential to values above the spiking threshold
after i refractory periods.
We can now continue this recursion up to an arbitrary
burst length and write
Pn = p(un < ϑ, un−1 ≥ ϑ, . . . , u1 ≥ ϑ|u0 = ϑ)
= p(u1 ≥ ϑ|u0 = ϑ)
× p(un < ϑ, un−1 ≥ ϑ, . . . , u2 ≥ ϑ|u1 ≥ ϑ, u0 = ϑ)
(Eq.47)

=

(1 − I1 ) p(u2 ≥ ϑ|u1 ≥ ϑ, u0 = ϑ)

× p(un < ϑ, un−1 ≥ ϑ, . . . , u3 ≥ ϑ|
|u2 ≥ ϑ, u1 ≥ ϑ, u0 = ϑ)
(Eq.50)

=

(1 − I1 )(1 − I2 ) p(u3 ≥ ϑ|u2 ≥ ϑ, u1 ≥ ϑ, u0 = ϑ)

× p(un < ϑ, un−1 ≥ ϑ, . . . , u4 ≥ ϑ|
|u3 ≥ ϑ, . . . , u1 ≥ ϑ, u0 = ϑ)
=

(1 − Ii ) p(un < ϑ|un−1 ≥ ϑ, . . . , u1 ≥ ϑ, u0 = ϑ)

i=1

(52)
n−1
X

= (1 −

Pi )

i=1

"Z

∞

×

#

ϑ

dun−1 p(un−1 |un−1 ≥ ϑ)

dun p(un |un−1 )
−∞

ϑ

{z

|

}

In

(53)
Z

∞

dun−1 p(un−1 |un−1 ≥ ϑ)

Tn =
ϑ

"Z

τkb (uk )

#

ϑ

×

dun p(un |un < ϑ, un−1 ) h T (un , ϑ) i

(54)

−∞

The transition from a product to a sum between Eq. (52)
and Eq. (53) requires the identity
n−1
Y

n−1
X

i=1

i=1

(1 − Ii ) = 1 −

Pi

,

(55)

which
can be easily shown by induction from Pn =
Qn−1
In i=1 (1 − Ii ) (Eq. 52) and P1 = I1 (Eq. 47). Since
limn→∞ Pn = 0, one can stop the recursion at some small
enough Pn .
What remains to be calculated is the average time-to
threshold τkb within a burst that follows the kth refractory period. Since we assume a HCS, we are looking at
a regime in which τeff  τsyn . Therefore, we can assume
ueff to be approximately unchanged during the short time
interval τkb (adiabatic approximation, see also [17]). For
a fixed uk , the jump time can be easily calculated from


= ln

% − uk
ϑ − uk


.

(56)

The average jump time can then be obtained by integrating over all suprathreshold values of uk , which in
turn have probabilities that follow from integrating over
all suprathreshold values of uk−1 :

Z ∞
Z ∞
% − uk
duk−1 p(uk |uk > ϑ, uk−1 ) .
τkb =
duk ln
ϑ − uk
ϑ
ϑ
(57)
With Eq. (46), (53), (54), (57) and (49), one could
now predict the activation function of an LIF unit in an
extreme high-noise regime (τeff → 0). We can, however,
generalize our approach by taking the finite nature of the
effective time constant into account.
If we go back to Eq. (21) and leave τeff = C/ h gtot i
small but finite, we can still perform all the remaining
approximations, but are required to modify Eq. (29):
τeff u̇(t) =

n−1
Y

Z

Eq. (21):

J syn (t)
I ext + gl El
+
− u(t) .
h gtot i
h gtot i

(58)

Together with Eq. (30), we now have a system of firstorder ODEs which can be solved analytically by standard
techniques (variation of constants). The PSPs are then
no longer a linear transformation of the exponentially
shaped PSCs, but rather alpha-shaped (more precisely, a
difference of exponentials):
 t−ts

− t−ts
−
e τeff − e τsyn
us (t) = Θ(t − ts )A
,
(59)
τeff − τsyn
w (E rev −h u

i)τ

with A = i i h gtoteffi syn i . This shape causes a lower
PSP peak than in the case of exponential PSPs, decreasing the overall width of the membrane potential distribution. Intuitively speaking, this results in a horizontal
shift and compression of the activation function.
More recently, analytical treatments of these phenomena have been proposed [33]. In these approaches, large
membrane time constants (equivalent to a long τeff ) and
small synaptic time constants are usually considered.
However, Eq. (59) is symmetric in τeff and τsyn , so the
same argument applies to our case as well, but the two
time constants need to be switched. It is, for example,
possible to correct the first passage time from the reset
p to
the threshold potential by using an expansion in τ 0 /τ
(with τ 0 and τ being the smaller and the larger of the
two time constants, respectively) [16]:
Z ϑeffσ−µ
√
h T (ϑ, u) i = τ π
dx exp(x2 )[erf(x)+1] , (60)
u−µ
σ

with µ and σ 2 the first two moments of the free membrane
potential distribution and an effective threshold
 r 0
1
τ
ϑeff ≈ ϑ − ζ
σ ,
(61)
2
2τ

11
in which ζ denotes the Riemann zeta function.
In our
p
particular case, the expansion is done in τeff /τsyn , so
τ 0 = τeff and τ = τsyn . With this approximation, we
assume that u converges from ρ to ueff in negligible time
after it is released from the refractory state. Afterwards,
its convergence to ueff is determined by Eq. (21). Note
how Eq. (60) is equivalent to a change of the integration
variable and limits in the original equation (Eq. 49) for
the first passage time.
Appendix V:
Translation of synaptic weights

A sufficient condition for a single neuron to sample
from the correct conditional distribution is given by its
activation function:
p(zk = 1|z\k ) = σ(vk ) .

(62)

The relationship between the abstract model and the LIF
implementation is defined by the lateral dilation α and
relative offset ū0k of the LIF activation function:


ūk − ū0k
.
(63)
p(zk = 1|z\k ) = σ
α
The parameters α and ū0k can be determined by fitting
Eq. (63) either to simulation results or to the theoretical prediction (Eq. 46). As a consequence, also synaptic
weights need to be rescaled by the factor α. Additionally, the difference in PSP shapes needs to be taken into
account.
We choose a translation rule in which we set the LIF
synaptic weights wij such that the area under a PSP
(Eq. 59) during the refractory state of the corresponding afferent neuron (i.e., for a duration τref ) is equal to
Wij τref α:
Z τref
rev
wkj (Ekj
− h ueff i)τsyn
Wkj τref α =
h gtot i
0




t−ts
s
exp − τeff − exp − t−t
τsyn
×
dt
τeff − τsyn


rev
wkj τsyn Ekj − µ
=
h gtot i τeff − τsyn
h
 τref
i

 τref
−
−
× τsyn e τsyn − 1 − τeff e τeff − 1

Additionally, depressing short-term plasticity [22] has
been applied to attenuate the amplitudes of consecutively arriving alpha-shaped PSPs from a network neuron
and emulate renewing synapses. In particular, within the
Tsodyks-Markram short-term plasticity model [22], the
synaptic efficacy parameter and recovery time constant
have been chosen as USE = 1 and τrec = τsyn , respectively.
Analogously to the weights, the biases
can be deter
mined from the condition σ(bk ) = σ( ūbk − ū0k /α) in the
absence of recurrent activity, i.e. Ikrec = 0.
Appendix VI:
Probabilistic model for demonstration of inference

We define a joint model p(y, z) over real-valued input nodes y = (y1 , . . . , yK ) and binary latent variables
z = (z1 , . . . , zK ) as sketched in Fig. 4A of the main
manuscript. The real-valued variables yk encode intensities of the input pixels. The latent variables zk correspond to neurons in the network. Our aim is to demonstrate that the network can sample from the posterior
distribution p(z | y).
For this example, we have chosen a particularly simple
likelihood function, namely a Gaussian emission model
with variance σ 2 = 1 and mean values µ = ± 21 . From
the graphical model we identify the structure of the joint
distribution
p(y, z) = p(z) ·

K
Y

p(yk | zk )

(66)

k=1

where p(z) is a Boltzmann distribution and the likelihood
is defined by
p(yk | zk ) = N (yk ; µ = 1/2)zk · N (yk ; µ = −1/2)1−zk ,
(67)
which is equivalent to
log p(yk | zk ) = zk · [log N (yk ; 1/2) − log N (yk ; −1/2)]
+ log N (yk ; −1/2)
= zk · yk + log N (yk ; −1/2) ,

(68)

using normal distributions
N (y; µ) := exp[−(y − µ)2 /2] /

√

2π

(69)

. (64)
with unit variance. The posterior of this model reads:

By setting τref = τsyn , we obtain the mapping between
the abstract and LIF synaptic weight domains:


rev
w
E
−
µ
kj
kj
1
Wkj =
τ
αCm
1 − τsyn
eff
h
 τsyn
i

−
−1
× τsyn e − 1 − τeff e τeff − 1
. (65)

p(y, z)
p(y)

.
1 T
= exp
z W z + z T (b + y)
Norm , (70)
2

p(z | y) =

where the normalizing constant depends on the input y,
but is independent of the network variables and thus

12
defines a Boltzmann distribution over z for any y. In
particular, we identify
the abstract membrane potential
P
vk = bk + yk + j Wkj zj for sampling from p(z | y) by
means of a spiking network.

N

Appendix VII:
Simulation parameters

0.1 nF
5 nS
-65 mV
-53 mV
0 mV
-90 mV
-52 mV
10 ms
10 ms

with pi =

Nispk τon
being the fraction
Tsim
spk
Ni representing the total

,

(73)

of time spent in

Appendix IX:
Sampling via recurrent networks of LIF neurons
(Fig. 3 and 4 in main manuscript)

TABLE I. Neuron parameters used for the simulations in the
main manuscript.

Synaptic noise was implemented as bombardment by
inhibitory and excitatory Poisson stimuli with rates
νinh = νexc = 5000 Hz. The excitatory synaptic weight
noise
= 0.0035 µS.
for the noise stimuli was set to wexc
noise
The inhibitory weight winh was adjusted as to yield
p(zk = 1) ≈ 0.5 with no current stimulus present. For
above parameters, this happens at an average free memnoise
brane potential of ū = −55 mV. This determines winh
according to
.

(72)

z = 1 and
number of spikes
in the ith out of N = 10 performed simulations. Since
the respective standard errors of the mean are smaller
than the size of the symbols, no error bars are shown.

membrane capacitance
leak conductance
leak potential
reset potential
excitatory reversal potential
inhibitory reversal potential
threshold voltage
synaptic time constant
refractory time constant

rev
Einh
− ū
wnoise
= exc
noise
rev
Eexc − ū
winh

1 X
pi ,
N i=1
v
u
N
u 1
X
2
s=t
·
[pi − p(z = 1)]
N − 1 i=1

p(z = 1) =

All simulations have been performed with the NEURON simulation package [34] and the PyNN API [35],
with a time step of dt = 0.01 ms. For the LIF neuron,
we have chosen the following parameters (compare with,
e.g., [36] for parameters fitted to experimental data):
Cm
gl
El
ρ
rev
Eexc
rev
Einh
ϑ
τsyn
τref

The deviations from the theoretical prediction (Eq. 46)
are smaller than the size of the symbols, therefore no
errorbars are shown.
The ordinate values and standard errors were calculated from the simulated spike train data according to

(71)

Appendix VIII:
The activation function of LIF neurons in a spiking
noisy environment (Fig. 2 in main manuscript)

In order to sweep through the activation function, the
external current I ext was varied. However, in order to
facilitate a comparison with the logistic activation function of the abstract model, we have represented p(z = 1)
as a function of ū instead. The latter is equivalent to the
mean µ of the corresponding Ornstein-Uhlenbeck process, with Eq. (42) allowing a direct translation between
I ext and ū.
The abscissa values in Fig. 2D in the main manuscript
represent averages of the free membrane potential obtained from 10 simulation runs with a total duration of
rev
Tsim = 100 s and firing threshold θ set to Eexc
= 0 mV.

The simulated network consists of K = 5 neurons
with a synaptic weight matrix W and a bias vector
b (both in the Boltzmann domain). All entries were
drawn from a beta distribution B(0.5, 0.5) and mapped
linearly to the interval [−0.6, 0.6]. More specifically,
bk , Wkj ∼ 1.2 · [B(0.5, 0.5) − 0.5]. The parameters and
mapping of the beta distribution were chosen with the
intent of generating diverse distributions, spanning multiple orders of magnitude. The bias bk , defined in the
Boltzmann domain, determines the probability p(zk =
1 | z \k = 0) for neuron k. In the LIF domain, the probability p(zk = 1 | z \k = 0) = 0.5 corresponds to the mean
free membrane potential ū0k . Then, a nonzero bias can
be described in the LIF domain as a linear shift from ū0k
to a mean membrane potential ūbk . This yields the linear
transformation
bk =(ūbk − ū0k )/α

,

(74)

where α represents the scaling factor between the two
domains. Both quantities ū0k and α can be determined
from the predicted activation function of a single LIF
unit. The first quantity constitutes the inflection point
of the activation function (at p(zk = 1 | z \k = 0) = 0.5),
the latter follows from the slope of the function.
By computing ūbk , we can map any bias bk of a single
unit of the Boltzmann machine onto a yet unconnected
LIF neuron. In simulations, ūbk was established by injectk
ing a temporally constant external current Iext
according
to
X
Ikext = (αbk +ū0k ) h gtot i−gl El − νi winoise Eirev τsyn . (75)
i

13
In order to achieve sampling network dynamics in the
LIF domain faithful to those displayed by an equivalent
Boltzmann machine, the Boltzmann weight matrix W
was translated into LIF network weights wij according
to Eq. (65). Thus, superposing PSPs saturate the membrane potential, approximating the constant amplitude
of a PSP in the abstract neuron model.
For Fig. 3B in the main manuscript, this setup of a
random Boltzmann machine was simulated N = 10 times
with different random seeds for the Poisson background
for a duration of Tsim = 10 s. The red bars show the
analytically computed target joint distribution pB (z).
The blue bars depict the network distribution pN (z),
calculated from the firing activity of the simulated LIF
network set up to match pB (z). The means and error
bars have been calculated as in Eq. (72) and (73),
respectively.
The above simulations were repeated with a significantly longer duration in order to study systematic deviations due to the LIF implementation. Fig. 3C in the main
manuscript shows the distance between the target distribution pB (z) and its LIF network representation pN (z)
in form of the Kullback-Leibler divergence
DKL (pN || pB ) = h log [ pN (z) / pB (z) ] ipN (z)

(76)

This estimate has been taken for one set of parameters
(W , b) for ten independent trials (thin lines) in an LIF
network at integration times T : 0 ≤ T ≤ Tsim = 106 ms.
The red dashed line displays the averaged DKL (pN || pB )
for the abstract network model with identical parameters
(W , b). The decrease of DKL (pN || pB ) for longer integration times indicates the increasing precision of the sampling network over time. Eventually, the DKL (pN || pB )
converges for the LIF network to a nonzero value, reflecting small systematic errors. Fig. 3D in the main
manuscript shows the distribution of DKL (pN || pB ) values for 100 randomly drawn Boltzmann machines emulated by LIF networks, evaluated from a single run of
Tsim = 106 ms each.
Fig. 4 was generated in the same way as Fig. 3, except
for the network size and weight distribution.
Appendix X:
Demostration of probabilistic inference
(Fig. 5 in main manuscript)

We trained a fully visible, fully connected Boltzmann
machine to store three hand-written digits (0, 3, 4) that
were taken from the MNIST data set [37] and were scaled
down to 12x12 pixels. The pixel intensities of these patterns (ranging from 0 to 1) were linearly scaled to an activation between 0.05 and 0.95, defining the target statistics h zk iT and h zk zj iT for the Boltzmann machine.
Then, external currents Ikb and synaptic weights wkj

were optimized via the update rules ∆Ikb ∝ h zk iT − zk
and ∆wkj ∝ h zk zj iT − zk zj , with samples z obtained
through sampling from an LIF network set up with
synaptic noise and neuron parameters as described above.
The recurrent connections, defined by wkj , induce
synaptic currents Ikrec in addition to the noise currents
Iknoise . Additionally, each neuron’s mean effective membrane potential Eq. (42) is shifted by external currents
of the trained quantities Ikb as well as input currents Iky .
The latter encode observations y which are defined in
the context of the probabilistic model described in AVI.
Hence, the total received current of neuron k in the network amounts to
Ik = Ikrec + Iknoise + Ikb + Iky .

(77)

The resulting samples z (t) displayed in subplots 5C
and 5D were taken for 4000 ms, after a burn-in time
that ensured that the network had converged to its
equilibrium distribution.
Figure 5C, main manuscript: Sampling from the prior.
From the probabilistic model it follows that p(z | y =
0) = p(z). This means that the LIF network will sample
from the prior when Iky = 0, ∀k. The figure makes use
of two projections of network states z (t) to illustrate the
sampled distribution:
1. Star plot: In order to illustrate the 12x12dimensional network states z (t) , a two-dimensional
linear projection in a star plot has been chosen
(blue dots). The axes indicate the three basis vectors B, representing pixel intensities h zk iT of the
digits (0, 3, 4),
034

h zk iT

= (B 0 , B 3 , B 4 )T

(78)

i
with
qP a total intensity normalization ||B || =
i 2
j |Bj | = 1.

The network states z (t) acquired from the simulation are projected onto this basis:
z 034 (t) = (B 0 · z (t) , B 3 · z (t) , B 4 · z (t) )T

.

(79)

This three-dimensional vector is then projected
onto a two-dimensional plane in coordinates


sin(φ0B ) sin(φ3B ) sin(φ4B )
proj
z
(t) =
z 034 (t) (80)
cos(φ0B ) cos(φ3B ) cos(φ4B )
with (φ0B , φ3B , φ4B ) = (0, 23π , 43π ) indicating the directions of the normalized basis vectors.
These linear projections z proj (t) of network states
z (t) are used to illustrate similarity of states as distance in the two-dimensional plane. In Fig. 5C in
the main manuscript, a network evolution time of

14
4000 ms is shown, samples z (t) taken every 2 ms,
i.e. in total 2000 projected states are displayed. We
ensured in longer simulations that the total simulation runtime is sufficient to represent the distribution under the mixing of the Markov chain.
The significant clustering of the network states
around the directions of the arrows indicates the
proximity to the target states (0, 3, 4) for a majority of time.
The transitions of the network states z (t) are depicted as a red trajectory in the star plot. This
trajectory connects 100 projected network states
within a time interval of 200 ms, demonstrating the
time evolution of z proj (t).
2. Snapshots: The squared color maps on the time
axis display the averaged pixel intensity of the network in a time window 25 ms:
Z t+12.5 ms
1
(t0 )
zk dt0 .
(81)
z̄k (t) =
25 ms t−12.5 ms
The equidistant time intervals between these snapshots were taken every 25 ms within the time frame
of the above mentioned red trajectory of 200 ms.
In summary, the network spends most of the time in
distinct modes (digits) and only little time in “blurred”
states that could not be clearly assigned to one digit.
Furthermore, it spends an approximately equal amount
of time in each mode, indicating that the prior p(z) does
not favor one of the three digits.
Figure 5D, main manuscript: Sampling from the posterior. Incomplete and ambiguous input y was provided
to the network by setting four of the inputs different from
zero: yk 6= 0 for k ∈ I with I denoting an index set and
|I| = 4. These inputs were chosen at the center of the
image where digits 3 and 4 both have black pixels, while
digit 0 is white. More precisely, I = {77, 78, 79, 80} when
pixels are indexed row-wise starting from the top-left corner. In the LIF implementation, positive values yk correspond to positive currents Iky , as specified in Eq. (75).
We set Iky = 0.831 nA for the four non-zero inputs, inducing an effective bias bk + yk by injecting a total current
Ikb + Iky . For the plot, the same projections of network
states were used as in Fig. 5C in the main manuscript.
Under the incomplete input, the network spends only
little time in the 0-mode or “blurred” states, while the
equilibrium distribution exhibits two distict modes in the
3- and 4-directions. Thus the posterior reflects both the
almost certain conclusion that “the input is not a zero”
and the uncertainty that “the input could either be a
three or a four”. In particular, the network response is
well-suited for further processing (e.g. by other cortical
populations or in a technical application). For instance,

the network states could be integrated by a linear classifier to recognize the digit class.

Appendix XI:
Prediction of the activation function for various
parameter sets

We have used Eq. (46) for predicting the activation
functions of LIF neurons in the HCS regime. In this
regime, the refractory time τref and the synaptic time
constant τsyn become the dominant time constants:
τeff  τref ≈ τsyn

.

(82)

So far, we have compared our prediction to simulation
data, as well as to two other predictions by [16] and [17]
in Fig. 2B of the main manuscript. Here, we depart from
these assumptions and show in Fig. 6 that our prediction
holds for several different parameter sets. In particular,
the predictions from [16] and [17] are only valid when either synaptic time constants (Fig. 6C) or refractory times
(Fig. 6E) become shorter.

Appendix XII:
Validity of the LIF sampling framework for various
parameter sets

For our various network simulations (Fig. 3-5 in the
main manuscript) we have used a set of biologically plausible parameters from [36] (see also “Appendix VII: Simulation parameters”). In particular, we have used high
input noise rates and weights in order to achieve the HCS,
as well as a small distance between firing threshold and
reset, in order to reduce τkb , during which a neuron falsely
encodes the state z = 0. Furthermore, we have assumed
τref = τsyn . Fig. 7 shows that these are by no means
strict constraints.

∗

[1]
[2]
[3]
[4]
[5]

The first two authors have contributed equally to this
work. We thank W. Maass for his essential support,
as well as S. Habenschuss, M. Brixner and the three
anonymous reviewers for their helpful comments. This
research was supported by EU grants #269921 (BrainScaleS), #237955 (FACETS-ITN), #604102 (Human
Brain Project), the Austrian Science Fund FWF #I753N23 (PNEUMA) and the Manfred Stärk Foundation.
K. Körding and D. Wolpert, Nature 427, 244 (2004).
J. Fiser, P. Berkes, G. Orbán, and M. Lengyel, Trends
in cognitive sciences 14, 119 (2010).
K. Friston, J. Mattout, and J. Kilner, Biological cybernetics 104, 137 (2011).
T. Yang and M. N. Shadlen, Nature 447, 1075 (2007).
P. Berkes, G. Orbán, M. Lengyel, and J. Fiser, Science
331, 83 (2011).

15

FIG. 6. Comparison of our prediction of the activation function (Eq. 46) to simulation data, as well as to the predictions
given by [16] and [17], for several parameter sets. (A) Standard parameter set as given in Tab. I. This figure shows the same
curves as Fig. 2B in the main manuscript. (B) Same as A, but with quadrupled membrane capacitance Cm and one quarter
of the leak conductance gl , input rates νexc,inh and weights wexc,inh . This parameter set is identical to the one used for the
top left corner of Fig. 7 and has the effect of slowing the membrane, i.e., increasing τeff by a factor of 16. (C) Same as A, but
with a decreased synaptic time constant τsyn = 3 ms. The prediction from [16] is improved, since the correlations in the preand post-refractory effective membrane potential are smaller in this scenario. (D) Same as A, but with an increased synaptic
time constant τsyn = 30 ms. The prediction from [16] deteriorates due to the longer-range membrane potential autocorrelation.
Conversely, the prediction from [17] improves, since the refractory time becomes less important. (E) Same as A, but with a very
short refractory time τref = 1 ms. Here, we enter the parameter range where [17] provides good predictions. (F) Same as A,
but with the input rates νexc,inh and weights wexc,inh decreased by a factor of 10, thereby slowing the membrane considerably
(imperfect HCS). Additionally, we have chosen a large reset-to-threshold distance of ϑ − % = 10 mV. In this scenario, the
τkb -term in Eq. (46) becomes dominant and the activation function departs from the logistic shape that it has in the HCS.

[6] S. Deneve, Neural computation 20, 91 (2008).
[7] L. Buesing, J. Bill, B. Nessler, and W. Maass, PLoS
computational biology 7, e1002211 (2011).
[8] R. P. Rao (2004) pp. 1113–1120.
[9] P. O. Hoyer and A. Hyvarinen, Advances in neural information processing systems , 293 (2003).
[10] Z. F. Mainen and T. J. Sejnowski, Science 268, 1503
(1995).
[11] M. A. Petrovici, J. Bill, I. Bytschok, J. Schemmel, and
K. Meier, arXiv preprint arXiv:1311.3211 (2013).
[12] A. Destexhe, J. Comp. Neurosci. 27, 493 (2009).
[13] A. Destexhe, M. Rudolph, and D. Pare, Nature Reviews
Neuroscience 4, 739 (2003).
[14] M. J. Richardson and W. Gerstner, Neural computation
17, 923 (2005).
[15] L. M. Ricciardi and L. Sacerdote, Biol. Cybern. 35, 1
(1979).
[16] N. Brunel and S. Sergi, Journal of Theoretical Biology
195, 87 (1998).
[17] R. Moreno-Bote and N. Parga, Physical Review Letters
92, 028102 (2004).
[18] L. M. Ricciardi and S. Sato, Journal of Applied Probability 25, 43 (1988).
[19] D. Pecevski, L. Buesing, and W. Maass, PLoS computational biology 7, e1002294 (2011).
[20] R. Salakhutdinov and G. E. Hinton (2009) pp. 448–455.

[21] A.-R. Mohamed, G. E. Dahl, and G. Hinton, IEEE
Trans. on Audio, Speech, and Lang. Proc. 20, 14 (2012).
[22] M. Tsodyks and H. Markram, Proceedings of the National Academy of Sciences of the United States of America 94, 719 (1997).
[23] G. Hinton, P. Dayan, B. Frey, and R. Neal, Science 268,
1158 (1995).
[24] D. Sulzer and S. Rayport, Amino acids 19, 45 (2000).
[25] W. Gerstner and J. L. van Hemmen, Network: Computation in Neural Systems 3, 139 (1992).
[26] R. Brette and W. Gerstner, J. Neurophysiol. 94, 3637
(2005).
[27] S. Mitra, S. Fusi, and G. Indiveri, IEEE Transactions on
Biomedical Circuits and Systems 3, 32 (2009).
[28] M. A. Petrovici, B. Vogginger, P. Müller, O. Breitwieser,
et al., PLOS ONE 9, e108590 (2014).
[29] D. Probst, M. A. Petrovici, I. Bytschok, et al., Frontiers
in Computational Neuroscience 9, 13 (2015).
[30] E. Neftci, S. Das, B. Pedroni, et al., Frontiers in Neuromorphic Engineering (2014).
[31] P. Lánskỳ, Physical Review E 55, 2040 (1997).
[32] M. U. Thomas, Journal of Applied Probability , 600
(1975).
[33] A. N. Burkitt, Biological cybernetics 95, 97 (2006).
[34] M. L. Hines and N. T. Carnevale, The NEURON Book
(Cambridge University Press, Cambridge, UK, 2006).

16

FIG. 7. Sampling with LIF neurons over a broad range of relevant model parameters. All plots depict the Kullback-Leibler
divergence between the sampled distribution (with a certain set of parameters) and the target distribution, which is identical to
the one used in Fig. 3B in the main manuscript. (A) Sweep over neuron size and leakage, as well as background input parameters.
The axes represent multiplicative scaling values for four parameters: background (Poisson) synaptic weights wexc,inh and firing
rates νexc,inh on the abscisa, neuron capacitance Cm and leak conductance gl on the ordinate. The parameter values used
throughout the main manuscript (see also “Appendix VII: Simulation parameters”) therefore have the coordinates (1, 1). The
network simulation runtimes were chosen as Tsim = 106 ms. As expected, the large neuron / weak noise scenario (top left
square) does not permit accurate sampling, as the activation function is no longer logistic (see also Fig. 6F). In general, the
plot shows that good sampling quality can be achieved for any neuron capacitance and leak as long as the background noise is
strong enough (HCS). (B) Sweep over the ratio of the synaptic and refractory time constants. The best sampling performance
is, indeed, achieved for τsyn ≈ τref , but the network still produces good approximations of the target distribution when the
two time constants are not precisely identical. All DKL data points result from 20 simulations with runtimes of Tsim = 105 ms
each. The error bars represent the standard error of the mean. (C) Effect of an increased distance from reset to threshold
potential. An increase in ϑ − % causes a gradual decay of the sampling quality, since the membrane requires additional time
to reach a suprathreshold ueff when the refractory period is over. This can, in principle, be accomodated by defining a larger
time window τon during which the neuron is considered to encode the state z = 1. Nevertheless, in an HCS, the effective time
constant can be low enough to render the threshold-to-reset distance irrelevant.

[35] A. P. Davison, D. Brüderle, J. Eppler, J. Kremkow,
E. Muller, D. Pecevski, L. Perrinet, and P. Yger, Front.
Neuroinform. 2 (2008).
[36] R. Naud, N. Marcille, C. Clopath, and W. Gerstner,
Biological Cybernetics 99, 335 (2008).

[37] Y. LeCun and C. Cortes, (1998).

