Pairwise Ising model analysis of human cortical neuron recordings

07/11/2017
Publication GSI2017
OAI : oai:www.see.asso.fr:17410:22430
contenu protégé  Document accessible sous conditions - vous devez vous connecter ou vous enregistrer pour accéder à ou acquérir ce document.
- Accès libre pour les ayants-droit
 

Résumé

During wakefulness and deep sleep brain states, cortical neural networks show a di erent behavior, with the second characterized by transients of high network activity. To investigate their impact on neuronal behavior, we apply a pairwise Ising model analysis by inferring the maximum entropy model that reproduces single and pairwise moments of the neuron's spiking activity. In this work we rst review the inference algorithm introduced in Ferrari, Phys. Rev. E (2016) [1]. We then succeed in applying the algorithm to infer the model from a large ensemble of neurons recorded by multi-electrode array in human temporal cortex. We compare the Ising model performance in capturing the statistical properties of the network activity during wakefulness and deep sleep. For the latter, the pairwise model misses relevant transients of high network activity, suggesting that additional constraints are necessary to accurately model the data.

Pairwise Ising model analysis of human cortical neuron recordings

Collection

application/pdf Pairwise Ising model analysis of human cortical neuron recordings Trang-Anh Nghiem, Olivier Marre, Alain Destexhe, Ulisse Ferrari
Détails de l'article
contenu protégé  Document accessible sous conditions - vous devez vous connecter ou vous enregistrer pour accéder à ou acquérir ce document.
- Accès libre pour les ayants-droit

Pairwise Ising model analysis of human cortical neuron recordings
application/pdf Pairwise Ising model analysis of human cortical neuron recordings

Média

Voir la vidéo

Métriques

0
0
446.81 Ko
 application/pdf
bitcache://aeba0dba9c1084149d3f80754b42272c8a7fde52

Licence

Creative Commons Aucune (Tous droits réservés)

Sponsors

Sponsors Platine

alanturinginstitutelogo.png
logothales.jpg

Sponsors Bronze

logo_enac-bleuok.jpg
imag150x185_couleur_rvb.jpg

Sponsors scientifique

logo_smf_cmjn.gif

Sponsors

smai.png
logo_gdr-mia.png
gdr_geosto_logo.png
gdr-isis.png
logo-minesparistech.jpg
logo_x.jpeg
springer-logo.png
logo-psl.png

Organisateurs

logo_see.gif
<resource  xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
                xmlns="http://datacite.org/schema/kernel-4"
                xsi:schemaLocation="http://datacite.org/schema/kernel-4 http://schema.datacite.org/meta/kernel-4/metadata.xsd">
        <identifier identifierType="DOI">10.23723/17410/22430</identifier><creators><creator><creatorName>Trang-Anh Nghiem</creatorName></creator><creator><creatorName>Olivier Marre</creatorName></creator><creator><creatorName>Alain Destexhe</creatorName></creator><creator><creatorName>Ulisse Ferrari</creatorName></creator></creators><titles>
            <title>Pairwise Ising model analysis of human cortical neuron recordings</title></titles>
        <publisher>SEE</publisher>
        <publicationYear>2018</publicationYear>
        <resourceType resourceTypeGeneral="Text">Text</resourceType><subjects><subject>Ising model</subject><subject>maximum entropy principle</subject><subject>natural gradient</subject><subject>human temporal cortex</subject><subject>multielectrode array recording</subject><subject>brain states</subject></subjects><dates>
	    <date dateType="Created">Sat 3 Mar 2018</date>
	    <date dateType="Updated">Sat 3 Mar 2018</date>
            <date dateType="Submitted">Tue 22 May 2018</date>
	</dates>
        <alternateIdentifiers>
	    <alternateIdentifier alternateIdentifierType="bitstream">aeba0dba9c1084149d3f80754b42272c8a7fde52</alternateIdentifier>
	</alternateIdentifiers>
        <formats>
	    <format>application/pdf</format>
	</formats>
	<version>37104</version>
        <descriptions>
            <description descriptionType="Abstract">During wakefulness and deep sleep brain states, cortical neural networks show a di erent behavior, with the second characterized by transients of high network activity. To investigate their impact on neuronal behavior, we apply a pairwise Ising model analysis by inferring the maximum entropy model that reproduces single and pairwise moments of the neuron's spiking activity. In this work we rst review the inference algorithm introduced in Ferrari, Phys. Rev. E (2016) [1]. We then succeed in applying the algorithm to infer the model from a large ensemble of neurons recorded by multi-electrode array in human temporal cortex. We compare the Ising model performance in capturing the statistical properties of the network activity during wakefulness and deep sleep. For the latter, the pairwise model misses relevant transients of high network activity, suggesting that additional constraints are necessary to accurately model the data.
</description>
        </descriptions>
    </resource>
.

Pairwise Ising model analysis of human cortical neuron recordings Trang-Anh Nghiem1 , Olivier Marre2 Alain Destexhe1 , and Ulisse Ferrari23 1 Laboratory of Computational Neuroscience, Unité de Neurosciences, Information et Complexité, CNRS, Gif-Sur-Yvette, France. 2 Institut de la Vision, INSERM and UMPC, 17 rue Moreau, 75012 Paris, France. 3 ulisse.ferrari@gmail.com Abstract. During wakefulness and deep sleep brain states, cortical neu- ral networks show a different behavior, with the second characterized by transients of high network activity. To investigate their impact on neu- ronal behavior, we apply a pairwise Ising model analysis by inferring the maximum entropy model that reproduces single and pairwise moments of the neuron’s spiking activity. In this work we first review the infer- ence algorithm introduced in Ferrari, Phys. Rev. E (2016) [1]. We then succeed in applying the algorithm to infer the model from a large en- semble of neurons recorded by multi-electrode array in human temporal cortex. We compare the Ising model performance in capturing the sta- tistical properties of the network activity during wakefulness and deep sleep. For the latter, the pairwise model misses relevant transients of high network activity, suggesting that additional constraints are necessary to accurately model the data. Keywords: Ising model, maximum entropy principle, natural gradient, human temporal cortex, multielectrode array recording, brain states Advances in experimental techniques have recently enabled the recording of the activity of tens to hundreds of neurons simultaneously [2] and has spurred the interest in modeling their collective behavior [3–9]. To this purpose, the pairwise Ising model has been introduced as the maximum entropy (most generic [10]) model able to reproduce the first and second empirical moments of the recorded neurons. Moreover it has already been applied to different brain regions in different animals [3, 5, 6, 9] and shown to work efficiently [11] . The inference problem for a pairwise Ising model is a computationally chal- lenging task [12], that requires devoted algorithms [13–15]. Recently, we pro- posed a data-driven algorithm and applied it on rat retinal recordings [1]. In the present work we first review the algorithm structure and then describe our successful application to a recording in the human temporal cortex [4]. We use the inferred Ising model to test if a model that reproduces empirical pairwise covariances without assuming any other additional information, also predicts empirical higher-order statistics. We apply this strategy separately to brain states of wakefulness (Awake) and Slow-Wave Sleep (SWS). In contrast to the former, the latter is known to be characterized by transients of high activity that modulate the whole population behavior [16]. Consistently, we found that the Ising model does not account for such global oscillations of the network dynamics. We do not address Rapid-Eye Movement (REM) sleep. 1 The model and the geometry of the parameter space The pairwise Ising model is a fully connected Boltzmann machine without hidden units. Consequently it belongs to the exponential family and has probability distribution: Pη X  = exp T (X) · η − log Z[η]  , (1) where X ∈ [0, 1]N is the row vector of the N system’s free variables and Z[η] is the normalization. η ∈ RD is the column vector of model parameters, with D = N(N + 1)/2 and T (X) ∈ [0, 1]D is the vector of model sufficient statistics. For the fully-connected pairwise Ising model the latter is composed of the list of free variables X and their pairwise products: {Ta(X)}D a=1 = { {Xi}N i=1 , {XiXj}N i=1,j=i+1 } ∈ [0, 1]D . (2) A dataset Ω for the inference problem is composed by a set of τΩ i.i.d. empirical configurations X: Ω = {X(t)}τΩ t=1. We cast the inference problem as a log- likelihood maximization task, which for the model (1) takes the shape: η∗ ≡ argmax η `[η] ; `[η] ≡ TΩ · η − log Z[η] , (3) where TΩ ≡ E  T (X) | Ω  is the empirical mean of the sufficient statistics. As a consequence of the exponential family properties, the log-likelihood gradient may be written as: ∇ `[η] = TΩ − Tη , (4) where Tη = E  T (X) η  is the mean of T (X) under the model distribution (1) with parameters η. Maximizing the log-likelihood is then equivalent to imposing TΩ = Tη: the inferred model then reproduces the empirical averages. Parameter space geometry In order to characterize the geometry of the model parameter space, we define the minus log-likelihood Hessian H[η], the model Fisher matrix J[η] and the model susceptibility matrix χ[η] as: χab[η] ≡ E  TaTb η  − E  Ta η  E  Tb η  (5) Jab[η] ≡ E  ∇a log Pη X  ∇b log Pη X  η  , (6) Hab[η] ≡ −∇a∇bl  η  , (7) As a property inherited from the exponential family, for the Ising model (1): χab[η] = Jab[η] = Hab[η] . (8) This last property is the keystone of the present algorithm. Moreover, the fact that the log-likelihood Hessian can be expressed as a co- variance matrix ensures its non-negativity. Some zero Eigenvalues can be present, but they can easily be addressed by L2-regularization [14, 1]. The inference prob- lem is indeed convex and consequently the solution of (3) exists and is unique. 2 Inference algorithm The inference task (3) is an hard problem because the partition function Z[η] cannot be computed analytically. Ref. [1] suggests applying an approximated natural gradient method to numerically address the problem. After an initial- ization of the parameters to some initial value η0, the natural gradient [17, 18] iteratively updates their values with: ηn+1 = ηn − αJ−1 [ηn] · ∇ `[ηn] . (9) For sufficiently small α, the convexity of the problem and the positiveness of the Fisher matrix ensure the convergence of the dynamics to the solution η∗ . As computing J[ηn] at each n is computationally expensive, we use (8) to approximate the Fisher with an empirical estimate of the susceptibility [1]: J[η] = χ[η] ≈ χ[η∗ ] ≈ χΩ ≡ Cov  T Ω  . (10) The first approximation becomes exact upon convergence of the dynamics, ηn → η∗ . The second assumes that (i) the distribution underlying the data belongs to the family (1), and that (ii) the error in the estimate of χΩ, arising from the dataset’s finite size, is small. We compute χΩ of Eq. (10) only once, and then we run the inference algo- rithm that performs the following approximated natural gradient: ηn+1 = ηn − αχ−1 Ω · ∇ `[ηn] . (11) Stochastic dynamics.4 The dynamics (11) require estimating ∇ `[η] and thus of Tη at each iteration. This is accounted by a Metropolis Markov-Chain Monte Carlo (MC), which collects Γη, a sequence of τΓ i.i.d. samples of the distribution (1) with parameters η and therefore estimates: TMC η ≡ E  T (X) Γη  . (12) This estimate itself is a random variable with mean and covariance given by: E  TMC η {Γη}  = Tη ; Cov  TMC η {Γη}  = J[η] τΓ , (13) where E  · {Γη}  means expectation with respect to the possible realizations Γη of the configuration sequence. 4 The results of this section are grounded on the repeated use of central limit theorem. See [1] for more detail. Data: TΩ,χΩ Result: η∗ ,TMC η∗ Initialization: set τΓ = τΩ,α = 1 and η0; estimate TMC η0 and compute 0 ; while  > 1 do ηn+1 ← ηn −αχ−1 Ω · ∇ l[ηn]; estimate TMC ηn+1 and compute n+1; if n+1 < n then increase α, keeping α ≤ 1 ; else decrease α and set ηn+1 = ηn; end n ← n + 1; end Fix α < 1 and perform several iterations. Algorithm 1: Algorithm pseudocode for the Ising model inference. For η sufficiently close to η∗ , after enough iterations, this last result allows us to compute the first two moments of ∇`MC η ≡ TΩ − TMC η , using a second order expansion of the log-likelihood (3): E  ∇`MC η {Γη}  = H[η] · (η − η∗ ) ; Cov  ∇`MC η {Γη}  = J[η∗ ] τΓ . (14) In this framework, the learning dynamics becomes stochastic and ruled by the master equation: Pn+1(η0 ) = Z dη Pn(η) Wη→η0 [η] ; Wη→η0 [η] = Prob ∇`MC η = η0 − η  , (15) where Wη→η[η] is the probability of transition from η to η0 . For sufficiently large τΓ and thanks to the equalities (8), the central limit theorem ensures that the unique stationary solution of (15) is a Normal Distribution with moments: E  η P∞(η)  = η∗ ; Cov  η P∞(η)  = α (2 − α)τΓ χ−1 [η∗ ] . (16) Algorithm. Thanks to (8) one may compute the mean and covariance of the model posterior distribution (with flat prior): E  η PPost (η)  = η∗ ; Cov  η PPost (η)  = 1 τΩ χ−1 [η∗ ] (17) where τΩ is the size of the training dataset. From (14), if η ∼ PPost we have: E  ∇`MC η {Γη∼P Post }  = 0 ; Cov  ∇`MC η {Γη∼P Post }  = 2χ[η∗ ] τΓ . (18) Interestingly, by imposing: 1 τΩ = α (2 − α)τΓ (19) the moments (16) equal (17) [1]. To evaluate the inference error at each iteration we define: n =

∇`MC ηn

χΩ = r τΩ 2D ∇`MC ηn ·χ−1 Ω · ∇`MC ηn . (20) Averaging  over the posterior distribution, see (18), gives  = 1. Consequently, if ηn 6= η∗ implies n > 1 with high probability, for ηn → η∗ thanks to (19) we expect n = p τΩ/τΓ /(2 − α) [1]. As sketched in pseudocode 1, we iteratively update ηn through (11) with τΓ = τΩ and α < 1 until n < 1 is reached. 3 Analysis of cortical recording 0.07 0.03 0.01 0.05 Cov[ X | DATA ] 0.07 0.03 0.01 0.05 Cov[ X | MODEL ] Awake 0.07 0.03 0.01 0.05 Cov[ X | DATA ] 0.07 0.03 0.01 0.05 Cov[ X | MODEL ] SWS Fig. 1. Empirical pairwise covariances against their model prediction for Awake and SWS. The goodness of the match implies that the inference task was successfully com- pleted. Note the larger values in SWS than Awake As in [4, 7], we analyze ∼ 12 hours of intracranial multi-electrode array recording of neurons in the temporal cortex of a single human patient. The dataset is composed of the spike times of N = 59 neurons, including NI = 16 inhibitory neurons and NE = 43 excitatory neurons. During the recording ses- sion, the subject alternates between different brain states [4]. Here we focused on wakefulness (Awake) and Slow-Wave Sleep (SWS) periods. First, we divided each recording into τΩ short 50ms-long time bins and encoded the activity of each neuron i = 1, . . . , N in each time bin t = 1, . . . , τΩ as a binary variable Xi(t) ∈ [0, 1] depending on whether the cell i was silent (Xi(t) = 0) or emit- ted at least one spike (Xi(t) = 1) in the time window t. We thus obtain one training dataset Ω = {{Xi(t)}N i=1}τΩ t=1 per brain state of interest. To apply the Ising model we assume that this binary representation of the spiking activity is representative of the neural dynamics and that subsequent time-bins can be con- sidered as independent. We then run the inference algorithm on the two datasets separately to obtain two sets of Ising model parameters ηAwake and ηSWS. Thanks to (4), when the log-likelihood is maximized, the pairwise Ising model reproduces the covariances E  XiXj | Ω  for all pairs i 6= j. To validate the inference method, in Fig. 1 we compare the empirical and model-predicted pair- wise covariances and found that the first were always accurately predicted by the second in both Awake and SWS periods. 0 2 4 6 8 10 12 14 K 10-5 10-4 10-3 10-2 10-1 100 Awake Independent p(K) data p(K) model 0 2 4 6 0 2 4 6 8 10 12 14 K 10-5 10-4 10-3 10-2 10-1 100 SWS Independent p(K) data p(K) model 0 2 4 6 Fig. 2. Empirical and predicted distributions of the whole population activity K = P i Xi. For both Awake and SWS periods the pairwise Ising model outperforms the independent model (see text). However, Ising is more efficient at capturing the popu- lation statistics during Awake than SWS, expecially for medium and large K values. This is consistent with the presence of transients of high activity during SWS. This shows that the inference method is successful. Now we will test if this model can describe well the statistics of the population activity. In particular, synchronous events involving many neurons may not be well accounted by the pairwise nature of the Ising model interactions. To test this, as introduced in Ref. [6], we quantify the empirical probability of having K neurons active in the same time window: K = P i Xi. In Fig. 2 we compare empirical and model prediction for P(K) alongside with the prediction from an independent neurons model, the maximum entropy model that as sufficient statistics has only the single variables and not the pairwise: {Ta(X)}N a=1 = {Xi}N i=1. We observed that the Ising model always outperforms the independent model in predicting P K  . Fig. 2 shows that the model performance are slightly better for Awake than SWS states. This is confirmed by a Kullback-Leibler divergence estimate: DKL PData Awake(K) PIsing Awake(K)  = 0.005; DKL PData SWS (K) PIsing SWS (K)  = 0.030 . This effect can be ascribed to the presence of high activity transients, known to modulate neurons activity during SWS [16] and responsible for the larger covariances, see Fig. 1 and the heavier tail of P(K), Fig. 2. These transients are know to be related to an unbalance between the contributions of excitatory and inhibitory cells to the total population activity [7]. To investigate the impact of these transients, in Fig. 3 we compare P(K) for the two populations with the cor- responding Ising model predictions. For the Awake state, the two contributions are very similar, probably in consequence of the excitatory/inhibitory balance [7]. Moreover the model is able to reproduce both behaviors. For SWS periods, instead, the two populations are less balanced [7], with the inhibitory (blue line) showing a much heavier tail. Moreover, the model partially fails in reproducing this behavior, notably strongly overestimating large K probabilities. 0 2 4 6 8 10 K 10-5 10-4 10-3 10-2 10-1 100 Awake p(K) E data p(K) E model p(K) I data p(K) I model 0 2 4 6 8 10 K 10-5 10-4 10-3 10-2 10-1 100 SWS p(K) E data p(K) E model p(K) I data p(K) I model Fig. 3. Empirical and predicted distributions of excitatory (red) and inhibitory (blue) population activity. During SWS, the pairwise Ising model fails at reproducing high activity transients, especially for inhibitory cells. Conclusions: (i) The pairwise Ising model offers a good description of the neural network activity observed during wakefulness. (ii) By contrast, tak- ing into account pairwise correlations is not sufficient to describe the statistics of the ensemble activity during SWS, where (iii) alternating periods of high and low network activity introduce high order correlations among neurons, especially for inhibitory cells [16]. (iv) This suggests that neural interactions during wake- fulness are more local and short-range, whereas (v) these in SWS are partially modulated by internally-generated activity, synchronizing neural activity across long distances [4, 16, 19]. Acknowledgments We thank B. Telenczuk, G. Tkacik and M. Di Volo for useful discussion. Research funded by European Community (Human Brain Project, H2020-720270), ANR TRAJECTORY, ANR OPTIMA, French State program Investissements dAvenir managed by the Agence Nationale de la Recherche [LIFESENSES: ANR-10-LABX-65] and NIH grant U01NS09050. References 1. U. Ferrari. Learning maximum entropy models from finite-size data sets: A fast data-driven algorithm allows sampling from the posterior distribution. Phys. Rev. E, 94:023301, 2016. 2. I.H. Stevenson and K.P. Kording. How advances in neural recording affect data analysis. Nature neuroscience, 14(2):139–142, 2011. 3. E. Schneidman, M. Berry, R. Segev, and W. Bialek. Weak pairwise correlations imply strongly correlated network states in a population. Nature, 440:1007, 2006. 4. A. Peyrache, N. Dehghani, E.N. Eskandar, J.R. Madsen, W.S. Anderson, J.A. Donoghue, L.R. Hochberg, E. Halgren, S.S. Cash, and A. Destexhe. Spatiotemporal dynamics of neocortical excitation and inhibition during human sleep. Proceedings of the National Academy of Sciences, 109(5):1731–1736, 2012. 5. L. S. Hamilton, J. Sohl-Dickstein, A. G. Huth, V. M. Carels, K. Deisseroth, and S. Bao. Optogenetic Activation of an Inhibitory Network Enhances Feedforward Functional Connectivity in Auditory Cortex. Neuron, 80:1066–76, 2013. 6. G. Tkacik, O. Marre, D. Amodei, E. Schneidman, W Bialek, and Berry M.J. Search- ing for collective behaviour in a network of real neurons. PloS Comput. Biol., 10(1):e1003408, 2014. 7. N. Dehghani, A. Peyrache, B. Telenczuk, M. Le Van Quyen, E. Halgren, S.S. Cash, N.G. Hatsopoulos, and A. Destexhe. Dynamic balance of excitation and inhibition in human and monkey neocortex. Scientific Reports, 6(23176), 2016. 8. C. Gardella, O. Marre, and T. Mora. A tractable method for describing complex couplings between neurons and population rate. eneuro, 3(4):0160, 2016. 9. G. Tavoni, U. Ferrari, S. Cocco, F.P. Battaglia, and R. Monasson. Functional cou- pling networks inferred from prefrontal cortex activity show experience-related ef- fective plasticity. Network Neuroscience, MIT press, 2017. 10. E.T. Jaynes. On The Rationale of Maximum-Entropy Method. Proc. IEEE, 70:939, 1982. 11. U. Ferrari, T. Obuchi, and T. Mora. Random versus maximum entropy models of neural population activity. Phys. Rev. E, 95:042321, 2017. 12. D. H. Ackley, G. E. Hinton, and T. J. Sejnowski. A learning algorithm for boltz- mann machines. Cognitive Science, 9:147–169, 1985. 13. T. Broderick, M. Dudik, G. Tkacik, R.E. Schapire, and W. Bialek. Faster solutions to the inverse pairwise Ising problem. Arxiv:0712.2437, 2007. 14. S. Cocco and R. Monasson. Adaptive cluster expansion for inferring Boltzmann machines with noisy data. Phys. Rev. Lett., 106:090601, 2011. 15. J. Sohl-Dickstein, P. B. Battaglino, and M. R. DeWeese. New Method for Param- eter Estimation in Probabilistic Models: Minimum Probability Flow. Phys. Rev. Lett., 107:220601, 2011. 16. A. Renart, J. De La Rocha, P. Bartho, L. Hollender, N. Parga, A. Reyes, and K.D. Harris. The asynchronous state in cortical circuits. Science, 327(5965):587– 590, 2010. 17. S. Amari. Natural Gradient Works Efficiently in Learning, Neural Computation. Neural Comput., 10:251–276, 1998. 18. S. Amari and S.C. Douglas. Why natural gradient? Proc. IEEE, 2:1213–16, 1998. 19. M. Le Van Quyen,L.E. Muller, B. Telenczuk, E. Halgren, S. Cash, N.G. Hatsopou- los, N. Dehghani, and A. Destexhe High-frequency oscillations in human and mon- key neocortex during the wake–sleep cycle. Proceedings of the National Academy of Sciences, 113, 33, 9363–68, 2016