## Synthetic Networks

Three modules in a network of 100 nodes were simulated as follows in 10 consecutive time windows. Nodes were initially assigned to the same module across all the time windows (30% of the nodes were assigned to module A, 50% to module B, 20% to module C). Then, 40% of the nodes were chosen at random to be fast/slow/random oscillators (15%, 15%, 10% respectively) to simulate diverse physiological properties of brain hubs. Fast/slow/random oscillators switched between their module and another random module every 5/3/1 time windows, respectively.

Once defined the modular time-varying structure, time series were created as follows. First, for each time window, a prototype matrix of synthetic functional connections was simulated as a 100x100 symmetric Lehmer matrix. Since in the Lehmer matrix the values of the functional connections increase with proximity to the diagonal and with progressive node numbers, such matrix allows the generation of modules with different degrees of intrinsic connections, thus, avoiding bias related to connection strengths. A time series was simulated for each node from the starting correlation structure and, for each module, the nodes' time series were selectively time-shifted to create modular structures.

Within each node, the time series of 100 voxels were simulated to simulate realistic voxels within a parcel. Initially, for each node, the node's time series was assigned to each voxel within the node. Then, two sources of noise were simulated to represent physiologically-unrelated noise (unstructured noise) and the "interference" of nearby biological units in the node signal (structured noise). Unstructured noise was simulated - for each voxel - as a random time series that was added to the voxel signal in a weighted manner. Structured noise was simulated by initially choosing a random [2 to 5] number of attractors on the external surface of the node. A random time series was simulated for each attractor. Then, each voxel in the node suffered each attractor's interference depending on its closeness to the attractor. Finally, to mirror common procedures in neuroimaging, each node time series was obtained by averaging the time series of the voxels within it. Different weights were considered for the two noises (low, medium, high) to investigate the impact of the signal-to-noise ratio to the module detection across the two methods.

Once noisy time series were synthesized, the generalized Louvain function was used to detect modules across time points (http://netwiki.amath.unc.edu/GenLouvain/GenLouvain). The generalized Louvain function was implemented both in the "standard" form (ω = 1) as well as using the probabilistically weighted method (PW). The whole procedure was repeated for 100 independent cycles so that, in each cycle, a new modular structure was randomly created. The three levels of structured/unstructured noises were introduced for each cycle. Other values controlling the procedure, such as γ and tau (τ), were set to default values (1) for simplicity. Analyses were performed using MatLab version 2019a (The MathWorks). The weighted approach is described in the next paragraphs. The procedure for obtaining synthetic networks is schematised in Fig. 1. To assess the independence of methods’ performances on the number of modules in the network, the whole procedure was repeated using five isomorphic modules of equal size (20% of the nodes assigned to each module).

## Probabilistically Weighted (PW) Multilayer networks

To implement control over nodal temporal features of the network, we introduced a matrix of probability weights W. Values in W were controlled by the temporal resolution parameter β and by the temporal features of the nodes. As probability weights, values in W basically represent the probability of nodes to be self-connected among different layers of the network. Thus, with temporal multilayer networks, W is a N by T-1 matrix in which T is the number of time windows and N is the number of nodes in the network. By contrast, if working with task multilayer networks, the matrix W would be a K by N matrix in which K represents within-node cross-condition interconnections and is calculated as K = Nc! In which Nc is the number of experimental conditions (e.g., in a classic two by two factorial psychology paradigm, there would be four conditions, and K = 4! = 24).

To control the strength of cross-layer connections and to mirror the weight of the null model for the structural resolution parameter, values of W were defined using a probabilistic approach. Initially, a set of Levy alpha-stable distributions, Ω, was generated. Stable distributions rely on four parameters, allowing complete control over the shape and the scale of the distribution. To reflect correlation (functional connectivity), the values of W (probabilities) were sampled within the interval [0 1]. Thus, the location parameter of the distributions, µ, was set to 0.50 and the scale parameter, σ, was set to 0.75. The first shape parameter, α, was empirically set to 0.40, while the second shape parameter, β, was modulated within the range [-1 1] to obtain a series of distributions in which the head and the tails are arranged progressively in the interval [0 1]. In other words, a set of probability distributions (Ω) was created in which the ratio of low versus high probabilities was controlled through the parameter β: with higher, positive values of β, the head and the fat tail of a specific distribution in Ω are located toward lower probability values (lower probabilities of cross-layer connections); instead, with lower, negative values of β, the head and the fat tail of a specific distribution in Ω are located toward higher probability values (higher probabilities of cross-layer connections). Values used for β in this study were: [-.75 − .25 0 .25 .75].

Then, the distributions in Ω were used to generate values in W based on the coherence between time intervals. More specifically, we used spectral coherence, also known as magnitude-squared coherence38,39, across consecutive time windows. Since coherence estimates the consistency of amplitude and phase between signals across frequencies, it is biologically appropriate for a node exhibiting coherent signal fluctuations in two adjacent time windows (or in two experimental conditions) *i* and *j* to have higher probability to be connected – to be “self-similar” – between the layers *i* and *j* (i.e., it occupied the right tail in the stable distribution). Using this procedure, the temporal resolution, β (the β shape parameter of the stable distributions), regulated the probabilities of connections among adjacent time windows in a way driven by the similarity of the node’s signals in the two time intervals. The procedure for estimating probabilistically weighted multilayer modules is illustrated in Fig. 2.

The performance of the two methods was assessed by calculating the degree of similarity between the true simulated modular structure and the estimated modular structure through two commonly used indexes for calculating correlations across modular structures, namely the adjusted mutual information index (AMI) and the Rand coefficient.

## Application on resting-state data

The procedure described above for detecting PW multilayer networks, together with the standard ω = 1 procedure was implemented also on resting-state data from a sample of 39 healthy participants (19 females and 20 males, aged 23 ± 2; 35 right-handed and 4 left-handed) without a history of psychiatric or neurological disease and contraindications for MRI scanning participated in the experiment. The experiment was approved by the local ethics committee (Comitato Etico per la Ricerca Biomedica delle province di Chieti e Pescara) and by the Department of Neuroscience, Imaging and Clinical Sciences (DNISC) of the G. d'Annunzio University of Chieti-Pescara.

All participants had a normal or corrected-to-normal vision and provided written informed consent before taking part in the study in accordance with the Declaration of Helsinki (2013).

## Resting-state Data acquisition

Each participant performed two consecutive task-free fMRI runs, each consisting of 376 volumes. The participants were instructed to watch a white fixation cross on a black screen without performing a cognitive task. Each run lasted approximately 7.5 minutes. Functional images were acquired using a Philips Achieva 3T scanner installed at the Institute for Advanced Biomedical Technologies (Gabriele d’Annunzio University, Chieti-Pescara, Italy). Whole-brain functional images were acquired with a gradient echo-planar sequence using the following parameters: repetition time (TR) = 1.2 s, echo time (TE) = 30 ms, field of view = 240x240x142.5 mm, flip angle = 65°, in-plane voxel size = 2.5 mm2, slice thickness = 2.5 mm. A high-resolution T1-weighted whole-brain image was also acquired after functional sessions using the following parameters: TR = 8 ms, TE = 3.7, FoV = 256x256x180 mm, flip angle = 8°, in-plane voxel size = 1 mm2, slice thickness = 1 mm.

## Resting-state Data analysis

Functional connectivity was calculated as the correlation among average parcels' timeseries for a total of 418 nodes using cortical and subcortical atlases from Joliot and colleagues45 plus the cerebellar atlas from Diedrichsen and colleagues46. These atlases were chosen because of the integrative method used to define both cortical and subcortical parcels without lateralization biases. The functional connectivity values, that is, the edges of the connectomes, were obtained using the z Fisher transform of the Pearson correlation among pre-processed time series and were used to create individual (subject-specific) weighted graphs. The modular architecture was visualized using BrainNet Viewer (www.nitrc.org/projects/bnv/). The PW method and the standard ω = 1 approach were both implemented to detect time-varying modular structures across 10 consecutive time windows of equal length.

We tested if brain-behavioral resting-state associations were differentially detected by the standard ω = 1 and by the PW method. Thus, nodal metrices of flexibility, recruitment, integration, and promiscuity were estimated for each subject using both the ω = 1 and the PW methods. The behavioral variables were represented by five latent, orthogonal psychosis-relevant factors that have also been used in previous works to study the sense of self and the sense of agency32,47,48. We investigated the relation between dynamic multilayer metrics and psychosis-related personality traits with a focus on factors related to psychotic disorders49,50. These five factors were reliably (Tucker congruence coefficient = 0.95) derived from a factor analysis with orthogonal (orthomax) rotation on a series questionnaires that were administered to a larger sample (N = 101) including the present fMRI sample (see references32,47,48 for the complete procedure) The questionnaires measured basic as well as psychosis-relevant personality traits: the Big-Five Questionnaire (BFQ short version51), the tolerance of uncertainty scale (IUS-1252), the Schizotypal Personality Questionnaire (SPQ53), the Community Assessment of Psychic Experience (CAPE54), and the State-Trait Anxiety Inventory (STAI255). The five factors showed the highest loadings in the following subscales: factor 1) STAI2, IUS-inhibitory, CAPE-depression, BFI-neuroticism (negative affect); factor 2) SPQ-cognitive perceptual, CAPE-positive, CAPE-negative (psychosis-like experiences); factor 3) IUS-prospective, BFI-extraversion, BFI-agreeableness, BFI-openness (sociality); factor 4) SPQ-interpersonal, SPQ-disorganizational (schizotypal); factor 5) BFI-conscientiousness.

Brain-behavioral relationships were investigated using whole brain, univariate mixed-effects models in which the five psychometric factor scores were used as continuous predictors for each multilayer metric of interest; random intercepts and slopes were added at the individual node level to allow for ROI-specific weighting of brain-behavior correlations; random slopes were added at the subject level. When significant associations were detected, we located the brain structural correlates using best linear unbiased predictors (BLUPs) to generate nodal conditional expectation (ICE) plots and to highlight nodes with the highest contribution.