Generative embeddings of brain collective dynamics using variational autoencoders
Yonatan Sanz Perl, Hernán Boccacio, Ignacio Pérez-Ipiña, Federico Zamberlán, Helmut Laufs, Morten Kringelbach, Gustavo Deco, Enzo Tagliazucchi
GGenerative embeddings of brain collective dynamics using variational autoencoders
Yonatan Sanz Perl,
1, 2
Hernn Boccacio, Ignacio Prez-Ipia, Federico Zamberln, Helmut Laufs, Morten Kringelbach, Gustavo Deco, and Enzo Tagliazucchi Universidad de San Andrs Physics Department (University of Buenos Aires) and Buenos Aires Physics Institute (CONICET) Department of Neurology, Christian-Albrechts-University Kiel Department of Psychiatry, University of Oxford Center for Brain and Cognition, Computational Neuroscience Group, Universitat Pompeu Fabra (Dated: July 6, 2020)We consider the problem of encoding pairwise correlations between coupled dynamical systems ina low-dimensional latent space based on few distinct observations. We used variational autoencoders(VAE) to embed temporal correlations between coupled nonlinear oscillators that model brain statesin the wake-sleep cycle into a two-dimensional manifold. Training a VAE with samples generatedusing two different parameter combinations resulted in an embedding that represented the wholerepertoire of collective dynamics, as well as the topology of the underlying connectivity network.We first followed this approach to infer the trajectory of brain states measured from wakefulness todeep sleep from the two endpoints of this trajectory; next, we showed that the same architecturewas capable of representing the pairwise correlations of generic Landau-Stuart oscillators coupledby complex network topology
Several biological systems can be understood in termsof simple dynamical rules coupled by heterogeneous con-nectivity patterns. Perhaps the most paradigmatic caseis the human brain, where complex collective behaviouremerges from the nonlinear dynamics of ≈ neuronsinteracting at ≈ synaptic connections [1]. In spite ofthis complexity at the microscopic scale, the brain spon-taneously self-organises into a reduced number of dis-crete states, such as those in the wake-sleep cycle, whichsuggests that a low-dimensional manifold is sufficient toencode its large-scale dynamics [2].The mechanisms underlying the emergence of differ-ent brain states can be probed using whole-brain modelsbased on conceptually simple local dynamical rules cou-pled according to empirical measurements of anatomicalconnectivity [3–5], for instance, by coupling nonlinear os-cillators with the long-range white matter tracts inferredfrom diffusion tensor imaging (DTI) [6]. After parameteroptimisation to reproduce neuroimaging data acquiredduring different brain states, the models can be used toexplore the interplay between local dynamics, long-rangestructural coupling, and the formation of large-scale ac-tivity patterns [7–9], and as methods for data augmen-tation to be combined with machine learning techniquesfor the purpose of brain state classification [10, 11].While whole-brain models can reproduce the func-tional connectivity of brain states such as those seen inthe progression from wakefulness to deep sleep [7, 12],it is unclear whether coupled dynamical systems canalso capture relationships between these states, encod-ing them into a low dimensional manifold that preservesthe ordering within progressions of brain states. Moregenerally, we consider a system of coupled units whosedynamics have been optimised to reproduce the second-order statistics (i.e. pairwise correlations) of a real-world system, and ask whether different discrete states of suchsystem can be efficiently represented by latent variablesthat are capable of reproducing the whole repertoire ofstates from a reduced number of representative examples.In the particular case of collective brain dynamics, thisis equivalent to asking whether the endpoint states of acertain progression, such as the descent from wakeful-ness into deep sleep, can be used to learn a latent repre-sentation which encodes all intermediate stages, and canbe extrapolated to produce correlations corresponding tostates beyond this progression.We used whole-brain models fitted to empirical datato generate pairwise correlation matrices for the differ-ent brain states that comprise the human wake-sleepcycle: wakefulness, N1, N2 and N3 sleep (N1 and N2are intermediate stages, while N3 is the deepest stage ofhuman sleep) [13]. Next, we trained a variational au-toencoder (VAE) with matrices corresponding to wake-fulness and N3 sleep, showing that intermediate (N1 andN2) sleep stages were embedded continuously in the la-tent space, and that the resulting two-dimensional man-ifold also extrapolated to capture known results con-cerning the structure-function relationship during uncon-sciousness [14–16]. Finally, we assessed the relationshipbetween latent space variables and the parameters ofgeneric coupled Stuart-Landau oscillators. Whole-brain model.
We start from a model constructedfrom 90 Stuart-Landau nonlinear oscillators, each repre-senting the dynamics within a macroscopic brain regionof interest [6]. The coupled dynamics are given by, a r X i v : . [ q - b i o . N C ] J u l dx j dt = ( a − x j − y j ) x j − ω j y j + G (cid:88) i J ij ( x i − x j ) + βη j (1) dy j dt = ( a − x j − y j ) y j + ω j x j + G (cid:88) i J ij ( y i − y j ) + βη j Where x j is the dynamical variable that simulates thefunctional magnetic resonance (fMRI) signal of region j and J ij represents the symmetrical coupling matrix thatweights the connectivity between regions i and j . Thismatrix is inferred from the diffusion of water moleculesin white matter from DTI recordings, and represents theempirical distribution of long-range axon bundles in thebrain. The bifurcation parameter ( a ) controls the prox-imity to oscillatory dynamics and G globally scales thecoupling between oscillators. Finally, η j is an additiveGaussian noise term that is scaled by β = 0 . a j ) into 6 parameters representing the contributionof each RSN to the local dynamics as a j = (cid:80) k ∆ k jk ,where jk equals 1 if the node j belongs to the k -th RSNand zero otherwise. We then applied a stochastic opti-misation algorithm to determine the ∆ k and G that bestreproduce the correlation matrix C ij of each state in theprogression from wakefulness to deep sleep. The C ij con-tains in its i, j entry the linear correlation between theempirical/simulated fMRI time series corresponding tonodes i and j [7]. Encoding the C ij with a VAE. We implemented aVAE to find a low-dimensional representation encodingthe progression of brain states. VAE are autoencoderstrained to map inputs to probability distributions in la-tent space, which can be regularised during the train-ing process to produce meaningful outputs after the de-coding step. The architecture of the implemented VAE(shown in Fig. 1) can be subdivided into three parts: theencoder network, the middle variational layer, and thedecoder network. The encoder consists of a deep neu-ral network with rectified linear units (ReLu) as activa-tion functions and two dense layers, which bottlenecksinto the two-dimensional variational layer, where units z and z span the latent space. The encoder networkapplies a nonlinear transformation to map the C ij intoGaussian probability distributions in latent space, andthe decoder network mirrors the encoder architecture toproduce reconstructed matrices C ∗ ij from samples of thesedistributions. FIG. 1. VAE architecture. The inputs are correlation matri-ces C ij obtained from the model (Eq. 1) fitted to wakefulnessand N3 sleep. The input layer has 8100 units, followed by anintermediate layer with 1028 neurons and a two-dimensionallatent space. The next two layers reverse the encoding pro-cess, yielding a matrix C ∗ ij for each z , z pair in the latentspace. The bottom panel presents input matrices C ij (abovediagonal) and their reconstructed versions C ij (below diago-nal) for the model fitted to wakefulness and N3 sleep. To train the network, the errors were backpropagatedvia gradient descent with the purpose of minimising a lossfunction composed of two terms: a standard reconstruc-tion error term (computed from the units in the outputlayer of the decoder), and a regularisation term computedas the Kullback-Leibler divergence between the distribu-tion in latent space and a standard Gaussian distribution.The regularisation term ensures continuity and complete-ness in the latent space, i.e. that similar values are de-coded into similar outputs, and that those outputs rep-resent meaningful combinations of the encoded inputs.[17].We generated 5000 correlation matrices C ij corre-sponding to wakefulness and N3 sleep using the modeldescribed in Eq. 1. We then created 80%/20% randomsplits to obtain training and test sets, and used the train-ing set to optimise the VAE parameters. The trainingprocedure consisted of batches with 128 samples and 50training epochs using an Adam optimiser and the loss FIG. 2. The latent space obtained from wakefulness and N3sleep contains the orderly progression of intermediate brainstates. A) Left: Latent state representation obtained by en-coding the test set (wakefulness and N3 sleep), and the en-coding obtained for the two intermediate states that were notused to train the VAE (N1 and N2 sleep). Right: Latent spacedivided into regions with maximal similarity to wakefulness,N1, N3 and N3 correlation matrices. B) Correlation matricesobtained by decoding an exhaustive exploration of the latentspace variables z and z . function described in the previous paragraph. The latent space encodes the progression of brain statesduring sleep.
The encoding process applied to the wake-fulness and N3 sleep data generated two distinct clus-ters in the latent space (Fig. 2A, left). The encodingof the correlation matrices corresponding to intermedi-ate sleep stages not used to train the VAE (N1 and N2sleep) resulted in separate clusters organised accordingto sleep depth (Fig. 2B). The emergence of a manifold inlatent space where the sequence of correlation matriceswas mapped preserving its continuity suggests that a low-dimensional representation can capture the signatures ofprogressively fading wakefulness.The latent space could be divided into regions cor-responding to wakefulness and all sleep stages, whilealso respecting the ordering of brain states in the de-scent to deep sleep (Fig. 2A, right). We applied thedecoder exhaustively throughout the latent space, ob-taining a pairwise correlation matrix for each z , z pair (Fig. 2B). Next, we computed the structural sim-ilarity index (SSIM) to compare each matrix obtainedfrom the latent space to the matrices corresponding towakefulness, N1, N2 and N3 sleep. SSIM is defined as( µ x µ y +0 . µ x + µ y +0 . )( σ x σ y +0 . σ x + σ y +0 . )( σ xy +0 . σ x σ y +0 . ), where x stands foreach C ij matrix shown in Fig. 2B and y is the average C ij computed for each brain state. The variables µ x , µ x , σ x , σ y and σ xy correspond to the local means, standard devi-ations and co-variances of matrices x and y respectively.SSIM has the advantage of simultaneously weighting theEuclidean and correlation distances between matrices [7].For each z , z pair, we determined the brain state withthe highest SSIM value and constructed the latent spaceparcellation shown in Fig. 2A (right panel). Again, weobserve that the latent space can be orderly divided intoregions corresponding to different sleep depth only fromthe model fitted to wakefulness and N3 sleep. Extreme latent space values predict collapse into struc-
FIG. 3. Latent space variables can be extrapolated to repro-duce increased structure-function coupling as a signature ofunconsciousness. A) An exhaustive exploration of the SSIMbetween the decoded correlation matrices and the empiricalstructural connectivity matrix. High z and low z maximisethis similarity. The red rectangle indicates the range of z and z reproduced in Fig. 2A. B) The empirical structural connec-tivity (left) and the best connectivity matrix reconstructedfrom the latent space (right) with the lower triangular partrepresenting the matrices thresholded at 0.2. FIG. 4. Relationship between latent space variables ( z , z )and the parameters of the homogeneous model ( a , G ). tural connectivity. After mapping the progression of brainstates during sleep into the latent space, we investigatedwhether the variables z , z could be extrapolated to re-produce signatures of other unconscious states. We hy-pothesised that moving past N3 sleep in the latent spacemanifold where the progression of brain states is rep-resented would increase the similarity between C ∗ ij (de-coded correlation matrices computed from the dynam-ics) and J ij (structural coupling matrix). As previouslyshown both in humans and non-human primates[14–16],states of deep unconsciousness are characterised by thecollapse of functional coupling to the underlying anatom-ical connectivity structure.We decoded a wider range of latent space variables andcomputed the SSIM between the output correlations andthe structural connectivity. As shown in Fig. 3A, movingbeyond the N3 region in Fig. 3A (high z , low z ) in-creased the similarity of the generated correlations withthe structural connectivity. Exploring a wider region ofthe latent space, we found the highest similarity betweenempirical ( J ij ) and reconstructed ( J ∗ ij ) structural con-nectivity given by SSIM( J ij , J ∗ ij ) = 0 .
81. Fig. 3A (left)shows the empirical J ij and Fig. 3A (right) shows thebest connectivity matrix reconstructed from the latentspace variables; in both cases the part below the diago-nal corresponds to the matrices thresholded at 0.2. Ashypothesised, moving past the N3 region in latent spacereproduced a well-known signature of deep unconscious-ness. This suggests that the latent space constructedfrom wakefulness and N3 sleep not only represented in-termediate stages, but also captured a manifold wherean ampler range of levels of consciousness can be repre-sented. Mapping the homogeneous model into the latent space.
To gain further understanding concerning how the VAE successfully captured the progression of brain states fromfew parameter combinations, we trained a VAE using anhomogeneous version of the nonlinear coupled oscillatorsin Eq. 1 (i.e. same a for all oscillators), and comparedthe latent space encoding in variables z , z with theparameters a and G [12]. While the resulting correlationmatrices do not reflect those obtained from the empiricaldata, the homogeneous model can be used to gain insighton the mapping performed by the VAE.We trained a VAE with 8000 correlation matrices ran-domly extracted from a set of 1000 matrices generatedwith the homogeneous model. Half of these matrices wasgenerated using a high coupling factor ( G = 0 .
8) and a bi-furcation parameter in the oscillatory regime ( a = 0 . G = 0 .
2) and a bifurcation parameter corresponding tofixed-point dynamics ( a = − . a between-0.02 and 0.02 and G between 0 and 1. For each pa-rameter combination, we found the combination of la-tent space variables that maximised the SSIM betweenboth matrices. In this way, we related each pair ( a, G )in the parameter space with each pair ( z , z ) in the la-tent space. We found that both sets of variables wererelated by approximately linear relationships ( G vs. z , r = − . p < . G vs. z , r = 0 . p < . a vs. z , r = − . p < . a vs. z , r = 0 . p < . a , the latent space approximates a linear trans-formation of the parameters governing the dynamics ofthe coupled oscillators. Discussion.
Several recent studies demonstrated thatlow-dimensional models suffice to capture the large-scalecorrelation structure of neural activity seen during dif-ferent brain states [7, 12]. We went a step further, show-ing that these models implicitly represent different brainstates as points in a low-dimensional manifold. Thiswas established following a constructive process that con-sisted of training a VAE with correlation matrices belong-ing to a reduced set of brain states, and showing that thelatent space represented intermediate states and couldbe extrapolated to produce hypothesised signatures ofdeeper unconsciousness. More generally, we showed thatcomplex nonlinear dynamics depending on two parame-ters could be represented by a latent space that approxi-mated a linear transformation of these parameters. Ourresults suggest that other (e.g. pathological) brain statescould be captured and understood in terms of trajecto-ries within a low-dimensional latent space, with poten-tial applications in diagnosis, prognosis, and data aug-mentation. Generally, we propose that whenever com-plex collective dynamics are suspected to emerge fromfew independent parameters, VAE can be applied to re-construct these parameters as a trajectory embedded ina low-dimensional latent space.Authors acknowledge funding from Agencia NacionalDe Promocion Cientifica Y Tecnologica (Argentina),grant PICT-2018-03103. [1] D. R. Chialvo, Nature physics , 744 (2010).[2] F. Cavanna, M. G. Vilas, M. Palmucci, and E. Tagli-azucchi, Neuroimage , 383 (2018).[3] G. Deco and V. K. Jirsa, Journal of Neuroscience ,3366 (2012).[4] M. Breakspear, Nature neuroscience , 340 (2017).[5] A. Haimovici, E. Tagliazucchi, P. Balenzuela, and D. R.Chialvo, Physical review letters , 178101 (2013).[6] G. Deco, M. L. Kringelbach, V. K. Jirsa, and P. Ritter,Scientific reports , 1 (2017).[7] I. P. Ipi˜na, P. D. Kehoe, M. Kringelbach, H. Laufs,A. Iba˜nez, G. Deco, Y. S. Perl, and E. Tagliazucchi,NeuroImage , 116833 (2020).[8] G. Deco, J. Cruzat, J. Cabral, G. M. Knudsen, R. L.Carhart-Harris, P. C. Whybrow, N. K. Logothetis, andM. L. Kringelbach, Current biology , 3065 (2018). [9] M. L. Kringelbach, J. Cruzat, J. Cabral, G. M. Knudsen,R. Carhart-Harris, P. C. Whybrow, N. K. Logothetis,and G. Deco, Proceedings of the National Academy ofSciences , 9566 (2020).[10] Y. S. Perl, C. Pallacivini, I. P. Ipina, M. L. Kringelbach,G. Deco, H. Laufs, and E. Tagliazucchi, Chaos, Solitonsand Fractals (2020).[11] D. M. Arbabyazd, K. Shen, Z. Wang, M. Hofman-Apitius, A. R. McIntosh, D. Battaglia, V. Jirsa, A. D. N.Initiative, et al. , BioRxiv (2020).[12] B. M. Jobst, R. Hindriks, H. Laufs, E. Tagliazucchi,G. Hahn, A. Ponce-Alvarez, A. B. Stevner, M. L. Kringel-bach, and G. Deco, Scientific reports , 1 (2017).[13] R. B. Berry, R. Budhiraja, D. J. Gottlieb, D. Gozal,C. Iber, V. K. Kapur, C. L. Marcus, R. Mehra,S. Parthasarathy, S. F. Quan, et al. , Journal of clinicalsleep medicine , 597 (2012).[14] J. L. Vincent, G. H. Patel, M. D. Fox, A. Z. Snyder, J. T.Baker, D. C. Van Essen, J. M. Zempel, L. H. Snyder,M. Corbetta, and M. E. Raichle, Nature , 83 (2007).[15] P. Barttfeld, L. Uhrig, J. D. Sitt, M. Sigman, B. Jarraya,and S. Dehaene, Proceedings of the National Academy ofSciences , 887 (2015).[16] E. Tagliazucchi, N. Crossley, E. T. Bullmore, andH. Laufs, Brain Structure and Function221