Probabilistic models for neural populations that naturally capture global coupling and criticality
We represent the response of a neural population with a binary vector s = {s1 , s2 , … , sN} ∈ {0 , 1}N identifying which of the N neurons elicited at least one action potential ( ‘1’ ) and which stayed silent ( ‘0’ ) during a short time window ., Our goal is to build a model for the probability distribution of activity patterns , p, ( s ) , given a limited number M of samples , D = { s ( 1 ) , … , s ( M ) } , observed in a typical recording session ., The regime we are mainly interested in is the one where the dimensionality of the problem is sufficiently high that the distribution p cannot be directly sampled from data , i . e . , when 2N ≫ M . Note that we are looking to infer models for the unconditional distribution over neural activity patterns ( i . e . , the population “vocabulary” ) , explored in a number of recent papers 8 , 9 , 11 , 13–18 , 24 , 34 , rather than to construct stimulus-conditional models ( i . e . , the “encoding models” , which have a long tradition in computational neuroscience 1–3 ) ., Previous approaches to modeling globally coupled populations focused on the total network activity , also known as synchrony , K ( s ) = ∑ i = 1 N s i ., The importance of this quantity was first analyzed in the context of probabilistic models in Ref 11 where the authors showed that a K-pairwise model , which generalizes a pairwise maximum entropy model by placing constraints on the statistics of K, ( s ) , is much better at explaining the observed population responses of 100+ salamander retinal ganglion cells than a pairwise model ., Specifically , a pairwise model assumes that the covariance matrix between single neuron responses , Cij = 〈sisj〉 , which can be determined empirically from data D , is sufficient to estimate the probability of any population activity pattern ., In the maximum entropy framework , this probability is given by the most unstructured ( or random ) distribution that reproduces exactly the measured Cij:, p ( s ; J ) = 1 Z ( J ) exp ( ∑ i , j = 1 N J i j s i s j ) , ( 1 ), where Z ( J ) is a normalization constant , and J is a coupling matrix which is chosen so that samples from the model have the same covariance matrix as data ., Note that because s i 2 = s i , the diagonal terms Jii of the coupling matrix correspond to single neuron biases , i . e . firing probabilities in the absence of spikes from other neurons ( previous work 11 used a representation si ∈ {−1 , 1} for which the single neuron biases need to be included as separate parameters and where Jii are all 0 ) ., A K-pairwise model generalizes the pairwise model and has the form, p ( s ; J , ϕ ) = 1 Z ( J , ϕ ) exp ( ∑ i , j = 1 N J i j s i s j + ∑ k = 0 N ϕ k δ k , K ( s ) ) ., ( 2 ), The coupling matrix J has the same role as in a pairwise model while the additional parameters ϕ are chosen to match the probability distribution of K, ( s ) under the model to that estimated from data ., The “potentials” ϕk introduced into the K-pairwise probabilistic model , Eq ( 2 ) , globally couple the population , and cannot be reduced to low-order interactions between , e . g . , pairs or triplets , of neurons , except in very special cases ., We will generically refer to probabilistic models that impose non-trivial constraints on population-level statistics ( of which the distribution of total network activity K is one particular example ) as “globally coupled” models ., Here we introduce new semiparametric energy-based models that extend the notion of global coupling ., These models are defined as follows:, p ( s ; α , V ) = e - V ( E ( s ; α ) ) Z ( α , V ) , ( 3 ), where E ( s; α ) is some energy function parametrized by α , and V is an arbitrary increasing differentiable function which we will refer to simply as the “nonlinearity . ”, The parametrization of the energy function should be chosen so as to reflect local interactions among neurons ., Crucially , while it is necessary to choose a specific parametrization of the energy function , we do not make any assumptions on the shape of the nonlinearity—we let the shape be determined nonparametrically from data ., Fig 1 schematically displays the relationship between the previously studied probabilistic models of population activity and two semiparametric energy-based models that we focus on in this paper , the semiparametric independent model ( which we also refer to as “V ( independent ) ” ) and the semiparametric pairwise model ( which we also refer to as “V ( pairwise ) ” ) ., Our motivation for introducing the global coupling via the nonlinearity V traces back to the argument made in Ref 11 for choosing to constrain the statistics of synchrony , K, ( s ) ; in short , the key intuition in earlier work has been that K, ( s ) is a biologically relevant quantity which encodes information about the global state of a population ., There are , however , many other quantities whose distributions could contain signatures of global coupling in a population ., In particular , while most energy functions—e . g . , the pairwise energy function , E ( s; J ) = −∑i , j Jijsisj—are defined solely in terms of local interactions between small groups of neurons , the statistics of these same energy functions ( for instance , their moments ) are strongly shaped by global effects ., Specifically , we show in Methods that the role of the nonlinearity in Eq ( 3 ) is precisely to match the probability density of the energy under the model to that estimated from data ., In other words , once any energy function for Eq ( 3 ) has been chosen , the nonlinearity V will ensure that the distributions of that particular energy in the model and over data samples agree ., Constraining the statistics of the energy E ( s; α ) is different from constraining the statistics of K, ( s ) , used in previous work ., First , the energy depends on a priori unknown parameters α which must be learned from data ., Second , while K, ( s ) is always an integer between 0 and N , the energy can take up to 2N distinct values; this allows for extra richness but also requires us to constrain the ( smoothed ) histogram of energy rather than the probability of every possible energy value , to prevent overfitting ., As we discuss next , the statistics of the energy are also closely related to criticality , a formal , model-free property distinguishing large , globally-coupled neural populations ., The notion of criticality originates in thermodynamics where it encompasses several different properties of systems undergoing a second-order phase transition 35 ., Today , many other phenomena , such as power-law distributed sizes of “avalanches” in neural activity , have been termed critical 20 ., Our definition , which we discuss below , is a restricted version of the thermodynamic criticality ., We consider a sequence of probability distributions { p N } N = 1 ∞ over the responses of neural populations of increasing sizes , N . These probability distributions define the discrete random variable s ( the population response ) , but they can also be thought of simply as functions which map a population response to a number between 0 and 1 ., Combining these two viewpoints , we can consider a real-valued random variable pN, ( s ) ∈ ( 0 , 1 ) which is constructed by applying the function pN to the random variable s ., The behavior of this random variable as N → ∞ is often universal , meaning that some of its features are independent of the precise form of pN ., As is conventional , we work with the logarithm of pN, ( s ) instead of the actual distribution ., We call a population “critical” if the standard deviation of the random variable log pN, ( s ) /N does not vanish as the population size becomes large , i . e ., 1 N σ ( log p N ( s ) ) ↛ 0 as N → ∞ ., ( 4 ), ( For completeness , we further exclude some degenerate cases such as when the probability density of log pN, ( s ) /N converges to two equally sized delta functions ., ) The above definition is related to criticality as studied in statistical physics ., In thermodynamics , σ ( log p N ( s ) ) / N is proportional to the square root of the specific heat , which diverges in systems undergoing a second-order phase transition ., While at a thermodynamical critical point σ ( log pN, ( s ) ) /N scales as N−γ with γ ∈ ( 0 , 1/2 ) , here we are concerned with the extreme case of γ = 0 ., Rather than being related to second-order phase transitions , this definition of criticality is related to the so-called Zipf law 23 ., A pattern s can be assigned a rank by counting how many other patterns have a higher probability ., In its original form , a probability distribution is said to satisfy Zipf law if the probability of a pattern is inversely proportional to its rank ., No real probability distribution is actually expected to satisfy this definition precisely , but there is a weaker form of Zipf law which concerns very large populations , and which is much less restrictive ., This weaker form can be stated as a smoothed version of the original Zipf law ., Consider patterns whose rank is in some small interval r , r + ΔN , and denote pN, ( r ) the average probability of these patterns ., We generalize the notion of Zipf law to mean that for very large populations pN, ( r ) ∝ r−1 ( ΔN is assumed to go to zero sufficiently quickly with N ) ., As shown in Ref 23 , a system is critical in the sense of Eq ( 4 ) precisely when it follows this generalized Zipf law ., Practically speaking , no experimentally studied population ever has an infinite size , and a typical way to check for signs of criticality is to see if a log-log plot of a pattern probability versus its rank resembles a straight line with slope −1 ., Most systems are not expected to be critical ., The simplest example is a population of identical and independent neurons ,, p N ( s ) = q ∑ i = 1 N s i ( 1 - q ) N - ∑ i = 1 N s i , ( 5 ), where q is the probability of eliciting a spike ., For such population ,, 1 N σ ( log p N ( s ) ) = 1 N q ( 1 - q ) log q 1 - q , ( 6 ), which vanishes for very large number of neurons , and so the system is not critical ., More generally , if pN, ( s ) can be factorized into a product of probability distributions over smaller subpopulations which are independent of each other and whose number is proportional to N , then log pN, ( s ) /N turns into an empirical average whose standard deviation is expected to vanish in the large N limit , and the population is not critical ., Reversing this argument , signatures of criticality can be interpreted as evidence that the population is globally coupled , i . e . that it cannot be decomposed into independent parts ., These preliminaries establish a direct link between criticality and semiparametric energy models of Eq ( 3 ) ., Nonlinearity in semiparametric energy models makes sure that the statistics of the energy E ( s; α ) , and , since V ( E ) is monotone , also the statistics of log p ( s; α , V ) are modeled accurately ( see Methods ) ., Because the behavior of log probability is crucial for criticality , as argued above , semiparametric energy models can capture accurately and efficiently the relevant statistical structure of any system that exhibits signs of criticality and/or global coupling ., To fully specify semiparametric energy models , we need a procedure for constructing the nonlinearity V ( E ) ., We cannot let this function be arbitrary because then the model could learn to assign nonzero probabilities only to the samples in the dataset , and hence it would overfit ., To avoid such scenarios , we will restrict ourselves to functions which are increasing ., We also require V ( E ) to be differentiable so that we can utilize its derivatives when fitting the model to data ., The class of increasing differentiable functions is very large ., It includes functions as diverse as the sigmoid , 1/ ( 1 + exp ( −E ) ) , and the square root , E ( for positive E ) , but we do not want to restrict ourselves to any such particular form—we want to estimate V ( E ) nonparametrically ., Nonparametric estimation of monotone differentiable functions is a nontrivial yet very useful task ( for example , consider tracking the height of a child over time—the child is highly unlikely to shrink at any given time ) ., We follow Ref 36 and restrict ourselves to the class of strictly monotone twice differentiable functions for which V′′/V′ is square-integrable ., Any such function can be represented in terms of a square-integrable function W and two constants γ1 and γ2 as, V ( E ) = γ 1 + γ 2 ∫ E 0 E exp ( ∫ E 0 E ′ W ( E ′ ′ ) d E ′ ′ ) d E ′ , ( 7 ), where E0 is arbitrary and sets the constants to γ1 = V ( E0 ) , γ2 = V′ ( E0 ) ., The function is either everywhere increasing or everywhere decreasing ( depending on the sign of γ2 ) because the exponential is always positive ., Eq ( 7 ) is easier to understand by noting that V ( E ) is a solution to the differential equation V′′ = WV′ ., This means , for example , that on any interval on which W = 0 , the equation reduces to V′′ = 0 , and so V ( E ) is a linear function on this interval ., If V ( E ) is increasing ( V′ > 0 ) , it also shows that the sign of W at a given point determines the sign of the second derivative of V at that point ., An advantage of writing the nonlinearity in the form of Eq ( 7 ) is that we can parametrize it by expanding W in an arbitrary basis without imposing any constraints on the coefficients of the basis vectors yet V ( E ) is still guaranteed to be monotone and smooth ., In particular , we will use piecewise-constant functions for W . This allows us to use unconstrained optimization techniques for fitting our models to data ., We start by considering one of the simplest models of the form Eq ( 3 ) , the semiparametric independent model:, p ( s ; α , V ) = e - V ( - ∑ i = 1 N α i s i ) Z ( α , V ) ., ( 8 ), If V were a linear function , the model would reduce to an independent model , i . e . a population of independent neurons with diverse firing rates ., In general , however , V introduces interactions between the neurons that may not have a straightforward low-order representation ., When fitted to our data , the nonlinearity V turns out to be a concave function ( see later sections on more complex models for a detailed discussion of the shape of the nonlinearity ) ., Note that if V had a simple functional form such as a low order polynomial , then the model Eq ( 8 ) would be closely related to mean field models of ferromagnetism with heterogenous local magnetic field studied in physics ., Our first goal is to use this simple model to verify our intuition that the nonlinearity helps to capture criticality ., Many population patterns are observed several times during the course of the experiment , and so it is possible to estimate their probability simply by counting how often they occur in the data 19 ., Given this empirical distribution , we construct a corresponding Zipf plot—a scatter plot of the frequency of a pattern vs its rank ., For systems which are close to critical , this should yield a straight line with slope close to −1 on a log-log scale ., We repeat the same procedure with samples generated from a semiparametric independent model as well as an independent model , which were both fitted to the responses of all 160 neurons ., Fig 2 shows all three scatter plots ., The independent model vastly deviates from the empirical Zipf plot; specifically , it greatly underestimates the probabilities of the most likely states ., In contrast , the learned semiparametric independent model follows a similar trend to that observed in data ., This does not mean that the semiparametric independent model itself is an excellent model for the detailed structure in the data , but it is one of the simplest possible extensions of the trivial independent model that qualitatively captures both global coupling and the signatures of criticality ., Since the semiparametric independent model is able to capture the criticality of the data distribution , we also expect it to accurately model other features of the data which are related to the globally coupled nature of the population ., To verify this , Fig 3A compares the empirical probability distribution of the total activity of the population K ( s ) = ∑i si to that predicted by the semiparametric independent model ., The match is very accurate , especially when compared to the same distribution predicted by the independent model ., This result goes hand in hand with the analysis in 39 which showed that interactions of all orders ( in our case mediated by the nonlinearity ) are necessary to model the wide-spread distribution of the total activity ., The independent model is a maximum entropy model which constrains the mean responses , 〈si〉 , of all neurons ., In other words , neurons sampled from the model would have the same firing rates as those in the data ( up to sampling noise ) ., Even though the semiparametric independent model is strictly more general , it does not retain this property when the parameters α and the nonlinearity V are learned by maximizing the likelihood of data ., Fig 3B demonstrates this point: although the predicted firing rates are approximately correct , there are slight deviations ., On the other hand , the nonlinearity induces pairwise correlations between neurons which is something the independent model by construction cannot do ., Fig 3C compares these predicted pairwise correlations to their data estimates ., While there is some correlation between the predicted and observed covariances , the semiparametric independent model often underestimates the magnitude of the covariances and does not capture the fine details of their structure ( e . g . the largest covariance predicted by the semiparametric independent model is about 5× smaller than the largest covariance observed in the data ) ., This is because a combination of independent terms and a single nonlinearity does not have sufficient expressive power , motivating us to look for a richer model ., One way to augment the power of the semiparametric independent model that permits a clear comparison to previous work is by means of the semiparametric pairwise model:, p ( s ; J , V ) = 1 Z ( J , V ) exp ( - V ( - ∑ i , j = 1 N J i j s i s j ) ) ., ( 9 ), We fit this model to the responses of the various subpopulations of the 160 neurons , and we compare the resulting goodness-of-fit to that of a pairwise ( Eq ( 1 ) ) , K-pairwise ( Eq ( 2 ) ) , and semiparametric independent model ( Eq ( 8 ) ) ., We measure goodness-of-fit as the improvement of the log-likelihood of data per neuron under the model relative to the pairwise model , as shown in Fig 4A ., This measure reflects differences among models rather than differences among various subpopulations ., The semiparametric pairwise model consistently outperforms the other models and this difference grows with the population size ., To make sure that this improvement is not specific to this particular experiment , we also fitted the models to two additional recordings from the salamander retina which were also collected as part of the study 11 ., One consists of 120 neurons responding to 69 repeats of a 30 second random checkerboard stimulus , and the other of 111 neurons responding to 98 repeats of a 10 second random full-field flicker stimulus ., As shown in Fig 4B , the improvements of individual models on these datasets are consistent with the ones observed for the population stimulated with a natural movie ., The advantage of using likelihood as a goodness-of-fit measure is its universal applicability which , however , comes hand-in-hand with the difficulty of interpreting the quantitative likelihood differences between various models ., An alternative comparison measure that has more direct relevance to neuroscience asks about how well the activity of a single chosen neuron can be predicted from the activities of other neurons in the population ., Given any probabilistic model for the population response , we use Bayes rule to calculate the probability of the ith neuron spiking ( si = 1 ) or being silent ( si = 0 ) conditioned on the activity of the rest of the population ( s−i ) as, p ( s i | s - i ; α ) = p ( s ; α ) p ( s i = 1 , s - i ; α ) + p ( s i = 0 , s - i ; α ) ., ( 10 ), We turn this probabilistic prediction into a nonrandom one by choosing whether the neuron is more likely to spike or be silent given the rest of the population , i . e ., s i ( s - i ; α ) = argmax s i ∈ { 0 , 1 } p ( s i | s - i ; α ) ., ( 11 ), In Fig 4C and 4D we compare such predictive single neuron models constructed from semiparametric pairwise , K-pairwise , pairwise , and semiparametric independent models learned from the data for populations of various sizes ., Specifically , we ask how often these models would make a mistake in predicting whether a chosen single neuron has fired or not ., Every population response in our dataset corresponds to 20 ms of an experiment and so we can report this accuracy as number of errors per unit of time ., Predictions based on the semiparametric pairwise model are consistently the most accurate ., Fig 5A shows the nonlinearities of the semiparametric pairwise models that we learned from data ., In order to compare the nonlinearities inferred from populations of various sizes , we normalize the domain of the nonlinearity as well as its range by the number of neurons ., Even though the nonlinearities could have turned out to have e . g . a sigmoidal shape , the general trend is that they are concave functions whose curvature—and thus departure from the linear V that signifies no global coupling—grows with the population size ., The shape of these nonlinearities is reproducible over different subnetworks of the same size with very little variability ., To further visualize the increasing curvature , we extrapolated what these nonlinearities might look like if the size of the population was very large ( the black curve in Fig 5A ) ., This extrapolation was done by subtracting an offset from each curve so that V ( 0 ) = 0 , and then fitting a straight line to a plot of 1/N vs . the value of V at points uniformly spaced in the function’s domain ., The plots of 1/N vs . V are only linear for N ≥ 80 , and so we only used these points for the extrapolation which is read out as the value of the fit when 1/N = 0 ., To quantify the increasing curvature , Fig 5B shows the average absolute value of the second derivative of V across the function’s domain ., The coupling matrix J of both the pairwise and the semiparametric pairwise models describes effective interactions between neurons , and so it is interesting to ask how the couplings predicted by these two models are related ., While Fig 5C shows a strong dependency between the couplings in a network of N = 160 neurons , the dependency is not deterministic and , moreover , negative couplings tend to be amplified in the semiparametric pairwise model as compared to the pairwise model ., Similarly to the semiparametric independent model , there is no guarantee that the semiparametric pairwise model will reproduce observed pairwise correlations among neurons exactly , even though pairwise model has this guarantee by virtue of being a maximum entropy model ., Fig 5D shows that despite the lack of such a guarantee , semiparametric pairwise model predicts a large majority of the correlations accurately , with the possible exceptions of several very strongly correlated pairs ., This is simply because the semiparametric paiwise model is very accurate–the inset of Fig 5D shows that it can also reproduce third moments of the responses ., A K-pairwise model also has this capability but , as shown in Ref 11 , a pairwise model systematically mispredicts higher than second moments ., Suppose we use the semiparametric pairwise model to analyze a very large population which is not globally coupled and can be divided into independent subpopulations ., The only way the model in Eq ( 9 ) can be factorized into a product of probability distributions over the subpopulations is if the function V is linear ., Therefore , the prior knowledge that the population is not globally coupled immediately implies the shape of the nonlinearity ., Similarly , a prior knowledge that the population is critical also carries a lot of information about the shape of the nonlinearity ., We show in Methods that if the parameters α are known , then the optimal nonlinearity in Eq ( 3 ) can be explicitly written as, V ( E ) = log ρ ¯ ( E ; α ) - log p ¯ ^ ( E ; α ) , ( 12 ), where ρ ¯ ( E ; α ) is the density of states which counts the number of patterns s whose energy is within some narrow range E , E + Δ ., The density of states is a central quantity in statistical physics that can be estimated also for neural activity patterns either directly from data or from inferred models 19 ., Similarly , p ¯ ^ ( E ; α ) is the empirical probability density of the energy E ( s; α ) smoothed over the same scale Δ ., Eq ( 12 ) follows from the relation p ¯ ^ ( E ; α ) ∝ ρ ¯ ( E ; α ) exp ( - V ( E ) ) , i . e . the probability of some energy level is just the number of states with this energy times the probability of each of these states ( see Methods ) ., We would like to establish a prior expectation on what the large N limit of the nonlinearites in Fig 5A is ., Adapting the same normalization as in the figure , we denote ϵ ( s; α ) = E ( s; α ) /N ., Changing variables and rewriting Eq ( 12 ) in terms of the empirical probability density of the normalized energy p ϵ ¯ ^ ( ϵ ) = N p ¯ ^ ( ϵ N ; α ) yields, V ( ϵ N ) = log ρ ¯ ( ϵ N ; α ) - log p ϵ ¯ ^ ( ϵ ) + log N ., ( 13 ), For a system where si can take on two states , the total number of possible activity patterns is 2N , and so we expect the log of the density of states to be proportional to N . If the system is critical , then by virtue of Eq ( 4 ) σ ( log pN ( s ) ) is proportional to N , and similarly we also expect σ ( E ( s; α ) ) ∝ N . This means that σ ( ϵ ( s; α ) ) = σ ( E ( s; α ) ) /N converges to some finite , nonzero number , and therefore log p ϵ ¯ ^ ( ϵ ) also stays finite no matter how large the population is ., Taken together , for large critical populations , the first term on the right hand side of Eq ( 13 ) is the only one which scales linearly with the population size , and hence it dominates the other terms:, V ( E ) ≈ log ρ ¯ ( E ; α ) ., ( 14 ), One of our important results is thus that for large critical populations , the nonlinearity should converge to the density of states of the inferred energy model ., In other words , for critical systems as defined in Eq ( 4 ) , there is a precise matching relation between the nonlinearity V ( E ) and the energy function E ( s; α ) ; in theory this is exact as N → ∞ , but may hold approximately already at finite N . To verify that this is the case for our neural population that has previously been reported to be critical , we compare in Fig 6A the nonlinearity inferred with the semiparametric pairwise model ( Fig 5A ) to the density of states estimated using a Wang and Landau Monte Carlo algorithm 40 for a sequence of subpopulations of increasing size ., As the population size increases , the nonlinearity indeed approaches the regime in which our prediction in Eq ( 14 ) holds ., This convergence is further quantified in Fig 6B which shows the average squared distance between the density of states and the nonlinearity ., The average is taken over the range of observed energies ., The nonlinearities are only specified up to an additive constant which we chose so as to minimize the squared distance between the density of states and the nonlinearity ., The link between global coupling and criticality is related to recent theoretical suggestions 28 , 29 , where global coupling between the neurons in the population emerges as a result of shared latent ( fluctuating ) variables that simultaneously act on extensive subsets of neurons ., In particular , Ref 28 theoretically analyzed models with a multivariate continuous latent variable h distributed according to some probability density q ( h ) , whose influence on the population is described by the conditional probability distribution, p N ( s | h ) = e - ∑ j h j O j ( N ) ( s ) Z N ( h ) , ( 15 ), where ZN ( h ) is a normalization constant , and O j ( N ) ( s ) are global quantities which sum over the whole population ., The authors showed that under mild conditions on the probability density q ( h ) of h , and the scaling of O j ( N ) ( s ) with N , the sequence of models, p N ( s ) = ∫ q ( h ) p N ( s | h ) d h ( 16 ), is critical in the sense of Eq ( 4 ) ., If the latent variable is one-dimensional , i . e . h = h , then the models in Eq ( 16 ) have exactly the form of models in Eq ( 3 ) with E ( s; α ) = O ( s ) , i . e . given a probability density q ( h ) of the latent variable , we can always find a nonlinearity V ( E ) such that, 1 Z ( α ) e - V ( E ( s ; α ) ) = ∫ 0 ∞ q ( h ) e - h E ( s ; α ) Z ( h ; α ) d h ., ( 17 ), The reverse problem of finding a latent variable for a given function V ( E ) such that this equation is satisfied does not always have a solution ., The condition for this mapping to exist is that the function exp ( −V ( E ) ) is totally monotone 41 , which , among other things , requires that it is convex ., While our models allow for more general nonlinearites , we showed in Fig 5A that the inferred functions V ( E ) are concave and so we expect this mapping to be at least approximately possible ( see below ) ., The mapping in Eq ( 17 ) is based on a Laplace transformation , a technique commonly used for example in the study of differential equations ., Laplace transformations are also often used in statistical physics where they relate the partition function of a system to its density of states ., While the mathematics of Laplace transformations yields conditions on the function V ( E ) so that it is possible to map it to a latent variable ( i . e . , exp ( −V ( E ) ) must be totally monotone ) , analytically constructing this mapping is possible only in very special cases ., We can gain a limited amount of intuition for this mapping by considering the case when the latent variable h is a narrow gaussian with mean h0 and variance σ2 ., For small σ2 , one can show that, V ( E ) ≈ h 0 E - σ 2 ( E - E 0 ) 2 , ( 18 ), where E0 is the average energy if σ2 = 0 , and the approximation holds only in a small neighborhood of E0 ( |E − E0| ≪ σ ) ., This approximation shows that the curvature of V ( E ) is proportional to the size of the fluctuations of the latent variable which , in turn , is expected to correlate with the amount of global coupling among neurons ., This relationship to global coupling can be understood from the right hand side of Eq ( 17 ) ., When the energy function is , for example , a weighted sum of individual neurons as in the semiparametric independent model of Eq ( 8 ) , then we can think of Eq ( 17 ) as a latent variable h ( perhaps reflecting the stimulus ) coupled to every neuron , and hence inducing a coupling between the whole population ., A non-neuroscience example is that of a scene with s representing the luminance in each pixel , and the latent h representing the lighting conditions which influence all the pixels simultaneously ., We used the right hand side of Eq ( 17 ) ( see Methods ) to infer the shapes of the probability densities of the latent variables which correspond to the nonlinearities in the semiparametric pairwise models learned from data ., These probability densities are shown in Fig 6C ., A notable difference to the formulation in Eq ( 16 ) is that the inferred latent variables scale with the population size; in particular , the inset to Fig 6C shows that the entropy of the inferred latent variable increases with the population size ., Entropy is a more appropriate measure of the “broadness” of a probability density than standard deviation when the density is multimodal ., Taken together with the results in Fig 4A , this suggests that global coupling is especially important for larger populations ., However , it is also possi
Advances in multi-unit recordings pave the way for statistical modeling of activity patterns in large neural populations ., Recent studies have shown that the summed activity of all neurons strongly shapes the population response ., A separate recent finding has been that neural populations also exhibit criticality , an anomalously large dynamic range for the probabilities of different population activity patterns ., Motivated by these two observations , we introduce a class of probabilistic models which takes into account the prior knowledge that the neural population could be globally coupled and close to critical ., These models consist of an energy function which parametrizes interactions between small groups of neurons , and an arbitrary positive , strictly increasing , and twice differentiable function which maps the energy of a population pattern to its probability ., We show that:, 1 ) augmenting a pairwise Ising model with a nonlinearity yields an accurate description of the activity of retinal ganglion cells which outperforms previous models based on the summed activity of neurons;, 2 ) prior knowledge that the population is critical translates to prior expectations about the shape of the nonlinearity;, 3 ) the nonlinearity admits an interpretation in terms of a continuous latent variable globally coupling the system whose distribution we can infer from data ., Our method is independent of the underlying system’s state space; hence , it can be applied to other systems such as natural scenes or amino acid sequences of proteins which are also known to exhibit criticality .
Populations of sensory neurons represent information about the outside environment in a collective fashion ., A salient property of this distributed neural code is criticality ., Yet most models used to date to analyze recordings from large neural populations do not take this observation explicitly into account ., Here we aim to bridge this gap by designing probabilistic models whose structure reflects the expectation that the population is close to critical ., We show that such principled approach improves previously considered models , and we demonstrate a connection between our models and the presence of continuous latent variables which is a recently proposed mechanism underlying criticality in many natural systems .
linguistics, social sciences, random variables, neuroscience, covariance, probability distribution, mathematics, statistics (mathematics), thermodynamics, entropy, animal cells, probability density, probability theory, physics, statistical models, cellular neuroscience, cell biology, neurons, biology and life sciences, cellular types, physical sciences, computational linguistics