Emotion brain network topology in healthy subjects following passive listening to different auditory stimuli

Study designs

The study was conducted using a cross over experimental design in which all subjects were exposed to a series of five auditory stimuli in specified orders. The order of sequences in which auditory stimuli were presented was generated in such a way that each stimuli must occur once within each subject, and each stimuli must have the same number of occurrence. This was done to control the period and order effects in the experiment. This study was conducted in accordance with the Helsinki declaration. It was reviewed and approved by the institutional Human Research Ethics Committee of Universiti Sains Malaysia (FWA Reg No: 00007718; IRB Reg No: 00004494). In addition, we followed the recommendations laid out by the guidelines for conducting and reporting MEG research (Gross et al., 2013; Keil et al., 2014). The sample size for this study was estimated using G* power and satisfied statistical power of 95% at a 0.05 two-sided significance level. The sample size was estimated based on the data observed in a previous study (Wu et al., 2012). The estimated total sample size for this study was 26. To mitigate potential attrition, we eventually decided to include 30 subjects in the study.

In this study, 32 subjects underwent screening. However, two individuals were excluded due to psychiatric illnesses, specifically, schizophrenia and generalized anxiety disorder. Thus, our final study cohort comprised 30 enrolled subjects (aged 21–35 years old). 15 of them were male and 15 were female. All of them were non-Arabic speakers. They have no known neurological disease, history of head trauma, hearing problem or psychiatric illnesses. None have history of illicit drugs abuse. All subjects provided written informed consent prior to the study. Each subject was advised to have an adequate sleep of 8 h every night for 1 week prior to the study. They were also advised not to consume any caffeinated foods and beverages, alcohols, and other psychoactive substances. Subjects with any metal containing implants were excluded from this study. Flow chart in Fig. 1 summarise the whole process of subject recruitment and MEG acquisition.

Stimuli

Five different auditory stimuli were considered in this study. The session starts with 3 min recording of resting state where the subjects were asked to sit comfortably with their eyes closed without being presented with any external stimuli or engaging in any activity. A specified order of five different auditory stimuli were subsequently administered for the total length of 3 min per stimuli. The stimuli were administered with Presentation software (Neurobehavioral Systems, Ltd., Berkeley, CA, USA). The sound was delivered through a pair of pneumatic headphones at individually adjusted loudness. In between stimuli, there was a 1-minute rest period where no auditory stimuli were given. The session ended with 3-minutes of resting state recording.

Both rhythmic and non-rhythmic auditory stimuli were considered in this study. Table 1 lists all the stimuli used in this study and their characteristics. Each stimulus contributes to the overall understanding of the neural mechanisms of music therapy. Monochord is an instrument having all of its 30 strings tuned to one base tone, producing varying overtones sounds that merge into one continuous sound without any specific scale and harmony. The produced sounds have been shown to have relaxation effects in children with anxiety (Goldbeck & Ellerkamp, 2012). Similar findings was also noted in patients who underwent chemotherapy (Lee et al., 2012). The monochord’s single-tone sound could help in understanding how simple, repetitive sounds affect the brain’s emotion network. Hare Krishna is one of the mantra recited in Hindu religion as part of meditation practice. Meditation has been shown to be an effective mean of emotional control (Braboszcz, Hahusseau & Delorme, 2010). Its inclusion helps in understanding how a melodious religious chant affects the neural network regulating emotion. Al-Kursi is the 225th verse of the 2nd chapter in the Quran. The verse was recited by a renown reciter and was recorded in a dedicated sound-proof Audio Room, School of Medical Sciences, USM Health Campus. This verse is used in Ruqyah which is a method of treatment using Quranic verses and supplications practised within the the muslim community (Haque & Keshavarzi, 2014; Abu-Rabia, 2005). Its inclusion in the study can shed light on how melodious recitation might influence emotion processing brain network informing the neural mechanism of its therapeutic effects.

Apart from the obvious reason of being composed in the same language, Arabic news and poems were included as the control of the Al-Kursi stimuli since their sound originated from the same instrument as the Quranic recitations, i.e., human vocal cord. Arabic poem titled “If my Lord asks me” was chosen for its rhythmic features which is comparable to the rhythmic features of the Quranic recitation, while Arabic news was chosen since in contrast to Quranic recitation, it is non-rhythmic in nature. By analyzing the different responses these stimuli evoke in the brain’s emotion network, the study aims to explore the neural underpinnings of music therapy. This could potentially lead to more targeted and personalized therapeutic interventions. The audio files of all stimuli used in the present study were available online (see Data Availability section).

Acquisition and preprocessing of MEG datasets

Simultaneous EEG and MEG data were recorded at the MEG laboratory of Universiti Sains Malaysia. The recording were conducted in an electrically and magnetically shielded room with Vectorview™ 306-channel MEG system (Elekta Neuromag; Elekta Oy, Helsinki, Finland) combined with a compatible 64-channels EEG system. The MEG system consist of 204 gradiometers and 102 magnetometers. Subjects were instructed to remove any metals, electronics or metal containing clothing before entering the room. The recording sample rate was set at 1000 Hz. The reference and ground electrode were put on the nose tip and on the right cheek respectively. Eye movements and blinks, were recorded with 2 electrooculogram (EOG) electrodes attached close to the external eye corners on both sides. Electrocardiogram (ECG) was also simultaneously recorded for the purpose of removing artifacts of cardiac origin. Four head position indicator (HPI) coils were embedded within the head cap and used to record head movement and subsequently used for movement artifact correction. These locations were determined relative to the nasion, two preauricular points and more than 100 additional points on the surface of scalp and nose which were recorded by the Polhemus Isotrak 3D digitizer. Through out experiment, subjects were comfortably seated with eyes closed. Overall, we have the recordings of 30 subjects for each stimuli except the 8th stimuli and post resting state where there were recordings from only 27 subjects.

MEG data was preprocessed using Brainstorm software (Tadel et al., 2011) and MNE software (Gramfort et al., 2014). A series of preprocessing steps were done in order to suppress or remove noise and artifacts. Bad channels were identified and excluded from analysis. We then applied notch filter to remove power line noise at 50 Hz, 100 Hz, and 150 Hz peaks. The temporal signal space separation (tSSS) were used to suppress noise originated from outside of the MEG dewar and to correct for head motion artifacts (Taulu & Simola, 2006). The following setting parameters were used: length of correlation window, 30 s; subspace correlation limit, 0.98; order of internal component of spherical expansion, 8; and order of external component of spherical expansion, 3. Physiological artifacts of cardiac origin, eye movements, eye blinks and muscular activity were isolated and removed using independent component analysis (ICA) method (Barbati et al., 2004). The resulting MEG were manually inspected to ensure satisfactory artifacts correction and also to identify any remaining or additional bad segments which were then excluded from subsequent analysis.

Source localization

Source localization was done using Brainstorm software (Tadel et al., 2011) by following the pipeline described in Niso et al. (2019). Specifically, we used fiducial and digitization points of each subjects to warp the template MRI of International Consortium for Brain Mapping (ICBM152) (Fonov et al., 2009) in order to create a psudo-individual anatomy. The warped MRIs were then coregistered with the MEG sensors for each recording sessions. We check the result to ensure satisfactory alignment and orientation of the MEG sensors and the warped MRI.

We modelled the remaining sensors and environmental noise as noise covariance matrix using 2 min of empty room MEG recordings. We then defined forward model for each recording sessions using the overlapping spheres method (Huang, Mosher & Leahy, 1999). The forward or head models describe how neuronal activity at the cortical sources is being transformed to the magnetic flux measured at the level of MEG sensors. In order to capture neuronal activity at cortical surfaces as well as at deep subcortical structures, we defined source space as a set of distributed dipoles of the whole brain volume with unconstrained dipoles orientation.

Once forward or head models were defined for each recording sessions, source estimation was carried out for each of them using depth weigthed minimum norm estimation (wMNE) (Lin et al., 2006). This method has been validated in many studies using either simulation or by comparing the estimation with intracranial EEG recording (Halder et al., 2018; Afnan et al., 2023; Mikulan et al., 2020; Pascarella et al., 2023). The depth weighting parameter was set at 0.75 as suggested by Lin et al. (2006) in order to minimize the localization error particularly for the deeper subcortical structures. The resulting inverse operators were then used to reconstruct time series at the level of brain regions of interest. We used automated anatomical labelling atlas 3 (AAL3) (Rolls et al., 2020) to segment volume source space and define brain regions of interest. The time series for each brain regions were computed as the mean of all voxels within each regions.

Brain regions of interest

We considered 61 brain regions that have been shown in the literature to be involved in emotional regulation (Banks et al., 2007; Cisler et al., 2013; Kohn et al., 2014). Among these are brain regions that were stipulated to be involved in experiencing musical emotions (Alluri et al., 2015). Dysfunction of these regions have also been shown to be associated with several neuropsychiatric pathologies including depression (Greicius et al., 2007; Veer et al., 2010; Lui et al., 2011), obsessive compulsive disorder (Harrison et al., 2009; Sakai et al., 2011; Hou et al., 2014; Takagi et al., 2017), and anxiety (Etkin et al., 2009; Kim et al., 2011; Hahn et al., 2011). We have also included brain regions that have been used in functional neuroimaging studies to monitor and assess treatment responses in these disorders (Crowther et al., 2015; Walsh et al., 2017; Scult et al., 2019; Walsh et al., 2019). Table 2 lists all the brain regions that were included in the analysis. These brain regions were specified by selecting and merging together brain regions defined in AAL3.

Connectivity analysis

For the purpose of our study, we applied amplitude envelope correlation (AEC) to extract bivariate connectivity strength between the pairs of all to all 61 × 61 brain regions (Bruns et al., 2000a; Bruns et al., 2000b). Connectivity strength were computed at the level of brain regions using time series estimated for each brain regions as described in source localization section.

The time series were first band pass filtered to frequency bands of interest: delta (δ) (2 − 4 Hz), theta (θ) (5 − 7 Hz), alpha (α) (8 − 12 Hz), beta (β) (15 − 29 Hz), and gamma (γ) (30 − 59 Hz). This frequency ranges were the convention used in Brainstorm software and have been used in many prior studies (Gehrig et al., 2012; Wiesman et al., 2022; Rempe et al., 2023; Bergwell et al., 2023; Nugent et al., 2020). In order to avoid spurious connectivity between the reconstructed source time series caused by signal leakage, we applied orthogonalized AEC where time series were orthogonalized prior to computation of AEC (Hipp et al., 2012). Orthogonalization of time series was performed in frequency domain before computing their power envelopes. A complex signal y_⊥x(t, f) is orthogonal to signal x(t, f) as defined by the following equation: (1)${\displaystyle y_{\perp x}\left(t,f\right)=\Im \left(\frac{x{\left(t,f\right)}^{\ast }}{\left|x\left(t,f\right)\right|}y\left(t,f\right)\right)}$

where ℑ denotes imaginary part of the complex signals and x(t, f)^∗ is the complex conjugate of x(t, f).

Analytical signal z(t) were then constructed using Hilbert transform,

(2)${\displaystyle z\left(t\right)=x\left(t\right)+iH\left[x\left(t\right)\right]}$ (3)${\displaystyle H\left[x\left(t\right)\right]=\frac{1}{\pi }PV{\int }_x^{\infty }\frac{x\left(\tau \right)}{t-\tau }d\tau }$

where i is the imaginary part, H[x(t)] is the Hilbert transform and PV is the Cauchy principle value. The amplitude envelope were then computed as follows: (4)${\displaystyle a\left(t\right)=\sqrt[]{{\left(x\left(t\right)\right)}^2+{\left(H\left[x\left(t\right)\right]\right)}^2}}$

With amplitude envelope of source time series obtained, connectivity (AEC) between amplitude envelope of ith brain region a_i and jth brain region a_j was computed as a Pearson correlation coefficient ρ between them: (5)${\displaystyle \rho \left(\mathbf{}{a}_{\mathbf{}i},\mathbf{}{a}_{\mathbf{}j}\right)=\frac{\mathbf{}{a}_{\mathbf{}i}{\mathbf{}{a}_{\mathbf{}j}}^{\mathbf{}T}}{\sqrt{\mathbf{}{a}_{\mathbf{}i}{\mathbf{}{a}_{\mathbf{}i}}^{\mathbf{}T}}\sqrt{\mathbf{}{a}_{\mathbf{}j}{\mathbf{}{a}_{\mathbf{}j}}^{\mathbf{}T}}}}$

where a_i and a_j are 1 × T vectors of mean centered amplitude envelope of ith and jth brain regions respectively. a_i^T and a_j^T are both transposes of a_i and a_j. Correlation coefficient of all to all brain regions were organized into matrix forming connectivity or weighted adjacency matrix A.

Network analysis

Once connectivity matrices A ∈ ℝ^61×61 with elements w_ij representing connectivity strength (AEC) between ith and jth brain regions were computed for each subjects, stimuli and frequency, we summarized the results using tools developed in graph theory. Graph theory offers several tools to characterize the topological features of complex brain networks. In order to ensure robust findings, the connectivity matrices were thresholded leaving out small and insignificant connectivity values which can be considered as noises.

We adapted the recommendations proposed by Cohen (2014) to threshold connectivity matrices using absolute thresholding method with median value of connectivity strength as the cut-off point. In order to ensure the threshold is independent of each stimuli, we determined the threshold value based on pooled values of connectivity strength from all stimuli. We determined specific threshold separately for each frequency bands since connectivity values differ across frequency bands.

In this study, we computed two topological features from the thresholded connectivity matrices using the Brain Connectivity Toolbox (BCT) (Rubinov & Sporns, 2010). These are weighted transitivity and weighted global efficiency. Weighted transitivity and global efficiency are particularly useful and relevant measures in this study since they allow for the assessment of the efficiency of both local information processing and global information transfer respectively.

Weighted transitivity is a variant of mean clustering coefficient which reflects the average prevalence of clustered connectivity around individual nodes. Weighted transitivity T is measured using the following formula: (6)${\displaystyle T=\frac{\sum _i^{N}2t_i}{\sum _i^N{k}_i\left(k_i-1\right)}}$

where k_i is the number of links connected to node i. k_i(k_i − 1) is thus the total number of possible links that could exist among the vertices. N is the total number of nodes in the network and t_i is the weighted geometric mean of triangle around node i computed as follows: (7)${\displaystyle t_i=\frac{1}{2}\sum _{j,h}^N{\left(w_{ij}{w}_{ih}{w}_{jh}\right)}^{\frac{1}{3}}}$

where w_ij is the connectivity strength from node i to node j. All self connection strength w_ii were set to zero (w_ii = 0 for all i).

Weighted global efficiency E measures the efficacy of the brain network in integrating information from distributed brain regions and is calculated using the following formula: (8)${\displaystyle E=\frac{1}{N}\sum _i^N\frac{\sum _{j,j\ne i}^N{\left(d_{ij}\right)}^{-1}}{N-1}}$

where d_ij is the shortest path length between nodes i and j calculated as follows: (9)${\displaystyle d_{ij}=\sum _{w_{uv}\in g_{i\to j}}f\left(w_{uv}\right)}$

where f is a map from connectivity weight to length and g_i→j is the shortest weighted path between node i and node j.

Statistical analysis

The analysis were done using R statistical programming language (R Core Team, 2020). We analysed the data using Bayesian linear mixed model. Brms package were used to fit the model (Bürkner, 2017). Brms is the R wrapper for Stan which is another probabilistic programming language that implement Hamiltonian Monte Carlo (HMC) sampling algorithm (Carpenter et al., 2017). Following the recommendation by Barr et al. (2013), the model include the maximal random effects structure justified by the experimental design, i.e., varying intercept and slope model with subjects as the grouping variable. Outliers were identified using Boxplot method where all data points that fall outside the threshold of 1.5 times the interquartile range below the first quartile or above the third quartile and were excluded before model fitting. As part of sensitivity analysis, we have also performed the statistical analysis while including the outliers. We have also conducted sensitivity analysis on strategies for handling missing data by using both only complete cases and multiple imputations to account for missing data. In line with our research objectives and experimental design, we fitted the following model for each network measures:

(10)${\displaystyle y_{i\mid s}\sim \mathit{Log}\mathcal{N}\left({\mu }_s,\sigma \right)}$ (11)${\displaystyle {\mu }_s={\alpha }_s+\sum _j{\beta }_{s,j}Stimul{i}_j}$ (12)${\displaystyle \left[\begin{array}{c}\hfill {\alpha }_s\hfill \\ \hfill {\beta }_{s,j}\hfill \end{array}\right]\sim \mathsf{MVNormal}\left(\left[\begin{array}{c}\hfill \bar{\alpha }\hfill \\ \hfill \bar{\beta }\hfill \end{array}\right],\mathbf{Q}\right)}$ (13)${\displaystyle \mathbf{Q}=\left[\begin{array}{cc}\hfill {\sigma }_{\alpha }\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\sigma }_{\beta }\hfill \end{array}\right]\mathbf{R}\left[\begin{array}{cc}\hfill {\sigma }_{\alpha }\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\sigma }_{\beta }\hfill \end{array}\right]}$ (14)${\displaystyle \bar{\alpha }\sim \mathcal{N}\left(\overline{{\mu }_k},0.2\right)}$ (15)${\displaystyle \bar{\beta }\sim \mathcal{N}\left(0,0.2\right)}$ (16)${\displaystyle {\sigma }_{\alpha }\sim \mathit{half}\mathcal{N}\left(0,0.1\right)}$ (17)${\displaystyle {\sigma }_{\beta }\sim \mathit{half}\mathcal{N}\left(0,0.1\right)}$ (18)${\displaystyle \mathbf{R}\sim \mathrm{LKJcorr(2)}}$ (19)${\displaystyle \sigma \sim \mathit{half}\mathcal{N}\left(0,0.1\right)}$

where y_i∣s refers to ith observation of dependent variables of interest (network measures - T and E) of subject s. $\mathit{Log}\mathcal{N}\left(\mu ,\sigma \right)$ is the notation for log-normal distribution with location µand scale σ. α_s is the specific intercept for subject s (pre-stimuli resting state). β_s,j is the coefficient parameter for jth stimulus (how much the stimulus differs from the pre-stimuli resting state recording; here the stimuli include post-stimuli resting state recordings) in subject s. Both α_s and β_s,j are jointly distributed with multivariate normal distribution around the overall intercept $\bar{\alpha }$ and slope $\bar{\beta }$. Q ∈ ℝ^2×2 is the covariance matrix of the multivariate Gaussian distribution defined as in (Eq.13). σ_α and σ_β are standard deviations of α_s and β_s,j respectively. R ∈ ℝ^2×2 is the correlation matrix. Both Eqs. (10) and (11) define the likelihood function and the rest (Eqs.12–19) define priors of our generative model. Additionally, as a component of our sensitivity analysis, we have also utilized the Student’s t-distribution to model the data. We have also tested the sensitivity of using non-informative priors where $\bar{\alpha }\sim \mathcal{N}\left(\overline{{\mu }_k},10\right),\bar{\beta }\sim \mathcal{N}\left(0,100\right),{\sigma }_{\alpha }\sim \mathit{half}\mathcal{N}\left(0,10\right),{\sigma }_{\beta }\sim \mathit{half}\mathcal{N}\left(0,10\right),$ and $\sigma \sim \mathit{half}\mathcal{N}\left(0,10\right)$.

We set weakly informative priors (McElreath, 2020) on the parameters $\bar{\beta }$ i.e the overall slope of the stimuli, the standard deviation of the intercept σ_α, the standard deviation of the slope σ_β, the correlation matrix R and the standard deviation of the models σ. For the purpose of our research, we set the prior for $\bar{\beta }$ so that the slopes values are symmetric around 0 with finite variance. We set priors for parameter $\bar{\alpha }$ to be normally distributed around each network measure mean $\overline{{\mu }_k}$ with finite variance. We checked the prior predictive distributions to ensure the chosen priors make sense and generate plausible and realistic data prior to observations.

Four sampling chains were used to explore and sample from the posterior distributions for each model. These chains ran for 14,000 iterations with 4,000 warm-up iterations yielding 40,000 samples for each parameters. Several HMC diagnostics were monitored to determine the convergence of numerical integration of posterior distribution. These include the potential scale reduction factor $\hat{R}$, traceplots and effective sample size (ESS). $\hat{R}$ of <1.1 indicates convergence of the sampling procedure and EES should be >10, 000 samples for stable posterior estimates (Kruschke, 2014). Once we fit the model, we checked the posterior predictive distributions to ensure that the model reflects the observed data.

We set equivalence regions based on principles that were laid out in Kruschke (2018); Kruschke & Liddell (2018); Kruschke (2011). We proceed by specifying the ROPE which is the range of values within which the difference between two conditions or stimuli is taken to be small enough to be practically equivalent. Specifically, we used the total credible interval (CI) of the difference between pre-stimuli resting state and post-stimuli resting state recordings to specify the ROPE. This ensures that the equivalence region covers the whole interval of the differences between the two resting states. The 90% CIs for each stimulus and their differences were defined using the 90% highest density intervals (HDI). These HDI were then used together with ROPE to test for the presence or absence of effects. The network measures for each stimulus were considered equivalence to the resting state recording (absence of effects) if the entire 90% CI of the difference (i.e., stimuli—pre-resting) lie within the equivalence region. If the entire 90% CI of the difference was outside the ROPE, the null hypothesis of equivalence or absence of effects should be rejected in favour of the alternate hypothesis of difference or presence of effects. If the 90% CI was neither completely lie within nor completely lie outside the ROPE, the presence or absence of effects was undecidable. We reported the means and 90% CI under posterior distributions for the differences between each parameters.