(C) PLOS One [1]. This unaltered content originally appeared in journals.plosone.org.
Licensed under Creative Commons Attribution (CC BY) license.
url:https://journals.plos.org/plosone/s/licenses-and-copyright

------------



Parametric Copula-GP model for analyzing multidimensional neuronal and behavioral relationships

['Nina Kudryashova', 'School Of Informatics', 'University Of Edinburgh', 'Edinburgh', 'United Kingdom', 'Theoklitos Amvrosiadis', 'Centre For Discovery Brain Sciences', 'Nathalie Dupuy', 'Nathalie Rochefort', 'Simons Initiative For The Developing Brain']

Date: 2022-02 

Each panel shows a scatter plot of samples drawn from a parametric copula family (named in the title of each plot) with a fixed parameter θ (shown in the bottom right corner). In total, we used 10 different copula elements—Gaussian + Frank + 4 × Clayton + 4 × Gumbel—to construct copula mixtures. The rightmost figure shows a mixture of two copulas from different copula families taken with equal mixing weights (0.5). Blue points here correspond to the samples drawn from a Clayton copula, orange points—to a 90°-rotated Gumbel copula. Note, that a mixture of copulas is a copula itself.

We are using Gaussian Processes (GP) to model the conditional dependencies of copula parameters (i.e. c(.|x)), following the approach by Hernández-Lobato et al. [ 31 ] which was originally developed for analysis of financial time-series. GPs are ideally suited for copula parametrization, as they make the dependencies explicit and provide an estimate of the uncertainty in model parameters, which we utilize in our model selection process (see Sec. 6 in Methods ). Similarly to previous computational studies [ 23 , 28 , 30 , 36 , 37 ], we use vine copula constructions to scale the method to higher dimensions (Sec. 7). Here, we extend the Copula-GP framework with more flexible copula mixture models (see Fig 1 and Methods ), suitable for inter-neuronal dependencies, and utilise the repeated trial structure in order to estimate conditional marginal distributions p i (y i |x) (e.g. distributions of single neuron responses at every moment x in the trial). We develop model selection algorithms, based on the fully-Bayesian Watanabe–Akaike information criterion (WAIC), that construct the best suited mixture models for a given data set. In the following section we show the utility of our extended Copula-GP model for analysis of neuronal and behavioral data.

The model in (1) has been previously applied to neuronal data [ 18 , 23 , 35 ]. In this paper, we extend that model by conditioning it on a continuous variable x: (2) where x, depending on the design of the experimental task, can represent any continuous task-related variable such as time, phase, coordinate in space, velocity, direction or orientation. In the further sections we will show that if this variable x also uniquely determines the stimulus, then the components of the model (2) can also be interpreted in terms of ‘noise correlations’ (i.e. joint trial-to-trial variability of responses to a given stimulus).

Many questions in neuroscience require accurate statistical models of neuronal responses to a certain stimulus. For instance, such models can be used for Bayesian decoding of the stimulus [ 20 ] or for information-theoretic analysis [ 33 ] in order to assess the coding precision of a neuronal population. Here we describe a general framework for constructing such statistical models, called Copula-GP; the technical details related to model implementation and fitting can be found in Methods.

2 Copula-GP separates noise correlations from stimulus statistics

We first examine a synthetic ‘toy’ model of calcium imaging data in order to illustrate and interpret the components of our Copula-GP model. We use a Generalized linear model (GLM) [38] and a model of the kinetics of the calcium indicator [39] in order to generate synthetic calcium imaging data of the activity of two neurons (see Methods). GLMs are widely used to model population spike trains based on spike history filters and coupling filters. While these models can be readily fit to small neuronal populations, they have parameter identifiability issues that lead to degraded performance in predicting neuronal activity patterns [40]. GLMs are also not directly applicable to calcium imaging data and require some additional spike inference algorithm (deconvolution, e.g. Deneux et al. 2016 [39]) to obtain spikes. So, instead of fitting GLMs, we propose to explicitly model the dynamic neuronal activity patterns (including the calcium dynamics) with Copula-GP. Our approach is complementary to GLMs: while GLMs consider time-lagged interactions via temporal spike history filters and cannot model instantaneous coactivation [41], our model focuses exclusively on the instantaneous correlations between calcium signals.

Calcium imaging in neurons is used as a proxy for monitoring their spiking activity since action potentials are tightly coupled to large increases in intracellular free calcium concentration [42]. Two-photon imaging is used to monitor fluorescence changes of calcium indicators in large population of neurons, with single cell resolution. When dense populations of neurons are imaged, the signals imaged in a given focal plane can be contaminated by signals from a volume around this plane. Different methods for correcting such out-of-focus contamination have been developed [43]. We assumed that each decontaminated recording y i (t) of calcium dynamics for each neuron is a function of its spiking activity n i (t), which is independent of the activity of other neurons (i.e. y i (n 1 …n N , t) = y i (n i , t), given time). This assumption is also implicitly made in deconvolution algorithms, which incorporate inductive biases for calcium dynamics and are also applied to individual neurons independently. We will consider the decontaminated recordings y i (t) to be a surrogate of neural spike trains n i (t), which can be approximated by blurring (convolving) the spikes with an exponential kernel (Deneux et al. 2016 [39]). The time constant in this kernel corresponds to calcium decay and also sets the time scale of interactions between spiking neurons, which would be captured by an instantaneous correlation of their calcium traces. Keeping these assumptions in mind, we aim to build instantaneous dependence models of calcium recordings for statistical and information-theoretic analysis.

We first consider a pair of neurons, which have the same tuning (receive the same input x(t) in Fig 2A), but are not coupled. Since the spiking activity of these neurons in response to a given stimulus is completely independent, we expect that their calcium traces, which are a temporally blurred versions of spike trains, to also be independent given the stimulus. A few simulated calcium transients (fluorescence across time) for one of these identical neurons are shown in Fig 2B.

PPT PowerPoint slide

PNG larger image

TIFF original image Download: Fig 2. Copula-GP finds that uncoupled neurons are independent given the stimulus. A GLM model of two identical uncoupled neurons that receive the same time-dependent input x(t); B simulated calcium transients (fluorescence across time) showing dynamic responses to the stimulus x(t) for one of the neurons; C calcium transients of two neurons (y 1 (t), y 2 (t)) projected onto a unit cube by the probability integral transform based on unconditional marginals; colored points show transformed samples (u 1 , u 2 ) corresponding to times t (color-coded). The clusters of similarly colored points (e.g. green) illustrate that the copula c(u) depends on time t; the particular shape and the location of the clusters depends on the function x(t); only 10% of data-points are shown (selected randomly). D same as C, but based on conditional marginals F i (y i |t). The resulting empirical copula describes ‘noise correlations’ between two neurons. The colored data-points ( , ) are uniformly distributed on the unit square, which suggests that there is no noise correlation between these neurons, the copula c(ut) is independent of time t, and the neurons are independent given the time-dependent stimulus. https://doi.org/10.1371/journal.pcbi.1009799.g002

By design, the activity of these neurons is independent given the stimulus (Fig 2A), but they do respond similarly to the stimulus. As a result of a naive analysis of their joint statistics, such as measuring correlation between recorded activities over time, their activity would appear dependent. The structure of this dependence can be visualized with a static copula model, as in (1).

To fit a copula-based model to a given dataset, one typically starts with the marginals (i.e. single variable distributions). Here, we model the single neuron marginals with their empirical distributions, estimated directly from the data in the form of empirical cumulative distribution functions (eCDFs). We then use these eCDFs to project the simulated neuronal recordings onto a unit cube using the probability integral transform: F(y i )→u i ∼ U [0,1] , such that each u i becomes uniformly distributed. Note, that both the neural responses y i and their transformed versions u i still depend on both the stimulus x(t) and time t. We can confirm that by visualizing the empirical copula c(u 1 …u N ) = c(F 1 (y 1 )…F N (y N )) (Fig 2C). The colors of the data points here indicate time t. Similarly colored clusters of points in copula space (e.g. green in Fig 2C) indicate that both the dependence and the marginals in such an unconditional model depend on time t and on the presented stimulus x(t). Hence, such dependence characterizes the joint statistics of neural responses to stimulus x(t), but not the interaction between neurons.

Spurious correlations between neuronal recordings are a well known problem when modeling responses to complex stimuli (e.g. time-dependent sequences of stimuli x(t)), which can be solved by distinguishing between noise and stimulus correlations [40, 44]. In our Copula-GP model, this separation is achieved by using conditional marginals and copulas:

Here, we only assume that the stimulus x(t) is fully determined by time t. In general, our framework (2) is applicable to any complex stimuli determined by any continuous variable (e.g. time, phase, coordinate in space, velocity, direction, orientation, etc.). The conditional marginals p i (y i |t) in this model account for the within-trial variability (e.g. dynamics of responses), while the copula model accounts for the trial-to-trial variability in neural responses.

The conditional marginals p i (y i |t) are estimated using the non-parametric fastKDE [45] algorithm. This approach is suitable for relatively large datasets, which have enough data points for direct estimation of the eCDF (see Sec D in S3 Text). If the number of datapoints is insufficient, one might consider using parametric marginal models instead. Note that the fastKDE also assumes smooth changes in the marginal distributions with respect to the conditioning variable. The corresponding conditional marginal CDFs F i (y i |t) are then used to map the data onto a unit hypercube using: , such that is uniformly distributed for any t. The resulting empirical copula correctly demonstrates the absence of dependence between y 1 (t) and y 2 (t), as the points are uniformly distributed on a unit square. The density of points in Fig 2D illustrates that for every given value of the transformed neural response , the conditional distribution of the other neural responses is the same ( ). This means that the variables and are independent , and their copula is the Independence copula. This result demonstrates that once we made the distinction between noise correlations and stimulus correlations, we could correctly identify that these two neurons are uncoupled (Fig 2A).

The noise correlations were previously linked to anatomical and functional relationships between cortical neurons [46, 47]. However, these relationships should be interpreted with caution, since correlation between neurons does not necessarily imply the mechanistic coupling between them [48]. Despite lack of mechanistic interpretability, accurate modelling of noise correlations is useful for understanding neural code from the information-theoretic perspective. Generally, noise correlations themselves can depend on the stimulus [49], and taking their stimulus dependence into account can improve the decoding accuracy [44, 50, 51]. In order to model such tuning of the noise correlations, not only the marginals but also the corresponding copula c(…|t) itself must be conditioned on t, as in (2).

We next consider two coupled neurons: one excitatory and one inhibitory (Fig 3A.i). They again receive the same input x(t) as in Fig 2A, in all trials. We added two time-dependent filters h 12 and h 21 , that couple the spike train history of each neuron to the other (Fig 3A.ii) [38]. The synthetic calcium traces (Fig 3B and 3C) demonstrate some non-trivial activity in both neurons after the stimulus presentation window, where all of the recurrent circuit dynamics unfolds. We expect that these dynamics will be reflected in some non-independent trial-to-trial co-variability of neural responses. Since one neuron inhibits the other, we also expect that their responses will be overall negatively correlated. While the unconditional copula analysis again demonstrates strong stimulus correlations (Fig 3D), the conditional copula (conditioned on the stimulus) reveals some structure in noise correlations (Fig 3E, see more points are concentrated in the upper-left corner).

PPT PowerPoint slide

PNG larger image

TIFF original image Download: Fig 3. Copula-GP describes the noise correlations between dynamically coupled neurons with a Clayton copula. A i. GLM model of two coupled neurons (excitatory and inhibitory) that receive the same time-dependent input x(t); ii. the spike history coupling filters h 12 and h 21 ; B-C simulated calcium transients (fluorescence across time) showing dynamic responses to the stimulus x(t) for excitatory and inhibitory neurons, respectively; D calcium transients of two neurons (y 1 (t), y 2 (t)) projected onto a unit cube by the probability integral transform based on unconditional marginals; colored points show transformed samples (u 1 , u 2 ) corresponding to times t (color-coded). The clusters of similarly colored points (e.g. green) illustrate that the copula c(u) depends on time t; the particular shape and the location of the clusters depends on the function x(t). E same as D, but based on conditional marginals F i (y i |t). The resulting copula describes ‘noise correlations’ between two neurons. The colored data-points ( , ) are not uniformly distributed on the unit square, which suggests that the noise correlation between these neurons and the copula c(ut) itself depends on time t. F Clayton copula parameter (θ) that characterizes the strength of the non-linear noise correlation between neurons (see Methods for details); G probability density plots illustrating the stimulus-dependent shape of noise correlations. The empirical dependence estimated from data samples is shown with black outlines, while the predictions of the Clayton copula model are shown in shades of blue. The proportion of the variance explained is indicated in the upper-right corner for each time interval. Orange circles indicate the heavy tail of the distribution, which can be best seen in the range t ∈ [0.6, 1.0] where the variables are stronger correlated. https://doi.org/10.1371/journal.pcbi.1009799.g003

In order to analyse the structure of the noise correlations in this example, we apply our Copula-GP model. For simplicity, here we use a single 90°-rotated Clayton copula, parameterized by a Gaussian Process (GP); for details regarding the inference scheme, see Methods. The inferred GP parameter reconstructs the stimulus-dependent changes in noise correlations (Fig 3F), which are most pronounced after the stimulation window. The corresponding Clayton copula model can accurately describe the shape of the conditional dependence, which we quantify with the proportion of variance explained in Fig 3G (see Sec. 10 in Methods).

This noise correlation model in Fig 3G can be interpreted as follows. The red point in the middle of each plot in Fig 3G corresponds to the median response for both neurons: 50% of recorded responses from a given neuron were higher, while 50% were lower. More generally, the values of the transformed neural responses ut correspond to the percentile scores of each response y in a marginal distribution of neural responses. The shades in Fig 3G correspond to the probability density, i.e. how likely it is to jointly observe a pair of responses and . Typically, most of the density is concentrated along one of the diagonals of the unit square: around if neurons are positively correlated, and around if neurons are negatively correlated. In this case, the density is mostly concentrated around the negative diagonal, suggesting that the responses of these two neurons are generally negatively correlated (as expected, given that one of the neurons inhibits the activity of the other). However, this dependence also has a heavy tail: a high degree of association between some of the extreme values of the variables. We can observe this association as a probability mass concentrated in the corner of the unit square (see orange circles in Fig 3G). It indicates that there is a high probability to observe a strongly activated inhibitory neuron together with the inactivated excitatory neuron. Such asymmetry in the joint trial-to-trial variability of two neurons suggests that their dependence is non-linear: their dependence is not well characterized by a single linear correlation coefficient as there is additional structure in the dependence beyond linear association. In our GLM model, this non-linearity stems from the different timescales of spike history filters: fast inhibitory and slow excitatory neurons.

The confidence intervals (95% CI, two standard deviations) shown in Fig 3F correspond to the uncertainty in model parameters, captured by the Gaussian process. The uncertainty in model parameters is not essential for this example, but necessary when modeling high-dimensional datasets with relatively low sample numbers, which is often the case in neuronal data with hundreds to thousands of recorded neurons and only few hundreds of trials. These uncertainties can be propagated through the model and produce the uncertainty measures for the decoded stimulus (when using Copula-GP for Bayesian decoding) or in estimated mutual information.

As we demonstrated in these synthetic examples, our Copula-GP model can separate the statistics corresponding to the network from the single unit statistics, or, in other words, distinguish (stimulus-dependent) noise correlations (i.e. shared neuronal trial-to-trial variability) from stimulus correlations. Moreover, it provides a visualization of the shape of noise correlations (Fig 3G), and reveals that the noise correlations are tuned to a particular time point (t ≈ 0.7) in a trial (Fig 3F). The main advantage of using copula models [18] is their expressive power in modeling the shape of the relationship between variables, which, in this case, is asymmetric due to the dynamic connections in the neuronal circuit being asymmetric (Fig 3). It was also shown previously that shared excitatory or inhibitory inputs can result in similar non-linear dependencies between neurons [24]. While copula models of the noise correlations do not provide any mechanistic interpretation of these dependencies (e.g. whether neurons are recurrently coupled or share an input), their expressive power allows one to build statistical models for information-theoretic analysis or Bayesian decoding. An additional explicit fully-Bayesian parameterization provided by a GP (Fig 3F) extends the copula model to time-dependent, location-dependent or otherwise continuously parameterized behavioral tasks.

[END]

[1] Url: https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1009799

(C) Plos One. "Accelerating the publication of peer-reviewed science."
Licensed under Creative Commons Attribution (CC BY 4.0)
URL: https://creativecommons.org/licenses/by/4.0/


via Magical.Fish Gopher News Feeds:
gopher://magical.fish/1/feeds/news/plosone/