Neural Dynamics as Sampling: A Model for Stochastic Computation in Recurrent Networks of Spiking Neurons

The organization of computations in networks of spiking neurons in the brain is still largely unknown, in particular in view of the inherently stochastic features of their firing activity and the experimentally observed trial-to-trial variability of neural systems in the brain. In principle there exists a powerful computational framework for stochastic computations, probabilistic inference by sampling, which can explain a large number of macroscopic experimental data in neuroscience and cognitive science. But it has turned out to be surprisingly difficult to create a link between these abstract models for stochastic computations and more detailed models of the dynamics of networks of spiking neurons. Here we create such a link and show that under some conditions the stochastic firing activity of networks of spiking neurons can be interpreted as probabilistic inference via Markov chain Monte Carlo (MCMC) sampling. Since common methods for MCMC sampling in distributed systems, such as Gibbs sampling, are inconsistent with the dynamics of spiking neurons, we introduce a different approach based on non-reversible Markov chains that is able to reflect inherent temporal processes of spiking neuronal activity through a suitable choice of random variables. We propose a neural network model and show by a rigorous theoretical analysis that its neural activity implements MCMC sampling of a given distribution, both for the case of discrete and continuous time. This provides a step towards closing the gap between abstract functional models of cortical computation and more detailed models of networks of spiking neurons.


Introduction
Attempts to understand the organization of computations in the brain from the perspective of traditional, mostly deterministic, models of computation, such as attractor neural networks or Turing machines, have run into problems: Experimental data suggests that neurons, synapses, and neural systems are inherently stochastic [1], especially in vivo, and therefore seem less suitable for implementing deterministic computations. This holds for ion channels of neurons [2], synaptic release [3], neural response to stimuli (trial-to-trial variability) [4,5], and perception [6]. In fact, several experimental studies arrive at the conclusion that external stimuli only modulate the highly stochastic spontaneous firing activity of cortical networks of neurons [7,8]. Furthermore, traditional models for neural computation have been challenged by the fact that typical sensory data from the environment is often noisy and ambiguous, hence requiring neural systems to take uncertainty about external inputs into account. Therefore many researchers have suggested that information processing in the brain carries out probabilistic, rather than logical, inference for making decisions and choosing actions [9][10][11][12][13][14][15][16][17][18][19][20][21][22]. Probabilistic inference has emerged in the 1960's [23], as a principled mathematical framework for reasoning in the face of uncertainty with regard to observations, knowledge, and causal relationships, which is characteristic for real-world inference tasks. This framework has become tremendously successful in real-world applications of artificial intelligence and machine learning. A typical computation that needs to be carried out for probabilistic inference on a high-dimensional joint distribution p(z 1 , . . . , z l ,z lz1 , . . . ,z K ) is the evaluation of the conditional distribution p(z 1 , . . . ,z l jz lz1 , . . . ,z K ) (or marginals thereof) over some variables of interest, say z 1 , . . . ,z l , given variables z lz1 , . . . ,z K . In the following, we will call the set of variables z lz1 , . . . ,z K , which we condition on, the observed variables and denote it by o.
Numerous studies in different areas of neuroscience and cognitive science have suggested that probabilistic inference could explain a variety of computational processes taking place in neural systems (see [10,11]). In models of perception the observed variables o are interpreted as the sensory input to the central nervous system (or its early representation by the firing response of neurons, e.g., in the LGN in the case of vision), and the variables z 1 , . . . ,z l model the interpretation of the sensory input, e.g., the texture and position of objects in the case of vision, which might be encoded in the response of neurons in various higher cortical areas [15]. Furthermore, in models for motor control the observed variables o often consist not only of sensory and proprioceptive inputs to the brain, but also of specific goals and constraints for a planned movement [24][25][26], whereas inference is carried out over the variables z 1 , . . . ,z l representing a motor plan or motor commands to muscles. Recent publications show that human reasoning and learning can also be cast into the form of probabilistic inference problems [27][28][29]. In these models learning of concepts, ranging from concrete to more abstract ones, is interpreted as inference in lower and successively higher levels of hierarchical probabilistic models, giving a consistent description of inductive learning within and across domains of knowledge.
In spite of this active research on the functional level of neural processing, it turned out to be surprisingly hard to relate the computational machinery required for probabilistic inference to experimental data on neurons, synapses, and neural systems. There are mainly two different approaches for implementing the computational machinery for probabilistic inference in ''neural hardware''. The first class of approaches builds on deterministic methods for evaluating exactly or approximately the desired conditional and/or marginal distributions, whereas the second class relies on sampling from the probability distributions in question. Multiple models in the class of deterministic approaches implement algorithms from machine learning called message passing or belief propagation [30][31][32][33]. By clever reordering of sum and product operators occurring in the evaluation of the desired probabilities, the total number of computation steps are drastically reduced. The results of subcomputations are propagated as "messages" or "beliefs" that are sent to other parts of the computational network. Other deterministic approaches for representing distributions and performing inference are probabilistic population code (PPC) models [34]. Although deterministic approaches provide a theoretically sound hypothesis about how complex computations can possibly be embedded in neural networks and explain aspects of experimental data, it seems difficult (though not impossible) to conciliate them with other aspects of experimental evidence, such as stochasticity of spiking neurons, spontaneous firing, trial-to-trial variability, and perceptual multistability.
Therefore other researchers (e.g., [16][17][18]35]) have proposed to model computations in neural systems as probabilistic inference based on a different class of algorithms, which requires stochastic, rather than deterministic, computational units. This approach, commonly referred to as sampling, focuses on drawing samples, i.e., concrete values for the random variables that are distributed according to the desired probability distribution. Sampling can naturally capture the effect of apparent stochasticity in neural responses and seems to be furthermore consistent with multiple experimental effects reported in cognitive science literature [17,18]. On the conceptual side, it has proved to be difficult to implement learning in message passing and PPC network models. In contrast, following the lines of [36], the sampling approach might be well suited to incorporate learning.
Previous network models that implement sampling in neural networks are mostly based on a special sampling algorithm called Gibbs (or general Metropolis-Hastings) sampling [9,17,18,37]. The dynamics that arise from this approach, the so-called Glauber dynamics, however are only superficially similar to spiking neural dynamics observed in experiments, rendering these models rather abstract. Building on and extending previous models, we propose here a family of network models, that can be shown to exactly sample from any arbitrary member of a well-defined class of probability distributions via their inherent network dynamics. These dynamics incorporate refractory effects and finite durations of postsynaptic potentials (PSPs), and are therefore more biologically realistic than existing approaches. Formally speaking, our model implements Markov chain Monte Carlo (MCMC) sampling in a spiking neural network. In contrast to prior approaches however, our model incorporates irreversible dynamics (i.e., no detailed balance) allowing for finite time PSPs and refractory mechanisms. Furthermore, we also present a continuous time version of our network model. The resulting stochastic dynamical system can be shown to sample from the correct distribution. In general, continuous time models arguably provide a higher amount of biological realism compared to discrete time models.
The paper is structured in the following way. First we provide a brief introduction to MCMC sampling. We then define the neural network model whose neural activity samples from a given class of probability distributions. The model will be first presented in discrete time together with some illustrative simulations. An extension of the model to networks of more detailed spiking neuron models which feature a relative refractory mechanism is presented. Furthermore, it is shown how the neural network model can also be formulated in continuous time. Finally, as a concrete simulation example we present a simple network model for perceptual multistability.

Recapitulation of MCMC sampling
In machine learning, sampling is often considered the ''gold standard'' of inference methods, since, assuming that we can sample from the distribution in question, and assuming enough computational resources, any inference task can be carried out with arbitrary precision (in contrast to some deterministic approximate inference methods such as variational inference). However sampling from an arbitrary distribution can be a difficult problem in itself, as, e.g., many distributions can only be evaluated modulo a global constant (the partition function). In order to circumvent these problems, elaborate MCMC sampling techniques have been developed in machine learning and statistics [38]. MCMC algorithms are based on the following idea: instead of producing an ad-hoc sample, a process that is heuristically comparable to a global search over the whole state space of the

Author Summary
It is well-known that neurons communicate with short electric pulses, called action potentials or spikes. But how can spiking networks implement complex computations? Attempts to relate spiking network activity to results of deterministic computation steps, like the output bits of a processor in a digital computer, are conflicting with findings from cognitive science and neuroscience, the latter indicating the neural spike output in identical experiments changes from trial to trial, i.e., neurons are ''unreliable''. Therefore, it has been recently proposed that neural activity should rather be regarded as samples from an underlying probability distribution over many variables which, e.g., represent a model of the external world incorporating prior knowledge, memories as well as sensory input. This hypothesis assumes that networks of stochastically spiking neurons are able to emulate powerful algorithms for reasoning in the face of uncertainty, i.e., to carry out probabilistic inference. In this work we propose a detailed neural network model that indeed fulfills these computational requirements and we relate the spiking dynamics of the network to concrete probabilistic computations. Our model suggests that neural systems are suitable to carry out probabilistic inference by using stochastic, rather than deterministic, computing elements. random variables, MCMC methods produce a new sample via a ''local search'' around a point in the state space that is already (approximately) a sample from the distribution.
More formally, a Markov chain M (in discrete time) is defined by a set S of states (we consider for discrete time only the case where S has a finite size, denoted by jSj) together with a transition operator T. The operator T is a conditional probability distribution T(sjs') over the next state s given a preceding state s'. The Markov chain M is started in some initial state s(0), and moves through a trajectory of states s(t) via iterated application of the stochastic transition operator T. More precisely, if s(t{1) is the state at time t{1, then the next state s(t) is drawn from the conditional probability distribution T(sjs(t{1)). An important theorem from probability theory (see, e.g., p. 232 in [39]) states that if M is irreducible (i.e., any state in S can be reached from any other state in S in finitely many steps with probability w0) and aperiodic (i.e., its state transitions cannot be trapped in deterministic cycles), then the probability p(s(t)~sjs(0)) converges for t?? to a probability p(s) that does not depend on the initial state s(0). This state distribution p is called the invariant distribution of M. The irreducibility of M implies that it is the only distribution over the states S that is invariant under its transition operator T, i.e.
Thus, in order to carry out probabilistic inference for a given distribution p, it suffices to construct an irreducible and aperiodic Markov chain M that leaves p invariant, i.e., satisfies equation (1).
Then one can answer numerous probabilistic inference questions regarding p without any numerical computations of probabilities. Rather, one plugs in the observed values for some of the random variables (RVs) and simply collects samples from the conditional distribution over the other RVs of interest when the Markov chain approaches its invariant distribution.
A convenient and popular method for the construction of an operator T for a given distribution p is looking for operators T that satisfy the following detailed balance condition, for all s, s'[S. A Markov chain that satisfies (2) is said to be reversible. In particular, the Gibbs and Metropolis-Hastings algorithms employ reversible Markov chains. A very useful property of (2) is that it implies the invariance property (1), and this is in fact the standard method for proving (1). However, as our approach makes use of irreversible Markov chains as explained below, we will have to prove (1) directly.

Neural sampling
Let p(z 1 , . . . ,z K ) be some arbitrary joint distribution over K binary variables z 1 , . . . ,z K that only takes on values w0. We will show that under a certain computability assumption on p a network N consisting of K spiking neurons n 1 , . . . ,n K can sample from p using its inherent stochastic dynamics. More precisely, we show that the stochastic firing activity of N can be viewed as a non-reversible Markov chain that samples from the given probability distribution p. If a subset o of the variables are observed, modelled as the corresponding neurons being ''clamped'' to the observed values, the remaining network samples from the conditional distribution of the remaining variables given the observables. Hence, this approach offers a quite natural implementation of probabilistic inference. It is similar to sampling approaches which have already been applied extensively, e.g., in Boltzmann machines, however our model is more biologically realistic as it incorporates aspects of the inherent temporal dynamics and spike-based communication of a network of spiking neurons. We call this approach neural sampling in the remainder of the paper.
In order to enable a network N of spiking neurons to sample from a distribution p(z 1 , . . . ,z K ) of binary variables z k , one needs to specify how an assignment (z 1 , . . . ,z K )[f0,1g K of values to these binary variables can be represented by the spiking activity of the network N and vice versa. A spike, or action potential, of a biological neuron n k has a short duration of roughly 1 ms. But the effect of such spike, both on the neuron n k itself (in the form of refractory processes) and on the membrane potential of other neurons (in the form of postsynaptic potentials) lasts substantially longer, on the order of 5 ms to 100 ms. In order to capture this temporally extended effect of each spike, we fix some parameter t that models the average duration of these temporally extended processes caused by a spike. We say that a binary vector (z 1 , . . . ,z K ) is represented by the firing activity of the network N at time t for k~1, . . . ,K iff: z k (t)~1un k has fired within the time interval (t{t, t: ð3Þ In other words, any spike of neuron n k sets the value of the associated binary variable z k to 1 for a duration of length t. An obvious consequence of this definition is that the binary vector (z 1 , . . . ,z K ) that is defined by the activity of N at time t does not fully capture the internal state of this stochastic system. Rather, one needs to take into account additional non-binary variables (f 1 , . . . ,f K ), where the value of f k at time t specifies when within the time interval (t{t, t the neuron n k has fired (if it has fired within this time interval, thereby causing z k~1 at time t). The neural sampling process has the Markov property only with regard to these more informative auxiliary variables f 1 , . . . ,f K . Therefore our analysis of neural sampling will focus on the temporal evolution of these auxiliary variables. We adopt the convention that each spike of neuron n k sets the value of f k to its maximal value t, from which it linearly decays back to 0 during the subsequent time interval of length t.
For the construction of the sampling network N , we assume that the membrane potential u k (t) of neuron n k at time t equals the log-odds of the corresponding variable z k to be active, and refer to this property as neural computability condition: where we write z k for z k (t) and z \k for the current values z i (t) of all other variables z i with i=k. Under the assumption we make in equation (4), i.e., that the neural membrane potential reflects the log-odds of the corresponding variable z k , it is required that each single neuron in the network can actually compute the right-hand side of equation (4), i.e., that it fulfills the neural computability condition.
A concrete class of probability distributions, that we will use as an example in the remainder, are Boltzmann distributions: with arbitrary real valued parameters b i , W ij which satisfy W ij~Wji and W ii~0 (the constant Z ensures the normalization of p(z)). For the Boltzmann distribution, condition (4) is satisfied by neurons n k with the standard membrane potential where b k is the bias of neuron n k (which regulates its excitability), W ki is the strength of the synaptic connection from neuron n i to n k , and W ki z i (t) approximates the time course of the postsynaptic potential in neuron n k caused by a firing of neuron n i with a constant signal of duration t (i.e., a square pulse). As we will describe below, spikes of neuron n k are evoked stochastically depending on the current membrane potential u k and the auxiliary variable f k . The neural computability condition (4) links classes of probability distributions to neuron and synapse models in a network of spiking neurons. As shown above, Boltzmann distributions satisfy the condition if one considers point neuron models which compute a linear weighted sum of the presynaptic inputs. The class of distributions can be extended to include more complex distributions using a method proposed in [40] which is based on the following idea. Neuron n k representing the variable z k is not directly influenced by the activities z \k of the presynaptic neurons, but via intermediate nonlinear preprocessing elements. This preprocessing might be implemented by dendrites or other (inter-) neurons and is assumed to compute nonlinear combinations of the presynaptic activities z \k (similar to a kernel). This allows the membrane potential u k , and therefore the log-odds ratio on the right-hand side of (4), to represent a more complex function of the activities z \k , giving rise to more complex joint distributions p(z). The concrete implementation of non-trivial directed and undirected graphical models with the help of preprocessing elements in the neural sampling framework is subject of current research. For the examples given in this study, we focus on the standard form of the membrane potential (6) of point neurons. As shown below, these spiking network models can emulate any Boltzmann machine (BM) [36].
A substantial amount of preceding studies has demonstrated that BMs are very powerful, and that the application of suitable learning algorithms for setting the weights W ij makes it possible to learn and represent complex sensory processing tasks by such distributions [37,41]. In applications in statistics and machine learning using such Boltzmann distributions, sampling is typically implemented by Gibbs sampling or more general reversible MCMC methods. However, it is difficult to model some neural processes, such as an absolute refractory period or a postsynaptic potential (PSP) of fixed duration, using a reversible Markov chain, but they are more conveniently modelled using an irreversible one. As we wish to keep the computational power of BMs and at the same time to augment the sampling procedure with aspects of neural dynamics (such as PSPs with fixed durations, refractory mechanisms) to increase biological realism, we focus in the following on irreversible MCMC methods (keeping in mind that this might not be the only possible way to achieve these goals).

Neural sampling in discrete time
Here we describe neural dynamics in discrete time with an absolute refractory period t. We interpret one step of the Markov chain as a time step dt in biological real time. The dynamics of the variable f k , that describes the time course of the effect of a spike of neuron n k , are defined in the following way. f k is set to the value t when neuron n k fires, and decays by 1 at each subsequent discrete time step. The parameter t is chosen to be some integer, so that f k decays back to 0 in exactly t time steps. The neuron can only spike (with a probability that is a function of its current membrane potential u k ) if its variable f k ƒ1. If however, f k w1, the neuron is considered refractory and it cannot spike, but its f k is reduced by 1 per time step. To show that these simple dynamics do indeed sample from the given distribution p(z), we proceed in the following way. We define a joint distribution p(f,z) which has the desired marginal distribution P f p(f,z)~p(z). Further we formalize the dynamics informally described above as a transition operator T operating on the state vector (f,z). Finally, in the Methods section, we show that p(f,z) is the unique invariant distribution of this operator T, i.e., that the dynamics described by T produce samples z from the desired distribution p(z). We refer to sampling through networks with this stochastic spiking mechanism as neural sampling with absolute refractory period due to the persistent refractory process.
Given the distribution p(z) that we want to sample from, we define the following joint distribution p(f,z) over the neural variables: This definition of p(f k jz k ) simply expresses that if z k~1 , then the auxiliary variable f k can assume any value in f1,2, . . . ,tg with equal probability. On the other hand f k necessarily assumes the value 0 if z k~0 (i.e., when the neuron is in its resting state).
The state transition operator T can be defined in a transparent manner as a composition of K transition operators, TT 1 0 . . . 0T K , where T k only updates the variables f k and z k of neuron n k , i.e., the neurons are updated sequentially in the same order (this severe restriction will become obsolete in the case of continuous time discussed below). We define the composition as (T k 0T l )( : )~(T k (T l ( : )), i.e., T l is applied prior to T k . The new values of f k and z k only depend on the previous value f' k and on the current membrane potential u k (z \k ). The interesting dynamics take place in the variable f k . They are illustrated in Figure 1 where the arrows represent transition probabilities greater than 0.
If the neuron n k is not refractory, i.e., f' k ƒ1, it can spike (i.e., a transition from f' k ƒ1 to f k~t ) with probability where s(x)~(1ze {x ) {1 is the standard sigmoidal activation function and the log denotes the natural logarithm. The term u k is the current membrane potential, which depends on the current values of the variables z i for i=k. The term log t in (8) reflects the granularity of a chosen discrete time scale. If it is very fine (say one step equals one microsecond), then t is large, and the firing probability at each specific discrete time step is therefore reduced. If the neuron in a state with f' k ƒ1 does not spike, f k relaxes into the resting state f k~0 corresponding to a non-refractory neuron. If the neuron is in a refractory state, i.e., f' k w1, its new variable f k assumes deterministically the next lower value f k~f ' k {1, reflecting the inherent temporal process: After the transition of the auxiliary variable f k , the binary variable z k is deterministically set to a consistent state, i.e., z k~1 if f k §1 and z k~0 if f k~0 .
It can be shown that each of these stochastic state transition operators T k leaves the given distribution p invariant, i.e., satisfies equation (1). This implies that any composition or mixture of these operators T k also leaves p invariant, see, e.g., [38]. In particular, the composition T~T 1 0 . . . 0T K of these operators T k leaves p invariant, which has a quite natural interpretation as firing dynamics of the spiking neural network N : At each discrete time step the variables f k ,z k are updated for all neurons n k , where the update of f k ,z k takes preceding updates for f i ,z i with iwk into account. Alternatively, one could also choose at each discrete time step a different order for updates according to [38]. The assumption of a well-regulated updating policy will be overcome in the continuous-time limit, i.e., in case where the neural dynamics are described as a Markov jump process. In the methods section we prove the following central theorem: is the unique invariant distribution of operator T, i.e., T is aperiodic and irreducible and satisfies The proof of this Theorem is provided by Lemmata 1 -3 in the Methods section. The statement that T (which is composed of the operators T k ) is irreducible and aperiodic ensures that p is the unique invariant distribution of the Markov chain defined by T, i.e., that irrespective of the initial network state the successive application of T explores the whole state space in a non-periodic manner.
This theorem guarantees that after a sufficient ''burn-in'' time (more precisely in the limit of an infinite ''burn-in'' time), the dynamics of the network, which are given by the transition operator T, produce samples from the distribution p(f,z). As by construction P f p(f,z)~p(z), the Markov chain provides samples from the given distribution p(z). Furthermore, the network N can carry out probabilistic inference for this distribution. For example, N can be used to sample from the posterior distribution p(z 1 . . . ,z l jz lz1 , . . . ,z K ) over z 1 . . . ,z l given z lz1 , . . . ,z K . One just needs to clamp those neurons n lz1 , . . . ,n K to the corresponding observed values. This could be implemented by injecting a strong positive (negative) current into the units with z j~1 (z j~0 ). Then, as soon as the stochastic dynamics of N has converged to its invariant distribution, the averaged firing rate of neuron n 1 is proportional to the following desired marginal probability p(z 1~1 jz lz1 , . . . ,z K )~X z 2 ,...,z l p(z 1~1 ,z 2 , . . . ,z l jz lz1 , . . . ,z K ): In a biological neural system this result of probabilistic inference could for example be read out by an integrator neuron that counts spikes from this neuron n 1 within a behaviorally relevant time window of a few hundred milliseconds, similarly as the experimentally reported integrator neurons in area LIP of monkey cortex [20,21]. Another readout neuron that receives spike input from n k could at the same time estimate p(z k~1 jz lz1 , . . . ,z K ) for another RV z k . But valuable information for probabilistic inference is not only provided by firing rates or spike counts, but also by spike correlations of the neurons n 1 , . . . ,n l in N . For example, the probability p(z 1~1 ,z 2~1 jz lz1 , . . . ,z K ) can be estimated by a readout neuron that responds to superpositions of EPSPs caused by near-coincident firing of neurons n 1 and n 2 within a time interval of length t. Thus, a large number of different probabilistic inferences can be carried out efficiently in parallel by readout neurons that receive spike input from different subsets of neurons in the network N .
Variation of the discrete time model with a relative refractory mechanism. For the previously described simple neuron model, the refractory process was assumed to last for t time steps, exactly as long as the postsynaptic potentials caused by each spike. In this section we relax this assumption by introducing a more complex and biologically more realistic neuron model, where the duration of the refractory process is decoupled from the duration t of a postsynaptic potential. Thus, this model can for example also fire bursts of spikes with an interspike interval vt. The introduction of this more complex neuron model comes at the price that one can no longer prove that a network of such neurons samples from the desired distribution p. Nevertheless, if the sigmoidal activation function s is replaced by a different activation function f , one can still prove that the sampling is ''locally correct'', as specified in equation (12) below. Furthermore, our computer simulations suggest that also globally the error introduced by the more complex neuron model is not functionally significant, i.e. that statistical dependencies between the RVs z are still faithfully captured.
The neuron model with a relative refractory period is defined in the following way. Consider some arbitrary refractory function g : ½0, . . . ,t?R with g(t)~0, g(0)~1, and g(l) §0 for l~1, . . . ,t{1. The idea is that g(f k ) models the readiness of the neuron to fire in its state f k . This readiness has value 0 when the neuron has fired at the preceding time step (i.e., f k~t ), and assumes the resting state 1 when f k has dropped to 0. In between, the readiness may take on any non-negative value according to the function g(f k ). The function g does not need to be monotonic, allowing for example that it increases to high values in between, yielding a preferred interspike interval of a oscillatory neuron. The firing probability of neuron n k in state f k is given by is an appropriate function of the membrane potential as described below. Thus this function g is closely related to the function g (called afterpotential) in the spike response model [5] as well as to the self-excitation kernel in Generalized Linear Models [42]. In general, different neurons in the network may have different refractory profiles, which can be modeled by a different refractory function for each neuron n k . However for the sake of notational simplicity we assume a single refractory function in the following.
In the presence of this refractory function g one needs to replace the sigmoidal activation function s(u k { log t) by a suitable function f (u k ) that satisfies the condition for all real numbers u. This equation can be derived (see Methods section Lemma 0) if one requires each neuron n k to represent the correct distribution p(z k jz \k ) over z k conditioned the variables z \k . One can show that, for any g as above, there always exists a continuous, monotonic function f which satisfies this equation (see Lemma 0 in Methods). Unfortunately (11) cannot be solved analytically for f in general. Hence, for simulations we approximate the function f for a given g by numerically solving (11) on a grid and interpolating between the grid points with a constant function. Examples for several functions g and the associated f are shown in Figure 2B and Figure 2C respectively. Furthermore, spike trains emitted by single neurons with these refractory functions g and the corresponding functions f are shown in Figure 2D for the case of piecewise constant membrane potentials. This figure indicates, that functions g that define a shorter refractory effect lead to higher firing rates and more irregular firing. It is worth noticing that the standard activation function s(u k { log t) is the solution of equation (11) for the absolute refractory function, i.e., for g(0)~g(1)~1 and g(l)~0 for 1vlƒt. The transition operator T k is defined for this model in a very similar way as before. However, for 1vf' k ƒt, when the variable f' k was deterministically reduced by 1 in the simpler model (yielding f k~f ' k {1), this reduction occurs now only with probability 1{g(f' k ) : f (u k ). With probability g(f' k ) : f (u k ) the operator T k sets f k~t , modeling the firing of another spike of neuron n k at this time point. The neural computability condition (4) remains unchanged, e.g., u k~bk z P K i~1 W ki z i for a Boltzmann distribution. A schema of the stochastic dynamics of this local state transition operator T k (f k jf' k ,z' \k ) is shown in Figure 2A.
This transition operator T k has the following properties. In Lemma 0 in Methods it is proven that the unique invariant distribution of T k , denoted as q Ã k (f k ,z k jf \k ,z \k ), gives rise to the correct marginal distribution over z k , i.e.
This means that a neuron whose dynamics is described by T k samples from the correct distribution p(z k jz \k ) if it receives a static input from the other neurons in the network, i.e., as long as its membrane potential u k is constant. Hence the ''local'' computa- . The x-axis is equivalent to the number of time steps since last spike (running from 0 to t from left to right). (C) Associated activation functions f according to (11). (D) Spike trains produced by the resulting three different neuron models with (hypothetical) membrane potentials that jump at time ½0:25s from a constant low value to a constant high value. Black horizontal bars indicate spikes, and the active states z k~1 are indicated by gray shaded areas of duration t : dt~20ms after each spike. It can be seen from this example that different refractory mechanisms give rise to different spiking dynamics. doi:10.1371/journal.pcbi.1002211.g002 tion performed by such neuron can be considered as correct. If however, several neurons in the network change their states in a short interval of time, the joint distribution over z is in general not the desired one, i.e., In the Methods section, we present simulation results that indicate that the error of the approximation to the desired Boltzmann distributions introduced by neural sampling with relative refractory mechanism is rather minute. It is shown that the neural sampling approximation error is orders of magnitudes below the one introduced by a fully factorized distribution (which amounts to assuming correct marginal distributions p(z k ) and independent neurons).
To illustrate the sampling process with the relative refractory mechanism, we examine a network of K~40 neurons. We aim to sample from a Boltzmann distribution (5) with parameters W ij , b i being randomly drawn from normal distributions. For the neuron model, we use the relative refractory mechanism shown in the mid row of Figure 2B. A detailed description of the simulation and the parameters used is given in the Methods section. A spike pattern of the resulting sampling network is shown in Figure 3A. The network features a sparse, irregular spike response with average firing rate of 13:9 Hz. For one neuron n 26 , indicated with orange spikes, the internal dynamics are shown in Figure 3B. After each action potential the neuron's refractory function g(f 26 ) drops to zero and reduces the probability of spiking again in a short time interval. The influence of the remaining network z \26 is transmitted to neuron n 26 via PSPs of duration t : dt~20 ms and sums up to the fluctuating membrane potential u 26 . As reflected in the highly variable membrane potential even this small network exhibits rich interactions. To represent the correct distribution p(z 26 jz \26 ) over z 26 conditioned on z \26 , the neuron n 26 continuously adapts its instantaneous firing rate. To quantify the precision with which the spiking network draws samples from the target distribution (5), Figure 3C shows the joint distribution of 5 neurons. For comparison we accompany the distribution of sampled network states with the result obtained from the standard Gibbs sampling algorithm (considered as the ground truth). Since the number of possible states z grows exponentially in the number of neurons, we restrict ourselves for visualization purposes to the distribution p(z 24 , . . . ,z 28 ) of the gray shaded units and marginalize over the remaining network. The probabilities are estimated from 10 7 samples, i.e., from 10 7 successive states z of the Markov chain. Stochastic deviations of the estimated probabilities due to the finite number of samples are quite small (typical errors and are comparable to systematic deviations due to the only locally correct computation of neurons with relative refractory mechanism. In the Methods section, we present further simulation results showing that the proposed networks consisting of neurons with relative refractory mechanism approximate the desired target distributions faithfully over a large range of distribution parameters. In order to illustrate that the proposed sampling networks feature biologically quite realistic spiking dynamics, we present in the Methods section several neural firing statistics (e.g., the interspike interval histogram) of the network model. In general, the statistics computed from the model match experimentally observed statistics well. The proposed network models are based on the assumption of rectangular-shaped, renewal PSPs. More precisely, we define renewal (or non-additive) PSPs in the following way. Renewal PSPs evoked by a single synapse do not add up but are merely prolonged in their duration (according to equation (6)); renewal PSPs elicited at different synapses nevertheless add up in the normal way. In Methods we investigate the impact of replacing the theoretically ideal rectangular-shaped, renewal PSPs with biologically more realistic alpha-shaped, additive PSPs. Simulation results suggest that the network model with alpha-shaped PSPs does not capture the target distribution as accurately as with the theoretically ideal PSP shapes, statistical dependencies between the RVs z are however still approximated reasonably well. Neural sampling in continuous time The neural sampling model proposed above was formulated in discrete time of step size dt, inspired by the discrete time nature of MCMC techniques in statistics and machine learning as well as to make simulations possible on digital computers. However, models in continuous time (e.g., ordinary differential equations) are arguably more natural and ''realistic'' descriptions of temporally varying biological processes. This gives rise to the question whether one can find a sensible limit of the discrete time model in the limit dt?0, yielding a sampling network model in continuous time. Another motivation for considering continuous time models for neural sampling is the fact that many mathematical models for recurrent networks are formulated in continuous time [5], and a comparison to these existing models would be facilitated. Here we propose a stochastically spiking neural network model in continuous time, whose states still represent correct samples from the desired probability distribution p(z) at any time t. These types of models are usually referred to as Markov jump processes. It can be shown that discretizing this continuous time model yields the discrete time model defined earlier, which thus can be regarded as a version suitable for simulations on a digital computer.
We define the continuous time model in the following way. Let t l k , for l~0,1, . . ., denote the firing times of neuron n k . The refractory process of this neuron, in analogy to Figure 1 and equation (8)-(9) for the case of discrete time, is described by the following differential equation for the auxiliary variable f k , which may now assume any nonnegative real number 0ƒf k ƒ1 : Here d(t{t l k ) denotes Dirac's Delta centered at the spike time t l k . This differential equation describes the following simple dynamics. The auxiliary variable f k (t) decays linearly with time constant t when the neuron is refractory, i.e., f k (t)w0. Once f k (t) arrives at its resting state 0 it remains there, corresponding to the neuron being ready to spike again (more precisely, in order to avoid point measures we set it to a random value in ½{2E,{E, see Methods). In the resting state, the neuron has the probability density 1 t exp (u k (t)) to fire at every time t. If it fires at t l k , this results in setting f k (t l k )~1, which is formalized in equation (12) by the sum of Dirac Delta's P l d(t{t l k ). Here the current membrane potential u k (t) at time t is defined as in the discrete time case, e.g., by u k~bk z P K i~1 W ki z i (t) for the case of a Boltzmann distribution (5). The binary variable z k (t) is defined to be 1 if f k (t)w0 and 0 if the neuron is in the resting state f k (t)~0. Biologically, the term W ki z i (t) can again be interpreted as the value at time t of a rectangular-shaped PSP (with a duration of t) that neuron n i evokes in neuron n k . As the spikes are discrete events in continuous time, the probability of two or more neurons spiking at the same time is zero. This allows for updating all neurons in parallel using a differential equation.
In analogy to the discrete time case, the neural network in continuous time can be shown to sample from the desired distribution p(z), i.e., p(z) is an invariant distribution of the network dynamics defined above. However, to establish this fact, one has to rely on a different mathematical framework. The probability distribution p t (f) of the auxiliary variables f 1 (t), . . . ,f K (t) as a function of time t, which describes the evolution of the network, obeys a partial differential equation, the so-called Differential-Chapman-Kolmogorov equation (see [43]): where the operator T, which captures the dynamics of the network, is implicitly defined by the differential equations (12) and the spiking probabilities. This operator T is the continuous time equivalent to the transition operator T in the discrete time case. The operator T consists here of two components. The drift term captures the deterministic decay process of f k (t), stemming from the term {1=t in equation (12). The jump term describes the noncontinuous aspects of the path at the time t l k when the neuron fires.
In the Methods section we prove that the resulting time invariant distribution, i.e., the distribution that solves L t p t (f)~0, now denoted p(f) as it is not a function of time, gives rise to the desired marginal distribution p(z) over z: where and 0 otherwise. For an explicit definition of T, a proof of the above statement, and some additional comments see the Methods section.
The neural samplers in discrete and continuous time are closely related. The model in discrete time provides an increasingly more precise description of the inherent spike dynamics when the duration dt of the discrete time step is reduced, causing an increase of t (such that t : dt is constant) and therefore a reduced firing probability of each neuron at any discrete time step (see the term log t in equation (8)). In the limit of dt approaching 0, the probability that two or more neurons will fire at the same time approaches 0, and the discrete time sampler becomes equal to the continuous time system defined above, which updates all units in parallel.
It is also possible to formulate a continuous time version of the neural sampler based on neuron models with relative refractory mechanisms. In the Methods section the resulting continuous time neuron model with a relative refractory mechanism is defined. Theoretical results similar to the discrete time case can be derived for this sampler (see Lemmata 9 and 10 in Methods): It is shown that each neuron ''locally'' performs the correct computation under the assumption of static input from the remaining neurons. However one can no longer prove in general that the global network samples from the target distribution p.
Demonstration of probabilistic inference with recurrent networks of spiking neurons in an application to perceptual multistability In the following we present a network model for perceptual multistability based on the neural sampling framework introduced above. This simulation study is aimed at showing that the proposed network can indeed sample from a desired distribution and also perform inference, i.e., sample from the correct corresponding posterior distribution. It is not meant to be a highly realistic or exhaustive model of perceptual multistability nor of biologically plausible learning mechanisms. Such models would naturally require considerably more modelling work.
Perceptual multistability evoked by ambiguous sensory input, such as a 2D drawing (e.g., Necker cube) that allows for different consistent 3D interpretations, has become a frequently studied perceptual phenomenon. The most important finding is that the perceptual system of humans and nonhuman primates does not produce a superposition of different possible percepts of an ambiguous stimulus, but rather switches between different selfconsistent global percepts in a spontaneous manner. Binocular rivalry, where different images are presented to the left and right eye, has become a standard experimental paradigm for studying this effect [44][45][46][47]. A typical pair of stimuli are the two images shown in Figure 4A. Here the percepts of humans and nonhuman primates switch (seemingly stochastically) between the two presented orientations. [16][17][18] propose that several aspects of experimental data on perceptual multistability can be explained if one assumes that percepts correspond to samples from the conditional distribution over interpretations (e.g., different 3D shapes) given the visual input (e.g., the 2D drawing). Furthermore, the experimentally observed fact that percepts tend to be stable on the time scale of seconds suggests that perception can be interpreted as probabilistic inference that is carried out by MCMC sampling which produces successively correlated samples. In [18] it is shown that this MCMC interpretation is also able to qualitatively reproduce the experimentally observed distribution of dominance durations, i.e., the distribution of time intervals between perceptual switches. However, in lack of an adequate model for sampling by a recurrent network of spiking neurons, theses studies could describe this approach only on a rather abstract level, and pointed out the open problem to relate this algorithmic approach to neural processes. We have demonstrated in a computer simulation that the previously described model for neural sampling could in principle fill this gap, providing a modelling framework that is on the one hand consistent with the dynamics of networks of spiking neurons, and which can on the other hand also be clearly understood from the perspective of probabilistic inference through MCMC sampling.
In the following we model some essential aspects of an experimental setup for binocular rivalry with grating stimuli (see Figure 4A) in a recurrent network of spiking neurons with the previously described relative refractory mechanism. We assigned to each of the 217 neurons in the network N a tuning curve V k (Q), centered around its preferred orientation Q Q k as shown in Figure 4B. The preferred orientations Q Q k of the neurons were chosen to cover the entire interval ½0,p) of possible orientations and were randomly assigned to the neurons. The neurons were arranged on a hexagonal grid as depicted in Figure 4F. Any two neurons with distance ƒ8 were synaptically connected (neighboring units had distance 1). We assume that these neurons represent neurons in the visual system that have roughly the same or neighboring receptive field, and that each neuron receives visual input from either the left or the right eye. The network connections were chosen such that neurons that have similar (very different) preferred orientations are connected with positive (negative) weights (for details see Methods section).
We examined the resulting distribution p(z) over the 217 dimensional network states. To provide an intuitive visualization of these high dimensional network states z, we resort to a 2dimensional projection, the population vector of a state z (see Methods for details of the applied population vector decoding scheme). Only the endpoints of the population vectors are drawn (as colored points) in Figure 4D,E. The orientation of the population vector is assumed to correspond to the dominant orientation of the percept, and its distance from the origin encodes the strength of this percept. We also, somewhat informally, call the strength of a percept its coherence and a network state which represents a coherent percept a coherent network state. A coherent network state hence results in a population vector of large magnitude. Each direction of a population vector is color coded in Figure 4D,E, using the color code for directions shown on the right hand side of Figure 4F. In Figure 4D the distribution p(z) of the network is illustrated by sampling of the network for ½ 20s, with samples z taken every millisecond. Each dot equals a sampled network state z. In a biological interpretation the spike response of the freely evolving network reflects spontaneous activity, since no observations, i.e., no external input, was added to the system. Figure 4D shows that the spontaneous activity of this simple network of spiking neurons moves preferably through coherent network states for all possible orientations due to the chosen recurrent network connections (being positive for neurons with similar preferred orientation and negative otherwise). This can directly be seen from the rare occurrence of population vectors with small magnitude (vectors close to the ''center'') in Figure 4D.
To study percepts elicited by ambiguous stimuli, where inputs like in Figure 4A are shown simultaneously to the left and right eye during a binocular rivalry experiment, we provided ambiguous input to the network. Two cells with preferred orientation Q Q k &45 0 and two cells with Q Q k &135 0 were clamped to 1. Additionally four neurons with Q Q k &0 0 resp. 90 0 were muted by clamping to 0. This ambiguous input is incompatible with a coherent percept, as it corresponds to two orthogonal orientations presented at the same time. The resulting distribution over the state of the 209 remaining neurons is shown for a time span of 20 s of simulated biological time (with samples taken every millisecond) in Figure 4E. One clearly sees that the network spends most of the time in network states that correspond to one of the two simultaneously presented input orientations (45 o and 135 o ), and virtually no time on orientations in between. This implements a sampling process from a bimodal conditional distribution. The black line marks a 500 ms trace of network states z around a perceptual switch: The network remained in one mode of high probability -corresponding to one percept -for some period of time, and then quickly traversed the state space to another mode -corresponding to a different percept.
Three of the states z around this perceptual switch (z(t 1 ), z(t 2 ) and z(t 3 ) in Figure 4E) are explicitly shown in Figure 4F. Neurons n k that fired during the preceding interval of 20 ms (marked in gray in Figure 4G) are drawn in the respective color of their preferred orientation. Inactive neurons are drawn in white, and clamped neurons are marked by a black dot (.). Figure 4G shows the action potentials of the 209 non-clamped neurons during the same 500 ms trace around the perceptual switch. One sees that the sampling process is expressed in this neural network model by a sparse, asynchronous and irregular spike response. It is worth mentioning that the average firing rate when sampling from the posterior distribution is only slightly higher than the average firing rate of spontaneous activity (16:1 Hz and 15:4 Hz respectively), which is reminiscent of related experimental data [7]. Thus on the basis of the overall network activity it is indistinguishable whether the network carries out an inference task or freely samples from its prior distribution. It is furthermore notable, that a focus of the network activity on the two orientations that are given by the external input can be achieved in this model, in spite of the fact that only two of the 217 neurons were clamped for each of them. This numerical relationship is reminiscent of standard data on the weak input from LGN to V1 that is provided in the brain [48,49], and raises the question whether the proposed neural sampling model could provide a possible mechanism (under the modelling assumptions made above) for cortical processing of such numerically weak external inputs.
The distribution of the resulting dominance durations, i.e., the time between perceptual switches, for the previously described setup with ambiguous input is shown for a continuous run of 10 4 s in Figure 4C (a similar method as in [18] was used to measure dominance durations, see Methods). This distribution can be approximated quite well by a Gamma distribution, which also provides a good fit to experimental data (see the discussion in [18]). We expect that also other features of the more abstract MCMC model for biological vision of [17,18], such as contextual biases and traveling waves, will emerge in larger and more detailed implementations of the MCMC approach through the proposed neural sampling method in networks of spiking neurons.

Discussion
We have presented a spiking neural network that samples from a given probability distribution via its inherent network dynamics. In particular the network is able to carry out probabilistic inference through sampling. The model, based on assumptions about the underlying probability distribution (formalized by the neural computability condition) as well as on certain assumptions regarding the underlying MCMC model, provides one possible neural implementation of the ''inference-by-sampling paradigm'' emerging in computational neuroscience.
During inference the observations (i.e., the variables which we wish to condition on) are modeled in this study by clamping the corresponding neurons by strong external input to the observed binary value. Units which receive no input or input with vanishing contrast (stimulus intensity) are treated as unobserved. Using this admittedly quite simplistic model of the input, we observed in simulations that our network model exhibits the following property: The onset of a sensory stimulus reduces the variability of the firing activity, which represents (after stimulus onset) a conditional distribution, rather than the prior distribution (see the difference between panels D and E of Figure 5. It is tempting to compare these results to the experimental finding of reduced firing rate variability after stimulus onset observed in several cortical areas [50]. We wish to point out however, that a consistent treatment of zero contrast stimuli requires more thorough modelling efforts (e.g., by explicitly adding a random variable for the stimulus intensity [35,51]), which is not the focus of the presented work.
Virtually all high-level computational tasks that a brain has to solve can be formalized as optimization problems, that take into account a (possibly large) number of soft or hard constraints. In typical applications of probabilistic inference in science and engineering (see e.g. [52,53]) such constraints are encoded in e.g., conditional probability tables or factors. In a biological setup they could possibly be encoded through the synaptic weights of a recurrent network of spiking neurons. The solution of such optimizations problems in a probabilistic framework via sampling, as implemented in our model, provides an alternative to deterministic solutions, as traditionally implemented in neural networks (see, e.g., [54] for the case of constraint satisfaction problems). Whereas an attractor neural network converges to one (possibly approximate) solution of the problem, a stochastic network may alternate between different approximate solutions and stay the longest at those approximate solutions that provide the best fit. This might be advantageous, as given more time a stochastic network can explore more of the state space and avoid shallow local minima. Responses to ambiguous sensory stimuli [44][45][46][47] might be interpreted as an optimization with soft constraints. The interpretation of human thinking as sampling process solving an inference task, recently proposed in cognitive science [28,55,56], further emphasizes that considering neural activity as an inferential process via sampling promises to be a fruitful approach.
Our approach builds on, and extends, previous work where recurrent networks of non-spiking stochastic neurons (commonly considered in artificial neural networks) were shown to be able to carry out probabilistic inference through Gibbs sampling [36]. In [57] a first extension of this approach to a network of recurrently connected spiking neurons had been presented. The dynamics of the recurrently connected spiking neurons are described as stepwise sampling from the posterior of a temporal Restricted Boltzmann Machine (tRBM) by introducing a clever interpretation of the temporal spike code as time varying parameters of a multivariate Gaussian distribution. Drawing one sample from the posterior of a RBM is, by construction, a trivial one-step task. In contrast to our model, the model of [57] does not produce multiple samples from a fixed posterior distribution, given the fixed input, but produces exactly one sample consisting of the temporal sequence of the hidden nodes, given a temporal input sequence. Similar temporal models, sometimes called Bayesian filtering, also underlie the important contributions of [58] and [32]. In [32] every single neuron is described as hidden Markov Model (HMM) with two states. Instead of drawing samples from the instantaneous posterior distribution using stochastic spikes, [32] presents a deterministic spike generation with the intention to convey the analog probability value rather than discrete samples. The approach presented here can be interpreted as a biologically more realistic version of Gibbs sampling for a specific class of probability distributions by taking into account a spike-based communication, finite duration PSPs and refractory mechanisms. Other implementations based on different distributions (e.g., directed graphical models) and different sampling methods (e.g., reversible MCMC methods) are of course conceivable and worth exploring.
In a computer experiment (see Figure 4, we used our proposed network to model aspects of biological vision as probabilistic inference along the lines of argumentation put forward in [16][17][18]. Our model was chosen to be quite simplistic, just to demonstrate that a number of experimental data on the dynamics of spontaneous activity [51,59,60] and binocular rivalry [44][45][46][47] can in principle be captured by this approach. The main point of the modelling study is to show that rather realistic neural dynamics can support computational functions rigorously formalized as inference via sampling.
We have also presented a model of spiking dynamics in continuous time that performs sampling from a given probability distribution. Although computer simulations of biological networks of neurons often actually use discrete time, it is desirable to also have a sound approach for understanding and describing the network sampling dynamics in continuous time, as the latter is arguable a natural framework for describing temporal processes in biology. Furthermore comparison to many existing continuous time neuron and network models of neurons is facilitated.
We have made various simplifying assumption regarding neural processes, e.g., simple symbolic postsynaptic potentials in the form of step-functions (reminiscent of plateau potentials caused by dendritic NMDA spikes [61]). More accurate models for neurons have to integrate a multitude of time constants that represent different temporal processes on the physical, molecular, and genetic level. Hence the open problem arises, to which extent this multitude of time constants and other complex dynamics can be integrated into theoretical models of neural sampling. We have gone one first step in this direction by showing that in computer simulations the two temporal processes that we have considered (refractory processes and postsynaptic potentials) can approximately be decoupled. Furthermore, we have presented simulation results suggesting that more realistic alpha-shaped, additive EPSPs are compatible with the functionality of the proposed network model.
Finally, we want to point out that the prospect of using networks of spiking neurons for probabilistic inference via sampling suggests new applications for energy-efficient spike-based and massively parallel electronic hardware that is currently under development [62,63].

Methods
We first provide details and proofs for the neural sampling models, followed by details for the computer simulations. Then we investigate typical firing statistics of individual neurons during neural sampling and examine the approximation quality of neural sampling with different neuron and synapse models.

Mathematical details
Notation. To keep the derivations in a compact form, we introduce the following notations. We define the function f w0 Details to neural sampling with absolute refractory period in discrete time. The following Lemmata 1 -3 provide a proof of Theorem 1. For completeness we begin this paragraph with a recapitulation of the definitions stated in Results. We then identify some central properties of the joint probability distribution p(f,z) and proof that the proposed network samples from the desired invariant distribution. For a given distribution p(z) over the binary variables z[f0,1g K with Vz[f0,1g K p(z)=0, the joint distribution over (f,z) with f[f0,1, . . . ,tg K is defined in the following way (see equation 7): The assumption p(z)=0 for all z is required to show the irreducibility of the Markov chain, a prerequisite to ensure the uniqueness of the invariant distribution of the MCMC dynamics. Furthermore, for the given distribution p(z) we define the functions u k : f0,1g K{1 ?R for k[f1, . . . ,Kg which map z \k .u k (z \k ): Instead of u k (z \k ) we simply write u k in the following. Lemma 1. The distribution p(f,z) has conditional distributions of the following form: These results can also be written more compactly in the following form: Here we use the fact that the logistic function s is the inverse of the logit function, i.e., p(z k~1 jz \k )~s(u k ).
This also shows that f k is independent from f \k given z \k , i.e., p(f k jf \k ,z \k )~p(f k jz \k ). Now we show the second relation using Bayes' rule: In order to facilitate the verification of the next two Lemmata, we first restate the definition of the operators T k in a more concise way: where u' k :~u k (z' \k )~logit(p(z k~1 jz' \k )). Lemma 2. For all k~1, . . . ,K the operator T k (f k jf' k ,z' \k ) leaves the conditional distribution p(f k jz' \k ) invariant.
Proof. For sake of simplicity, denote T k (f k~i jf' k~j ,z' \k )~T k ij for i,j[f0,1, . . . ,tg and p(f k~i jz' \k )~p i . We have to show p i! P t j~0 T k ij p j for i[f0,1, . . . ,tg. First we show p t~P t j~0 T k tj p j using p 0~1 {s(u k ) and p 1~p2~. . .~p t~s (u k )t {1 (which results from Lemma 1): Here we used the definition of the logistic function Here we used 1{s(x)~s({x). It is trivial to show p i~P t j~0 T k ij p j for i~1, . . . ,t{1 as P t j~0 T k ij p j~T k i,iz1 p iz1~piz1~pi . Here we used the facts that T k i,iz1~1 and p i~piz1 for i~1, . . . ,t{1 by definition. Lemma 3. For all k~1, . . . ,K the operator T k (z,fjf',z') leaves the distribution p(f,z) invariant.
Proof. We start from Lemma 2, which states that Here we used the relations d(z k ,f w0 k )~p(z k jf k ) and p(f k ,z k jz \k )~p(z k jf k )p(f k jz \k ) as well as p(f k jz \k )~p(f k jf \k ,z \k ) which directly follow from the definitions of T k (f,z,jf',z') and p(f,z).
Finally, we can verify that the composed operator T~T 1 0 . . . 0T K samples from the given distribution p.
Theorem 1. p(f,z) is the unique invariant distribution of operator T.
Proof. As all T k leave p(f,z) invariant, so does the concatenation T~T 1 0 . . . 0T K . To ensure that p(f,z) is the unique invariant distribution, we have to show that T is irreducible and aperiodic. T is aperiodic as the transition probabilities T k 00~1 { s(u k { log t)w0 and T k 00 v1 (this follows from the assumption Vz p(z)=0 made above).
The operator T is also irreducible for the following reason. First we see that from any state (f',z') in at most t steps we can get to the zero-state (f,z)~0 2K (and stay there) with non-zero probability, as T k i,iz1~1 for i~1, . . . ,t{1 and T k 01~1 {s(u k { log t)w0. Furthermore, it can be seen that any state (f f,ẑ z) can be reached from the zero-state (f,z)~0 2K in at most t steps since T k N0~s (u k { log t)w0 for any value of u k . Hence every final state (f f,ẑ z) can be reached from every starting state (f',z') in at most 2t steps with non-vanishing probability.
Details to neural sampling with a relative refractory period in discrete time. We augment the neuron model with a relative refractory period described by a function g(f k ). We first ensure existence of the corresponding function f (u k ). Based on these functions we then introduce the transition operator T of the Markov chain. This operator is shown to entail correct ''local'' computations.  (g 1 , . . . ,g t )[(R z 0 ) t be a tuple of non-negative real numbers, with g t~0 and at least one element g i §1. This defines the refractory function via g(f k ) :~g f k . There exists a unique C ? function f : R?(0,1) with the following property Vu[R: Furthermore, the function f has the property: Ai[f1, . . . ,tg Vu[R : 0vg i f (u)v1: Proof. Let g max :~max j[f1,...tg g j ; we know that g max §1. We define the function F : (0,1=g max )?R z : We can see that F is a positive C ? function on (0,1=g max ). Furthermore, F (x)=x is defined as a sum of functions of the form . Each factor 1=(1{g j x) is positive and strictly monotonous. Therefore, F is strictly monotonous on (0,1=g max ) with the limits: From here on, with the letter f we will denote the function characterized by the above Lemma for the given tuple g (which denotes the chosen refractory function). Definition 1. Define g 0~1 . The transition operator T k is defined in the following way for all k~1, . . . ,K: with u k~uk (z' \k ). Lemma 5. For all k~1, . . . ,K the unique invariant distribution q Ã (z k ,f k jf' \k ,z' \k ) of the operator T k (z k ,f k jf',z') fulfills P f k q Ã (z k ,f k jf' \k ,z' \k )~p(z k jz' \k ). This means, for a constant configuration z' \k , the operator T k produces samples z Ã k from the correct conditional distribution p(z k jz' \k ).
Proof. We define: where the function h(f k jz' \k ) is defined as: It is trivial to see that q Ã has the correct marginal distribution over z k : We now show that q Ã is the unique invariant distribution of T k . Because of the definition of T k , we only have to show that We denote q Ã (f k~i jz' \k )~: q i and T k (f k~i jf' k~j ,z' \k )~: T ij , i.e., we have to show Vi[f0,1, . . . ,tg : q i~P j T ij q j .
It is trivial to show q i~P j T ij q j for 1ƒiƒt{1, as there is only one non-vanishing element of transition operator, namely T i,iz1 : Here we used q i~h (f k~i jz \k )s(u k ) for iw0 and the definition of h(f k jz \k ). Now we show q 0~P j T 0j q j starting from equation (15) and additionally using the relations exp (u k )~s(u k )=(1{s(u k )) and q 0~1 {s(u k ) as well as the definition of q 1 . We define for the sake of simplicity y :~P t a~1 P t j~az1 (1{g j f (u k )): (1{f (u k ))(1{s(u k ))zf (u k )(1{s(u k ))! q 0 : We finally show q t~P j T tj q j , using the definition of q t~s (u k )h(f k~t jz \k )~s (u k ) y : The argument that the transition operator T k is aperiodic and irreducible is similar to the one presented in Lemma 1.
Details to neural sampling with an absolute refractory period in continuous time. In contrast to the discrete time model we define the state space of f k to be R z |½{2E,{E for Ew0, i.e., as the union of the positive real numbers and a small interval ½{2E,{E. We will define the sampling operator in such a way that after neuron k was refractory for exactly its refractory period t, its refractory variable f k is uniformly placed in the small interval ½{2E,{E, which represents now the resting state and replaces f k~0 . This avoids point measures (Dirac's Delta) on the value f k~0 . This system is still exactly equivalent to the system discussed in the main paper, as all spike-transition probabilities of T for f k v0 are constant. Hence, it does not matter which values f k assumes with respect to the spike mechanism during its nonrefractory period as long as f k v0. Definition 2. For a given distribution p(z) over the binary variables z[f0,1g K with Vz[f0,1g K p(z)=0, we define a joint distribution over (f,z) with f[R K in the following way: where I E :~½{2E,{E is the refractory resting state interval. In accordance with this definition we can also write The distribution p(f,z) has the following marginal distribution: Definition 3. For k[f1, . . . ,Kg and x[R the operator T k x is defined in the following way for a function q : R?R: where the functional F is defined as the one-sided limit from above at 0: The operator T is defined in the following way for a probability distribution q(f) on R K : where q(f 1 , . . . ,f k{1 , : ,f kz1 ,f K ) : R?R denotes the function q(f) of f k where f \k is held constant and u k :~u k (f w0 \k ). The transition operator T defines the following Fokker-Planck equation for a time-dependent distribution q t (f): The jump and drift functions W k (fjf') and A k (f) associated to the operator T are given by: Lemma 7. The operator T k uk leaves the conditional distribution p(f k jf \k ) invariant with u k~uk (f w0 \k ), i.e.: Proof. This is easy to proof using calculus and the relations L f k x R z (f k )~d(f k ) and F (p( : jf \k ))~s(u k )~exp (u k )(1{s(u k )). The Lemma follows then from the definition of T :~P k T k uk .
Details to neural sampling with a relative refractory period in continuous time. As already assumed in the case of the absolute refractory sampler in continuous time, we define the state space of f k to be R z |½{2E,{E for Ew0. Lemma 9. Let g be a continuous, non-negative function g : ½0,1?R z 0 with g(f k )~1 for f k ƒ0. There exists a unique C ? function f : R?R z with the following property Vu[R: Proof. We define the function F : R z 0 ?R in the following way: where a(r) :~Ð Hence F is strictly monotonously increasing. Furthermore, the following relations hold: Therefore the equation: has exactly one solution f (u) with F (f (u))~exp (u) in R z . From applying the implicit function theorem to F (x,u) :~F (x){ exp (u) it follows that f is C ? .
Definition 4. For all k[f1, . . . ,Kg and x[R the operator T k x is defined in the following way for a function q : R?R: The transition operator T k x defines the following Fokker-Planck equation for a time-dependent distribution q t (f k ): The jump and drift functions W k (f k jf k' ) and A k (f k ) associated to the operator T k x are given by: We define the distribution q Ã (f k jz \k ) as: where a(f k ) :~Ð 1 0 g(f' k )df' k . By applying the operator T k u k to q Ã one can verify that T k u k q Ã~0 holds using the definition of f (u k ) given in (16). Furthermore we can compute the ratio:

Details to the computer simulations
The simulation results shown in Figure 2, Figure 3 and Figure 4 used the biologically more realistic neuron model with the relative refractory mechanism. During all experiments the first second of simulated time was discarded as burn-in time. The full list of parameters defining the experimental setup is given in Table 1. All occurring joint probability distributions are Boltzmann distributions of the form given in equation (5). Example Python [64] scripts for neural sampling from Boltzmann distributions are available on request and will be provided on our webpage. The example code comprises networks with both absolute and relative refractory mechanism. It requires standard Python packages only and is readily executable.  Table 1. Panel (C) shows the corresponding functions f (u), which result from numerically solving equation (11). The spike patterns in panel (D) show the response of the neurons when the membrane potential is low (u k~{ 1 for 0vtv250 ms) or high (u k~z 2 for 250 msvtv 500 ms). These membrane potentials encode p(z k~1 )~0:269 and p(z k~1 )~0:881, respectively according to (3) and (4). The binary state z k~1 is indicated by gray shaded areas of duration t : dt~20 ms after each spike. Figure 3: Sampling from a Boltzmann distribution by spiking neurons with relative refractory mechanism. We examined the spike response of a network of 40 randomly connected neurons which sampled from a Boltzmann distribution. The excitabilities b k as well as the synaptic weights W ki (~W ik ) were drawn from Gaussian distributions (with diagonal elements W ii~0 ). For the full list of parameters please refer to Table 1 As stated before, the neuron model with relative refractory mechanism g k (f) does not entail the correct overall invariant distribution p(z). To estimate the impact of this approximation on the joint network dynamics, we compared the distribution p(z 24 , . . . ,z 28 ) over five neurons (indicated by gray background in A) in the spiking network with the correct distribution obtained from Gibbs sampling. The probabilities were estimated from 10 7 samples. A more quantitative analysis of the approximation quality of neural sampling with a relative refractory mechanism is provided below. Figure 4: Modeling perceptual multistability as probabilistic inference with neural sampling. We demonstrate probabilistic inference and learning in a network of orientation selective neurons. As a simple model we consider a network of 217 neurons on a hexagonal grid as shown in panel (F). Any two neurons with distance ƒ8 were synaptically connected (neighboring units had distance 1). For the remaining parameters of the network and neuron model please refer to Table 1. Each neuron featured a p-periodic tuning curve as depicted in panel (B): with base sensitivity v 0 , contrast C, peakedness k and preferred orientation Q Q k . The preferred orientations Q Q k of the neurons were chosen to cover the entire interval ½0,p) of possible orientations with equal spacing and were randomly assigned to the neurons.

Details to
For simplicity we did not incorporate the input dynamics in our probabilistic model, but rather trained the network directly like a fully visible Boltzmann machine. We used for this purpose a standard Boltzmann machine learning rule known as contrastive divergence [41,65]. This learning rule requires posterior samplesz z, i.e., network states under the influence of the present input, and approximate prior samples z ? , which reflect the probability distribution of the network in the absence of stimuli. The update rules for synaptic weights and neuronal excitabilities read: g ki~g if n k and n i are connected While more elaborate policies can speed up convergence, we simply used a global learning rate g which was constant in time.
The values of W ki and b k were initialized at 0. We generated binary training patterns in the following way: 1. A global orientation Q was drawn uniformly from ½0,p), 2. each neuron was independently set to be active with probability p(z k~1 )~V k (Q), 3. the resulting network statez z was taken as posterior sample.
To obtain an approximate prior sample z ? we let the network run for a short time freely starting from (f f,z z). The variablesz z were also assumed to be observed withf f k * iid. uniformly in f1, . . . ,tg ifz z k~1 andf f k~0 otherwise. After evolving freely for 20 time steps, the resulting network state z ? was taken as approximate prior sample and W and b were updated according to (19). This process was repeated N train~1 0 5 times. As a result, neurons with similar preferred orientations featured excitatory synaptic connections (W ki~6 :4 : 10 {3 +6:7 : 10 {3 = mean + standard deviation of weight distribution), those with dissimilar orientations maintained inhibitory synapses (W ki~{ 4:9 : 10 {3 +5:2 : 10 {3 ). Here, preferred orientations Q Q i and Q Q j are defined as similar if V i ( Q Q j ){v 0~Vj ( Q Q i ){v 0 w0:5 C, otherwise they are dissimilar. Neuronal biases converged to b k~{ 0:08+0:03.
We illustrate the learned prior distribution p(z) of the network through sampled states when the network evolved freely. As seen in panel (D), the population vector -a 2-dimensional projection of the high dimensional network state -typically reflected an arbitrary, yet coherent, orientation (for the definition of the population vector see below). Each dot represents a sampled network state z.
To apply an ambiguous cue, we clamped 8 out of 217 neurons: Two units with Q Q k &p=4 and two with Q Q k &3p=4 were set active, two units with Q Q k &0 and two with Q Q k &p=2 were set inactive. This led to a bimodal posterior distribution as shown in panel (E). The sampling network represented this distribution by encoding either global perception separately: The trace of network states z(t) roamed in one mode for multiple steps before quickly crossing the state space towards the opposite percept.
We define the population vector x of a network state z as a function of the preferred orientations of all active units: This definition of x is not based on the preferred orientations Q Q k which are used for generating external input to the network from a given stimulus with orientation Q. It is rather based on the preferred orientationsQ Q k measured from the network response. We used population vector decoding based on the measured valuesQ Q k , as they are conceptually closer to experimentally measurable preferred orientations, and this decoding hence does not require knowledge of the (unobservable) Q Q k . For every neuron n k the preferred orientationQ Q k was measured in the following way. We estimated a tuning curveṼ V k (Q) by a van-Mises fit (of the form (18)) to data from stimulation trials in which neuron n k was not clamped, i.e., where n k was only stimulated by recurrent input (feedforward input was modeled by clamping 8 out of 217 neurons as a function of stimulus orientation Q as before). Due to the structured recurrent weights, the experimentally measured tuning curvesṼ V k (Q) were found to be reasonably close to the tuning curves V k (Q) used for external stimulation.Q Q k was set to the preferred orientation ofṼ V k (Q) (localization parameter of the van-Mises fit). The measured valuesQ Q k turned out to be consistent with the preferred orientations Q Q k ( Q Q k {Q Q k~6 : 10 {4 +8:3 : 10 {3 averaged over all K neurons). The mean and standard deviation of the remaining parameter values v 0 , C and k of the fitted tuning curves V V k (Q) are listed in Table 1 next to the ones used for stimulation.
The population vector x was defined in (20) with the argument 2Q Q k (instead ofQ Q k ) as orthogonal orientations should cancel each other and neighborhood relations should be respected. For example neurons withQ Q k~E andQ Q k~p {E contribute similarly to the population vector for small e. But counter to intuition the population vector of a state z with dominant orientation Q z will point into direction Q x~2 Q z . For visualization in panel (D) and (E) we therefore rescaled the population vector: If (x 0 ,x p=4 ).(r x ,Q x ) in polar coordinates, then the dot is located at (r x ,Q x =2) in accord with intuition. The black semicircles equal jxj~r x~4 5.
The population vector (x 0 ,x p=4 )[R 2 was also used for measuring the dominance durations shown in panel (C). To this R 2 was divided into 3 areas: (a) x p=4 v{35, (b) {35ƒx p=4 ƒ35, (c) 35vx p=4 . We detected a perceptual switch when the network state entered area (a) or (c) while the previous perception was (c) or (a), respectively.
In panel (F) neurons n k with z k~1 are plotted with their preferred orientation color code, inactive neurons are displayed in white. Cells marked by a dot (.) were part of the observed variables o. The three network states correspond to z(t i ) with t 1~1 00 ms, t 2~2 50 ms and t 3~4 00 ms in the spike pattern in panel (G). The spike pattern shows the response of the freely evolving units around a perceptual switch during sampling from the posterior distribution. The corresponding trace of the population vector is drawn as black line in panel (E). The width of the light-gray shaded areas in the spike pattern equals the PSP duration t : dt, i.e., neurons that spiked in these intervals were active in the corresponding state in (F).

Firing statistics of neural sampling networks
In previous sections it was shown that a spiking neural network can draw samples from a given joint distribution which is in a welldefined class of probability distributions (see the neural computability condition (4)). Here, we examine some statistics of individual neurons in a sampling network which are commonly used to analyze experimental data from recordings. The spike trains and membrane potential data are taken from the simulation presented in Figure 3. Figure 5A,B exemplarily show the distribution of the membrane potential u k and the interspike interval (ISI) histogram of a single neuron, namely neuron n 26 which was already considered in Figure 3B. The responses of other neurons yield qualitatively similar statistics. The bell-shaped distribution of the membrane potential is commonly observed in neurons embedded in an active network [66]. The ISI histogram reflects the reduced spiking probability immediately after an action potential due the refractory mechanism. Interspike intervals larger than the refractory time constant t : dt~20 ms roughly follow an exponential distribution. Similar ISI distributions were observed during invivo recordings in awake, behaving monkeys [67]. Figure 5C shows a scatterplot of the coefficient of variation (CV) of the ISIs versus the average ISI for each neuron in the network. The neurons exhibited a variety of average firing rates between 3:5 Hz and 31:5 Hz. Most of the neurons responded in a highly irregular manner with a CV &1. Neurons with high firing rates had a slightly lower CV due to the increased influence of the refractory mechanism The dashed line marks the CV of a Poisson process, i.e., a memoryless spiking behavior. The CV of neuron n 26 is marked by a cross. The structure of this plot resembles, e.g., data from recordings in behaving macaque monkeys [68] (but note the lower average firing rate).

Approximation quality of neural sampling with different neuron and synapse models
The theory of the neuron model with absolute refractory mechanism guarantees sampling form the correct distribution. In contrast, the theory for the neuron model with a relative refractory mechanism only shows that the sampling process is ''locally correct'', i.e., that it would yield correct conditional distributions p(z k jz \k ) for each individual neuron if the state of the remaining network z \k stayed constant. Therefore, the stationary distribution of the sampling process with relative refractory mechanism only provides an approximation to the target distribution. In the following we examine the approximation quality and robustness of sampling networks with different refractory mechanisms for target Boltzmann distributions with parameters randomly drawn from different distributions. Furthermore, we investigate the effect of additive PSP shapes with more realistic time courses.
We generated target Boltzmann distributions with randomly drawn weights W ki and biases (excitabilities) b k and computed the similarity between these reference distributions and the corresponding neural sampling approximations. The setup of these simulations is the same as for the simulation presented in Figure 3. As we aimed to compare the distribution q Ã (z) sampled by the network with the exact Boltzmann distribution p(z), we reduced the number of neurons per network to K~10. This resulted in a state space of 2 10 possible network states z for which the normalization constant for the target Boltzmann distribution could be computed exactly. The weight matrix W was constraint to be symmetric with vanishing diagonal. Off-diagonal elements were drawn from zero-mean normal distributions with three different standard deviations s~0:03, s~0:3 and s~3, whereas the b k were sampled from the same distribution as in Figure 3. For every value of the hyperparameter s we generated 100 random distributions. For Boltzmann distributions with small weights (s~0:03), the RVs are nearly independent, whereas distributions with intermediate weights (s~0:3) show substantial statistical dependencies between RVs. For very large weights (s~3), the probability mass of the distributions is concentrated on very few states (usually 90% on less than 10 out of the 2 10 states). Hence, the range of the hyperparameter 0:03ƒsƒ3 considered here covers a range a very different distributions.
The approximation quality of the sampled distribution was measured in terms of the Kullback-Leibler divergence between the target distribution p and the neural approximation q Ã We estimated q Ã from 10 7 samples for each simulation trial using a Laplace estimator, i.e., we added a priori 1 to the number of occurrences of each state z. Table 2 shows the means and the standard deviations of the Kullback-Leibler divergences between the target Boltzmann distributions and the estimated approximations stemming from neural sampling networks with three different neuron and synapse models: the exact model with absolute refractory mechanism and two models with different relative refractory mechanisms shown in the bottom and middle row in Figure 2B. Additionally, as a reference, we provide the (analytically calculated) Kullback-Leibler divergences for fully factorized distributions, i.e., q Ã (z)~P k q Ã (z k ) with correct marginals q Ã (z k )~p(z k ) but independent variables z i , z j for i=j.
The absolute refractory model provides the best results as we expected due to the theoretical guarantee to sample from the correct distribution (the non-zero Kullback-Leibler divergence is caused by the estimation from a finite number of samples). The models with relative refractory mechanism provide faithful approximations for all values of the hyperparameter s considered here. These relative refractory models are characterized by the theory to be ''locally correct'' and turn out to be much more  Mean and standard deviation of the Kullback-Leibler divergence D KL (pjjq Ã ) between reference Boltzmann distributions p and neural sampling approximations q Ã for three different neuron models (corresponding to columns) and three different values for the reference distribution hyperparameter s (corresponding to rows). The parameter s controls the standard deviation of the weights of the reference distributions p(z). In case of very strong synaptic interactions (leading to sharply peaked distributions, s~3) the approximation quality of the spiking network degrades, if the neurons feature a relative refractory mechanism. The data was computed from 100 randomly generated Boltzmann distributions and their neural approximations for each value of s. doi:10.1371/journal.pcbi.1002211.t002 accurate approximations than fully factorized distributions if substantial statistical dependencies between the RVs are present (i.e., s~0:3, s~3). As expected, a late recovery of the refractory function g(f) is beneficial for the approximation quality of the model as it is closer to an absolute refractory mechanism. Figure 6 shows the full histograms of the Kullback-Leibler divergences for the intermediate weights group (s~0:3). Systematic deviations due to the relative refractory mechanism are on the same order as the effect of estimating from finite samples (as can be seen, e.g., from a comparison with the absolute refractory model which has 0 systematic error). For completeness, we mention that the divergences of the fully factorized distributions of 2 out of the 100 networks with D KL w0:1 are not shown in the plot. The theorems presented in this article assumed renewed (i.e., non-additive), rectangular PSPs. In the following we examine the effect of additive PSPs with more realistic time courses. We define additive, alpha-shaped PSPs in the following way. The influence Du ki of each presynaptic neuron n i on the postsynaptic membrane potential u k is modeled by convolving the input spikes with a kernel k: where k(s)~l : (e {s=t z {e {s=t{ ) for s §0 and k(s)~0 for sv0, and t f i for f [N are the spike times of the presynaptic neuron n i . The time constant governing the rising edge of the PSPs was set to t {~3 ms. The time constant controlling the falling edge was chosen equal to the duration of rectangular PSPs, t z~t : dt~20 ms. The scaling parameter l was set such that the time integral over a single PSP matches the time integral over the theoretically optimal rectangular PSP, i.e., lt : dt=(t z {t { )~20=17. These parameters display a simple and reasonable choice for the purpose of this study (an optimization of l, t z and t { is likely to yield an improved approximation quality). Figure 7A shows the resulting shape of the non-rectangular PSP. Furthermore the time course of the function g(f k (t)) caused by a single spike of neuron n k is shown in order to illustrate that the time constants of g and of a PSP are closely related due to the Figure 6. Comparison of neural sampling with different neuron and synapse models. The figure shows a histogram of the Kullback-Leibler divergence between 100 different Boltzmann distributions over K = 10 variables (with parameters randomly drawn, see setup of Figure 3) and approximations stemming from different neural sampling networks. Networks with absolute refractory mechanism provide the best approximation (as expected from theoretical guarantees). Networks consisting of neurons with relative refractory mechanisms, with only ''locally'' correct sampling, also provide a close fit to the true distribution (see inset) compared to a fully factorized approximation (assuming correct marginals and independent variables). Furthermore, it can be seen that sampling networks with more realistic, alpha-shaped, additive PSPs still fit the true distribution reasonably well. doi:10.1371/journal.pcbi.1002211.g006 Figure 7. Sampling from a Boltzmann distribution with more realistic PSP shapes. (A) The upper panel shows the shape of a single PSP elicited at time t~0. The lower panel shows the time course of the refractory function g(f k (t)) caused by a single spike of neuron n k at t~0. The greyshaded area of length t : dt~20ms indicates the interval of neuron n k being active (i.e., z k~1 ) due to a single spike of neuron n k at time t~0. (B) Shown is the probability distribution of 5 out of 40 neurons. The plot is similar to Figure 3C, however it is generated with a sampling network that features alpha-shaped, additive PSPs. It can be seen that the network still produces a reasonable approximation to the true Boltzmann distribution (determined by Gibbs sampling). doi:10.1371/journal.pcbi.1002211.g007 assumption t z~t : dt made above. Preliminary and non-exhaustive simulations seem to suggest that the choice t z~t : dt yields better approximation quality than setting t z &t : dt or t z %t : dt; however it is very well possible that a mismatch between t z and t : dt can be compensated for by adapting other parameters, e.g., the PSP magnitude or a specific choice of the refractory function g. Figure 7B shows the results of an experiment, similar to the one presented in Figure 3C , with additive, alpha-shaped PSPs and relative refractory mechanism. While differences to Gibbs sampling results are visible, the spiking network still captures dependencies between the binary random variables quite well.
For a quantitative analysis of the approximation quality, we repeated the experiment of Figure 6 with additive, alpha-shaped PSPs (shown as green bars). The Kullback-Leibler divergence D KL (pjjq Ã ) to the true distribution is clearly higher compared to the case of renewed, rectangular PSPs. Still networks with this more realistic synapse model account for dependencies between the random variables z and yield a better approximation of p(z) than fully factorized distributions.