Interlayer connectivity reconstruction for multilayer brain networks using phase oscillator models

Large-scale neurophysiological networks are often reconstructed from band-pass filtered time series derived from magnetoencephalography (MEG) data. Common practice is to reconstruct these networks separately for different frequency bands and to treat them independently. Recent evidence suggests that this separation may be inadequate, as there can be significant coupling between frequency bands (interlayer connectivity). A multilayer network approach offers a solution to analyze frequency-specific networks in one framework. We propose to use a recently developed network reconstruction method in conjunction with phase oscillator models to estimate interlayer connectivity that optimally fits the empirical data. This approach determines interlayer connectivity based on observed frequency-specific time series of the phase and a connectome derived from diffusion weighted imaging. The performance of this interlayer reconstruction method was evaluated in-silico. Our reconstruction of the underlying interlayer connectivity agreed to very high degree with the ground truth. Subsequently, we applied our method to empirical resting-state MEG data obtained from healthy subjects and reconstructed two-layered networks consisting of either alpha-to-beta or theta-to-gamma band connectivity. Our analysis revealed that interlayer connectivity is dominated by a multiplex structure, i.e. by one-to-one interactions for both alpha-to-beta band and theta-to-gamma band networks. For theta–gamma band networks, we also found a plenitude of interlayer connections between distant nodes, though weaker connectivity relative to the one-to-one connections. Our work is an stepping stone towards the identification of interdependencies across frequency-specific networks. Our results lay the ground for the use of the promising multilayer framework in this field with more-informed and justified interlayer connections.


Introduction
Human brain functioning is widely believed to emerge from neuronal network activity operating at distinct spatiotemporal scales. At the macroscopic level, these functional brain networks may be derived from functional MRI (fMRI), electroencephalography and magnetoencephalography (MEG) [1]. The topology of these networks can be characterized using metrics from the field of network science [2]. This approach has provided a plenitude of new insights into temporal fluctuations of brain states during e.g. cognitive tasks, and revealed common pathways in several neurological disorders [3,4]. Admittedly, the application of network science to EEG and MEG data comes with the challenge of reconstructing frequency-specific functional networks. Once defined, these arguably distinct networks are usually analyzed in isolation, despite the observations that oscillations in different frequency bands might have common neuronal sources and that distinct band limited oscillations may show functional interactions [5][6][7][8]. This probable interdependency calls for studying the frequency-specific functional networks in unison [9,10].
Several studies have demonstrated advantages of multilayer networks [11][12][13] to integrate multivariate biological information [14], especially in the context of neuroscience [9,10,[15][16][17][18][19][20][21][22]. Here, functional networks have been considered as interconnected networks, in which different frequency-specific networks make up different layers with an identical number of nodes and a certain connectivity pattern between layers. However, there is no consensus on how to deal with, or even determine, the interlayer connectivity. This lack of consensus is unfortunate since the choice for interlayer connectivity topology may have a significant impact on the properties of a multilayer network [23,24]. In the case of encephalography, interlayer connectivity can be regarded as a proxy for cross-frequency coupling [8,[25][26][27][28][29][30][31]. Estimation of the cross-frequency coupling with existing metrics remains difficult in practice and it remains an open question whether these can be observed in non-invasive resting-state data [27,32]. Recent studies hint at the presence of cross-frequency coupling in resting-state data [25,31,33]. Here, we push this notion further by incorporating a new data driven approach for the estimation of cross-frequency coupling within the framework of multilayer networks. This approach eventually allows for integration of within band functional connectivity and between band cross-frequency coupling into a single framework.
In the current study, we reconstructed the interlayer connectivity structure for empirical multilayer MEG networks using a recently introduced network reconstruction approach. This approach has been originally developed for a different area of network science, namely the modelling of epidemic outbreaks in a population network [34,35]. As an advantage, the proposed method does not require prior information about interlayer connectivity and estimates interlayer connectivity directly from observed nodal activities using a quantitative description of the nodal activities. In brief, one identifies the sparsest interlayer connectivity matrix given the observed time series of the phase at every node and an a priori defined structural network. Hence, interlayer connectivity is viewed in terms of cross-frequency phase synchronization [36], similarly as in [37], where we extend previous work by considering a multilayer network framework. Our approach requires a quantitative phase description of MEG data. Here, we used two phase oscillator models to provide such a description: a Kuramoto-like network model [38], and a phase oscillator network model derived from the Jansen-Rit (JR) neural mass model [39]. Both models reflect some characteristics of our empirically observed data [39][40][41] where the phase dynamics of the JR model arguably includes a realistic and neurobiologically informed phase interaction function.
Given a lack of ground truth in empirical data, we first evaluated the performance of the interlayer reconstruction approach using in silico data with known ground truth for interlayer connectivity. We subsequently applied the reconstruction approach to, and inferred interlayer connectivity from, empirical MEG data using the multilayer Kuramoto (Kur)-like model. To test the null hypothesis that the reconstructed interlayer connectivity was obtained from a system without underlying interlayer connectivity, we compared the reconstructed interlayer connectivity from genuine empirical data with interlayer connectivity reconstructed from phase randomized surrogates. Apart from inference of interlayer connectivity from empirical data, we also evaluated the reverse of the previous step. We tested whether using the interlayer connectivity as input to the multilayer Kur-like model would result in patterns of intra-layer functional connectivity and layer dependencies similar to those observed in the empirical data. Furthermore, we also reconstructed the interlayer connectivity from the empirical data using the neurobiologically informed phase oscillator network model. The similarity between the reconstructed interlayer connectivity from the Kur-like model and the neurobiologically informed model was estimated in order to evaluate the stability and generalizability of the interlayer connectivity solutions across the different models. The reconstruction approach was applied to experimental resting-state MEG data in order to estimate interlayer connectivity for a two-layered alpha and beta band network, and for a two-layered theta and gamma band network.

Interlayer network reconstruction for multilayer brain networks
An illustration of the general framework is shown in figure 1. While a detailed description of the network reconstruction can be found as supplementary material (https://stacks.iop.org/NJP/23/063065/mmedia), here we briefly sketch the main steps. We consider a two-layered (M = 2) network with identical number of nodes N per layer to model MEG networks, e.g. one layer resembles the phase connectivity in the alpha and the other that of the beta band. Moreover, coupling between nodes within a layer is informed by an N × N anatomical connectivity matrix A with elements a jl , with j, l ∈ {1, . . . N}. Similarly, we considered the symmetric N × N interlayer connectivity matrix B with elements b jl , again with j, l ∈ {1, . . . N}. Recall that this interlayer connectivity characterizes the cross-frequency coupling. This structure gives rise to a multilayer phase oscillator network model, in which the phase dynamics θ L j ∈ [0, 2π) of node j ∈ {1, . . . , N} in layer L ∈ {1, 2}evolve according to Here, H corresponds to the phase-interaction function between the nodes, and the integer scalars n, m represent the frequency ratio between different layers, and c to a global coupling strength. The integer scalars n, m are required since phase locking between oscillators with different intrinsic frequencies occurs in an n:m ratio, i.e. for every n cycles of one oscillator co-occurs with m cycles of the other oscillator. We assume that H is the same for within and between layers. We refer the reader to [37] for an extensive overview on phase interaction functions (or also called coupling functions). The two layers differ in their mean and standard deviation of the respective distributions of natural frequencies g ω L , which we both considered to be Gaussian. By approximating the temporal derivative dθ L j /dt as a finite difference, one can rewrite equation (1) as with discrete time point k and sampling step length T. For every time step k, we can further define the vectors and the scalars that we stack across layers and time steps k. By this, equation (2) obeys the form where B j denotes the jth column of the interlayer connectivity matrix B. Estimation of the interlayer connectivity B j can, hence, be based on solving the set of linear equations in equation (5). However, in the presence of dynamics, noise and/or measurement errors, a mere QR decomposition of Υ j might be inappropriate. Therefore, following [35], we estimated B j by considering the interlayer reconstruction problem as a constrained least absolute shrinkage and selection operator problem (LASSO problem [42,43] Here, ρ j > 0 refers to a regularization parameter and u to the all-one vector, . . . || 1 and . . . || 2 denote the one-norm and two-norm, respectively. The first term on the right-hand side of equation (6) corresponds to the mismatch between empirically derived phases and the multilayer phase oscillator network model, while the second one tunes the sparsity of the solutionB j by the regularization parameter ρ j . We apply κ-fold cross-validation using training and test sets (split data into half) to estimate the value of the parameter ρ j within a pre-defined range. This cross-validation yields the minimum mean squared error between the phases obtained for the phase oscillator network model and the measured data. This implementation allows for non-zero elements on the diagonal of the estimated B matrix, i.e. it may include one-to-one interlayer coupling. General framework and phase interaction function Jansen-Rit model. Phase data from two layers, be it empirical data or simulated data θ L j , is used as input together with a given structural network A. Here, we assume that the quantitative model (Kuramoto-like or Jansen-Rit) provides a good description of the phase data that is used as input, or in other words, that the phase interaction function H is accurate. This leaves the interlayer connectivity B as the only unknown term in the equation, for which the equation is subsequently solved. We either use the Kuramoto-like description (H is a sine function), or we derive H from the Jansen-Rit model using a phase reduction formalism. Panel (B) shows the phase interaction function H derived from the Jansen-Rit (JR) model with the external input P = 150, i.e. with the model operating in the oscillatory regime.

Phase oscillator models
We consider two seminal phase oscillator models to express the phase interaction function H in equation (1). For the first model, we examined resulting in a multilayer Kur-like model [44,45]. Since this 'simple' form arguably lacks a proper neurobiological underpinning, we also consider a phase interaction function based on neural mass dynamics. In more detail, we derive a phase interaction function from the JR model [39], which is known to generate realistic MEG/EEG oscillations [46][47][48]. This model describes the evolution of the synaptic activity and the firing rate for three interconnected neuronal populations (i.e. an inhibitory, excitatory and pyramidal neuronal population) in terms of six coupled first-order differential equations dx/dt = F (x), with x = [x 1 , . . . , x 6 ] T provided in the supplementary material. The phase interaction function can be derived via the phase reduction technique based on weakly coupled oscillator theory described by Ashwin and coworkers [49], see also [50] for recent literature on phase reduction techniques. The two main assumptions of weakly coupled oscillator theory are: (1) weak coupling (i.e. a perturbed system stays close to the intrinsic limit cycle); (2) (nearly) identical oscillators. One first computes the phase response curve for a single population that displays limit-cycle oscillation. The phase response curve Q = q 1 , . . . , q 6 describes the response of a limit-cycle to a small perturbation and can be given by the periodic solution of the adjoint equation [49] dQ dt where DF corresponds to the Jacobian of F, evaluated along the limit-cycle x. The expression · · · denotes the Euclidean inner product, ω is the frequency (which here is equal to 2π divided by the period of the limit-cycle). Equation (8) can be solved by backward integration in time [51]. Rewriting the system in terms of phases by introducing x (t) = x θ/ω and considering that all oscillators lie on the same limit-cycle of a system, one can treat their interactions as small perturbations μ (x, t) = [μ 1 (x 1 , t) , . . . , μ 6 (x 6 , t)] acting on the intrinsic dynamics of the oscillators that hence becomes  We restrict our analysis to phase response components of the excitatory population q 2 , and, hence, set μ (x) = [0, μ 2 (x 2 ) , 0, 0, 0, 0]. This choice is justified by the fact that interaction between neuronal populations is mediated by excitatory neurons. By introducing a rotating phase ψ = θ − Tt/Δ, and when assuming coupling c to be very small, ψ evolves very slowly, dψ/dt is approximately zero, and averaging over one period T of the rotating phase gives [49] and δ = T/Δ − ω. H (ψ) corresponds to the phase interaction function, which has to be evaluated numerically. For this purpose, one can express H as a Fourier series (up to the 3rd term) where d n and e n are the Fourier cosine and sine coefficients, respectively, and ε denotes the remainder.
Using parameters as in [52], with the external input P = 150 (see supplementary material), the JR model operates in the oscillatory regime, for which we illustrate the phase interaction function in figure 1(B). As can be observed from figure 1(B), the derived phase interaction function will not necessarily result in a higher order periodic function, but could also be a shifted sine function, much like the so-called Kuramoto-Sakaguchi model [53]. In addition, the shape of the obtained phase interaction function is very similar to a phase interaction function that resulted from a data driven approach using EEG data [54]. The bifurcation diagrams for the two employed phase interaction functions can be found in the supplementary material.

Simulations with ground truths for interlayer connectivity
An overview of our analysis is provided in figure 2. We evaluated the interlayer network reconstruction by providing input from simulations with ground truth for interlayer connectivity. We generated phase time series using the network dynamics described in equation (1) with the phase interaction function given either by equation (7) or (11). We also added uncorrelated white Gaussian fluctuations w L j with zero mean and variance σ 2 to the right-hand side of (1) to account for dynamic noise; see, e.g., [55]. We simulated Equation (12) is solved using an Euler-Maruyama scheme with a sampling time T = 0.01. For the Kur-like model, every node in layer L had a natural frequency ω L j randomly drawn from a Gaussian distribution, centred around ω L 0 , with ω 1 0 = 1 and ω 2 0 = 2. In general, two coupled phase oscillator networks with (symmetric) unimodal frequency distributions are equivalent to a single phase oscillator network with a bimodal frequency distribution [45]. The ratio ω 2 0 /ω 1 0 = 2 was also used for the neurobiologically informed (JR) model. For the JR model we set ω L j = ω L 0 as we assumed identical node dynamics for every layer in this case; otherwise all settings agreed with the Kur-like model. These identical node dynamics were the result of the underlying assumptions of weakly coupled oscillator theory, not applicable to the Kur-like model.
For a first set of simulations (figure 3(A)), we tested the estimation accuracy versus network sparsity or link density (i.e. the ratio of the number of links to the maximum possible number of links given the network size). We generated an unweighted and connected structural connectivity matrix A using the Erdös-Rényi random model with an approximate link density of 0.2 and only included network realisations that were connected [56]. The network size agreed with our empirical data (N = 78; see below in section 'reconstruction of interlayer connectivity for empirical MEG networks'). In order to avoid that the network reconstruction of B was biased by network A, we constructed interlayer connectivity B with a different topology. Hence, We computed B using a range-dependent random graph model [57] using the Contest toolbox for Matlab [58]. The probability of two nodes i and j in different layers being connected was given as p ij = αλ |j−i|−1 ; j, i ∈ {1, . . . , N}, with λ = 0.9 and α was varied over the interval [0.5, 1.75], resulting in networks with different link densities. Equal increases in α do not result in a linear increase in link densities (figures 3(A)-(D)). Furthermore, the interlayer coupling matrix was constrained to be symmetric by setting B → B + B T /2. In a second set of simulations, we tested the estimation accuracy versus the variance σ 2 of the Gaussian white noise in the dynamics of equation (12). We used the same range-dependent random graph model to generate B, with α = 0.9. We varied the variance of the noise in the interval [0.05, 0.5]. The other parameters for the network models for all Kuramoto-like simulations were: c = 1, m = 1 and n = 2, sampling time T = 0.01 and the number of observations K = 10 × N 2 .
In order to evaluate the accuracy of our interlayer connectivity reconstruction algorithm, we employed three error metrics. The first one was the false positive rate (FPR) for the interlayer links, and equals the fraction of node pairs i, j for which b ij = 0 butb ij > 0, whereb ij refers to a link between node pairs i, j in the estimated interlayer connectivity matrixB. The second one was false negative rate (FNR) for the interlayer links, which equals the fraction of node pairs i, j for which b ij > 0 butb ij = 0. FPR and FNR were normalised (divided) by the link density of the ground truth B matrix in order to provide a sense of the relative error compared to the number of links in the network. The third metric was the relative true positive deviation (RTPD), which measures the deviation on the link weights and equals b ij −b ij /b ij averaged over all node pairs i, j for which b ij > 0 andb ij > 0. The RTPD captures different information compared to FNR and FPR and estimates whether the link weights of the reconstructed links are similar to the corresponding link weights in the ground truth network. Figure 3(A) shows the effect of different link densities on FPR and FNR for both the multilayer Kur-like model and the JR model. There was an increase in FPR for the Kur-like model with increasing link density, whereas the JR model stayed close to zero false positives for increasing link density ( figure 3(A)). The reverse was true for FNR, i.e. FNR for the Kuramoto model remained close to zero, while for the JR model there was an increase in the FNR for increasing link densities. RTPD showed a different picture ( figure 3(B)): for both models the RTPD remained stable for different link densities, with lower errors for the Kur-like model than for the JR model, i.e. there was (in this case) a systematic underestimation (figure 3(E)) of the link weights for the JR model (see figure 3(E)). When increasing the noise in the system with a fixed link density for the ground truth interlayer connectivity B of 0.14 (α = 0.9), FPR for the Kuramoto based simulations slightly increased ( figure 3(C)). Again, RTPD (figure 3(D)) was hardly affected, with the Kur-like model yielding a smaller error than the JR model over the whole range of noise levels. Figure 3(E) shows an example of a ground truth interlayer connectivity matrix and the reconstructed B matrices for the Kur-like and JR model, which show quite high accuracy of estimation of the ground truth links. On average, smaller link weights (systematic underestimation) with more false negatives could be observed for the JR model, while more false positives and a more accurate estimation of the link weights were present for the Kur-like model.

Reconstruction of interlayer connectivity for empirical MEG networks
Empirical MEG and diffusion weighted imaging data We used MEG data from the Human Connectome Project [59,60], consisting of resting-state MEG data from 89 healthy subjects. Every subject underwent three separate recording sessions. The three separate recording sessions were used as training and validation datasets for separate analysis steps (see figure 2).
We refer to [61] for details of the pre-processing pipeline for this dataset. The data have partly been provided pre-processed [59], after passing through a pipeline to remove any artefactual segments of time from the recordings. We performed additional processing steps for source localization. An atlas-based beamforming approach was adopted to project MEG sensor level data into source-space [62]. The cortex was parcellated into 78 cortical regions according to the automated anatomical labelling (AAL) atlas [63] and the centroid voxel for every region of interest was extracted to serve as representative voxel for every region [64]. Pre-computed single-shell source models are provided by the HCP at multiple resolutions [65], registered into the standard co-ordinate space of the Montreal Neuroimaging Institute. Data were beamformed with depth normalization onto centroid voxels using normalised lead fields and estimates of the data covariance. Covariance was computed based on broad band data with a time window spanning the whole experiment [66]. Regularization was applied to the data covariance matrix using the Tikhonov method with a regularization parameter equal to 5% of the maximum eigenvalue of the unregularized covariance matrix. Dipole orientation was determined using a non-linear search for optimum signal to noise ratio [67]. This complete process resulted in N = 78 electrophysiological timecourses, each representative of a separate AAL region.
The anatomical network data was also obtained from the Human Connectome Project. Details of the data collection and processing pipeline for this diffusion weighting-based anatomical network can be found in [61].
Reconstruction of interlayer connectivity from empirical data using the multilayer Kuramoto model Since neither model clearly outperformed the other one for the numerically simulated data (figure 3), we first used the arguably simpler multilayer Kur-like model for the following analyses. An overview of the analysis for this part of the study can be found in figure 2(B). We first considered interlayer connectivity between the alpha and beta band. Concatenated phase data from all subjects from session one together with the empirical connectome A were fed into equations (5) and (6) with c = 1, m = 1 and n = 2. The choice for c = 1 is justified since we have no knowledge of the optimal value beforehand and a too small value for c could erroneously result in the absence of phase synchronization (see supplementary material on phase interaction functions. Thus, we assumed that our quantitative description of the phases is an appropriate model for the empirical phase data. We also used surrogate phase data to test the outcome in case of no  genuine underlying phase synchronization in the data. Surrogate time series were reconstructed by non-uniform phase randomisation of the original empirical source-reconstructed time series, i.e. for every nodal time series a different sample of Gaussian white noise was added to the Fourier phases, followed by the inverse Fourier transform. The instantaneous phases were subsequently extracted (as described above) before the interlayer network reconstruction via equations (5) and (6). The output of the algorithm applied to empirical phase data yielded a very sparse interlayer connectivity matrix, with predominant, and strong, one-to-one connections (diagonal of B, see figure 4(A)), and a few weak connections between distant nodes. Application of our interlayer network reconstruction to surrogate data resulted in an even sparser and almost empty B matrix ( figure 4(B)), which did not resemble the interlayer connectivity as reconstructed from genuine experimental data. This result indicates that one-to-one coupling as interlayer connectivity is not necessarily the minimal or sparsest solution of the interlayer network reconstruction problem, and that reconstructed interlayer connectivity for the empirical data could not have been obtained from a system without underlying interlayer connectivity. Application of our interlayer network reconstruction to a different session for all subjects revealed similar interlayer connectivity with strong one-to-one connections (diagonal of B, see figure S2), however, with a slightly different set of much weaker connections between distant nodes, i.e. off-diagonal elements in the B matrix. Thus, these weaker estimated connections show lower between session consistency.
However, one-to-one coupling for alpha to beta is influenced by spectral leakage between frequency filtered time series due to the small overlapping decaying tails after band pass filtering for these adjacent frequency bands (see supplementary material). Nevertheless, even in the absence of spectral leakage between frequency filtered time series, one-to-one coupling remained nonzero for alpha to beta band connectivity (see supplementary material), and therefore interlayer connectivity is not merely an artefact of spectral leakage. We adhered to the use of classical frequency bands since these are widely accepted in the field. Although readers should be aware of this potential bias due to band pass filtering and the alternative is to use adjusted frequency bands.
Simulations of multilayer networks informed by reconstructed interlayer connectivity to explain empirical intra-layer MEG networks An overview of the analysis steps for this section is given in figure 2(C). Apart from inferring the interlayer connectivity from the empirical phase data in the previous section, it remains to be elucidated whether this interlayer connectivity would result in empirically observed intra-layer MEG connectivity and correlations between intra-layers. We therefore simulated equation (12) using the Kur-like model with the reconstructed interlayer connectivity B from the previous section and the empirical connectome A, and analysed how similar the simulated intra-layer functional connectivity was to the empirically determined functional connectivity. In the previous section we used phase information from MEG data as input variable to find B. Now we use the reconstructed B matrix as input in order to simulate the phases θ L j (that is, we are not solving equations (5) and (6), but just simulating equation (12)). We simulated equation (12) using the Kur-like model for a range of c. For every value of c, we computed the (i) Kuramoto order parameter for each layer, r L = 1 N N j=1 e iθ L j ; (ii) and evaluated the fit between simulated and empirical functional connectivity for every layer (W L,sim vs W L,emp ) in terms of the variance of the empirical data explained by the simulations (expressed as adjusted R 2 ); (iii) we estimated functional connectivity for each layer using the PLV, resulting in W L,sim , and subsequently computed a correlation between the functional connectivity patterns of the two intra-layers, ρ = corr W 1,sim , W 2,sim using the Spearman correlation coefficient with the upper-triangular part of these symmetrical matrices.
As shown in figure 5(A), there was a rapid transition for the Kuramoto order parameter from weak overall phase synchronization to strong overall phase synchronization without a clear plateau for lower coupling values. The absence of such a plateau for smaller coupling values, as typically seen for single layer Kuramoto models [41], might be due to the non-zero interlayer connectivity. The best fit of simulated PLV connectivity matrices with group averaged empirical PLV connectivity matrices was observed for intermediate coupling values (around c = 0.5, figure 5(B)), i.e. at the transition between high and low r L . Adjusted R 2 values of around 0.3 to 0.4 indicate that to some extent multilayer Kuramoto network models can approximate patterns of functional connectivity as observed in empirical MEG data. There was a slightly better fit for the alpha band intra-layer network than for the beta band intra-layer network. Correlation between intra-layers for empirical data was ρ = 0.62. The coupling for which we obtained the best match of simulated and empirical ρ was found for c = 0.6 ( figure 5(C)).
Validation of empirical interlayer connectivity reconstruction using the Jansen-Rit model We next applied our network reconstruction algorithm to concatenated phase data from all healthy adults, this time from session three, together with the empirical connectome A (see figure 2(D)), i.e. empirical data was fed into equations (5) and (6). The phase interaction function derived from the JR model operating in the oscillatory regime was now used for equation (1). Analysis was otherwise identical to that outlined in figure 2(B). Again, we also reconstructed interlayer connectivity for phase randomised surrogate data. As shown in figure 4(C), we found a strong one-to-one interlayer coupling similar as for the Kur-like model (figure 4(A)) but with fewer off-diagonal elements in B. The link weights on the diagonal of the reconstructed interlayer connectivity matrix B were close to one, as was the case with the Kur-like model. When applied to surrogate data, a very sparse interlayer coupling matrix with a few off-diagonal non-zero elements was obtained ( figure 4(D)). This B matrix obtained for surrogate data did not resemble interlayer connectivity from genuine phase data. In addition, application of our interlayer network reconstruction to a different session for all subjects revealed similar interlayer connectivity with strong one-to-one connections (diagonal of B, see figure S2). However, the analysis resulted in a slightly different set of much weaker connections between distant nodes, i.e. off-diagonal elements in the B matrix. Again, these weaker estimated long distance connections show lower between session consistency.

Reconstruction of multilayer networks using a different set of frequency bands
So far, we reconstructed a two-layer network based on the alpha and beta band. Previous work has also demonstrated cross-frequency coupling between theta and gamma bands [25]. Similar to figure 4, we reconstructed the interlayer connectivity matrix based on a two-layered network consisting of theta and gamma band networks. Concatenated phase data, obtained from band-pass filtering into the theta and gamma bands, from all subjects from session one together with the empirical connectome A were fed into equations (5) and (6) with c = 1, m = 1 and n = 6. As outlined above, we also used surrogate phase data to test the outcome in the absence of genuine underlying phase synchronization in the data. Results for both the multilayer Kur-like model and the JR model are illustrated in figure 6. The estimated interlayer connectivity matrices between theta and gamma band layers show many more off-diagonal elements compared to the alpha-beta band case (compare figures 4(A) and (C) with 6(A) and (D)), and as well in comparison to surrogate data. This indicates the presence of long-range interlayer connections between regions. For the mean connectivity for every brain region (figures 6(B) and (E)), we found the main off-diagonal connections to involve occipital and fronto-parietal areas. The Spearman correlation between the estimated interlayer connectivity matrix from the Kur-like model and the JR model was R = 0.9, p < 0.001. Link weights from the interlayer connectivity matrix obtained from the Kur-like model were on average higher than link weights in the interlayer connectivity matrix obtained from the JR model (Mann-Whitney test, p < 0.001, Z = 4.4). The latter results are in line with the plots obtained from the simulations (figure 3(E)), where the JR based model was shown to underestimate the link weights from ground truth networks, and the Kur-like model prone to false positives.

Discussion
We reconstructed the interlayer connectivity for MEG networks using a recently developed network reconstruction algorithm for epidemic spreading models [35]. We used two types of phase interaction functions, namely a sine function for the Kur-like model, and a phase interaction function obtained from weakly coupled JR neural mass models operating in the oscillatory regime. We demonstrated that the network reconstruction approach accurately captures simulated interlayer connectivity for both types of phase oscillator network models, and is robust to different levels of dynamic noise in the phase data and increasing levels of link density. Application to empirical MEG data revealed that, when alpha and beta bands were considered, empirical interlayer connectivity was dominated by one-to-one connectivity between layers, which was consistent for both phase oscillator models for different MEG recording sessions. However, when theta and gamma band layers were considered, there were also widespread long distance connections between regions, i.e. strong off-diagonal elements in the B matrix.
The main result of the current analysis is that the topology of interlayer coupling strongly depends on the combination of frequency bands. While alpha-to-beta band interlayer connectivity was dominated by one-to-one interlayer coupling, theta-to-gamma band interlayer connectivity, in contrast, also showed strong connectivity between multiple distant nodes. This result is in line with previous work that has reported on long-distance cross-frequency connections [70]. The current dominant connections involving the parieto-occipital areas in the theta-gamma band were also observed for cross-frequency amplitude-amplitude coupling in resting-state data MEG data [9,10]. For example, recent studies reported on the presence of hippocampal-prefrontal theta-gamma coupling in local field potentials recorded in animals during a spatial working memory task [71] and on the presence of frontotemporal theta-gamma coupling during working memory tasks in adults [70]. These findings are assumed to relate to the tendency of interaction between neural network motifs that generate theta and gamma oscillations [5]. Our results appeared consistent for both phase oscillator network models, with the difference that for theta-gamma coupling the use of the Kur-like model resulted in more and stronger off-diagonal elements compared to the JR model. Results from our simulations suggest that the observation of fewer connections between distant nodes for the JR model in the empirical theta-gamma coupling could be due to false negatives for this model, an underestimation of link weights for this model, and due to false positives for the multilayer Kur-like model. By using surrogate data, we also showed that a one-to-one coupling matrix for the alpha and beta band is not the minimal solution of the algorithm and does not emerge from a system with no underlying functional connectivity.
In the current paper, we pave the way towards standardization of the usage of multilayer brain networks, which is important for comparison between studies. Although the use of multilayer brain networks is still in its infancy, many different methodological options have already been published [14,17], yielding promising results. Several studies choose one-to-one connectivity for the multilayer network analysis [20,72] and were able to show disease induced effects in centrality of brain regions and differences indisrupted core-periphery structure in Alzheimer's disease. Also a portion of literature opted for all-to-all connectivity [10] and demonstrated disease induced effects for interlayer coupling for patients suffering from major depression [73,74].
Though the estimated interlayer connectivity matrix is considered to be an approximation for cross-frequency functional connectivity, we stress that the estimated interlayer connectivity and cross-frequency coupling are not equivalent. Cross-frequency functional connectivity is usually the result of estimation of a pairwise statistical dependency, whereas the current interlayer connectivity estimation is the result of a multivariate assessment and optimization using an a priori defined structural network and phase interaction function. We therefore cannot draw firm conclusions about the existence or non-existence of cross-frequency coupling between pairwise nodes for two reasons: (i) we only considered phase connectivity [31,75], whereas cross-frequency coupling is usually assessed in terms of phase-amplitude coupling [6,[76][77][78][79] or sometimes even amplitude-amplitude coupling [9,10]; (ii) our algorithm merely estimates the sparsest interlayer connectivity matrix that can explain the data, it does not assess whether a specific cross-frequency connection is neurobiologically plausible or not [25,26]. More importantly, the goal of the current work was to reconstruct interlayer network topology rather than the estimation of pairwise strength of cross-frequency coupling. The advantage of the current approach is that it makes use of a priori information of the structural network and that its multivariate assessment of interlayer connectivity is also less influenced by several factors (e.g. non-sinusoidal signal properties) known to give rise to spurious estimates of pairwise cross-frequency coupling [27].
There was a difference in reconstruction accuracy in terms of FPR and FPN for the JR and Kur-like model. The trade-off between decreasing FPR and increasing FNR translates into setting the regularization parameter ρ j : a greater value of ρ j results in more zero elements of B j and, hence, a lower FPR and higher FNR; a smaller value of ρ j results in more non-zero elements of B j and, hence, a lower FNR and higher FPR. We set the value for ρ j by cross-validation with the main objective to maximise the fit of the phase oscillator model (1) to the data. Thus, the trade-off between FPR and FNR is done implicitly. In fact, figure 3 shows that for the Kur model a smaller value for ρ j (which results in a higher FPR and a lower FNR) has more predictive power, as assessed by cross-validation, than a larger value for ρ j . For the JR model, the opposite holds.
Some other aspects of our work also warrant further discussion. Firstly, we only included two layers in our multilayer network model, whereas usually more than two layers can be estimated from MEG data. For example, delta (1-4 Hz) to low gamma (30)(31)(32)(33)(34)(35)(36)(37)(38)(39)(40)(41)(42)(43)(44)(45)(46)(47)(48) networks are typically estimated, which would entail at least five layers. However, analytical work and bifurcation analysis on phase oscillator networks has incorporated not more than three multimodal frequency distributions [45]. Therefore, without theoretical background it is difficult to predict model behaviour of phase oscillator models given the width and distance between the multimodal frequency distributions. In addition, it is difficult to assess a priori whether the solution to the multilayer phase oscillator model equations is expected to be stable or unstable. We therefore restricted our analysis to two layers. Secondly, previous literature demonstrated that the inclusion of conduction delays increases the power of neuronal models to explain empirical patterns of intra-layer functional connectivity [40,80,81]. However, we neglected conduction delays in the dynamics of our model. Even without delays, the match with empirical MEG data was fairly good. Yet, future work may examine if the inclusion of delays would indeed lead to even better descriptions of multilayer functional networks. Thirdly, since our aim was to test the null hypothesis of no underlying interlayer connectivity, we did not to apply a uniform phase randomization to preserve the linear correlations or static connectivity between timecourses [82,83]. Fourthly, we restricted our analysis to resting-state (task-free) data, implying that no conclusions can be drawn from the current analyzes with regards to the topology of interlayer connections in task-based MEG networks. Fifthly, the presence of periodic signals in the alpha and beta band in empirical data [84] justifies the treatment of phase synchronization in these canonical frequency bands and further treat these as layers in a multilayer network framework. However, recent interest in treating neural power spectra in terms of periodic and aperiodic components challenges the use of (only) canonical frequency bands [85]. It remains an open question how this new approach would influence connectivity estimation, and hence how to reconstruct interdependent networks. An alternative way to treat interdependent electrophysiological networks is to consider connectivity for the aperiodic part of the spectrum as one layer, and consider connectivity for all periodic components on top of the aperiodic part of the spectrum as separate layers. Sixthly, it remains an open question whether pairwise interactions form the adequate building blocks to assess (multilayer) network topology and recent studies have demonstrated the potential role of higher order coupling in this context [86][87][88]. Lastly, it is an open question whether the assumption of identical phase interaction functions for within and between layer interactions is justified.
We demonstrated the robustness of an interlayer network reconstruction algorithm in simulated brain networks. Application to empirical multilayer brain networks revealed that interlayer connectivity is dominated by one-to-one coupling for a two-layered alpha and beta band network, and revealed additional widespread long distance interlayer connections for a two-layered theta and gamma band network. Therefore, in future resting-state empirical multilayer network analyzes, a one-to-one coupling, i.e. a multiplex network description, is only justified for specific combinations of frequency bands. For other scenarios, one may require a complete multilayer network description.