Effects of an Oral Contraceptive on Dynamic Brain States and Network Modularity in a Serial Single-Subject Study

Hormonal contraceptive drugs are used by adolescent and adult women worldwide. Increasing evidence from human neuroimaging research indicates that oral contraceptives can alter regional functional brain connectivity and brain chemistry. However, questions remain regarding static whole-brain and dynamic network-wise functional connectivity changes. A healthy woman (23 years old) was scanned every day over 30 consecutive days during a naturally occurring menstrual cycle and again a year later while using a combined hormonal contraceptive. Here we calculated graph theory-derived, whole-brain, network-level measures (modularity and system segregation) and global brain connectivity (characteristic path length) as well as dynamic functional brain connectivity using Leading Eigenvector Dynamic Analysis and diametrical clustering. These metrics were calculated for each scan session during the serial sampling periods to compare metrics between the subject’s natural and contraceptive cycles. Modularity, system segregation, and characteristic path length were statistically significantly higher across the natural compared to contraceptive cycle scans. We also observed a shift in the prevalence of two discrete brain states when using the contraceptive. Our results suggest a more network-structured brain connectivity architecture during the natural cycle, whereas oral contraceptive use is associated with a generally increased connectivity structure evidenced by lower characteristic path length. The results of this repeated, single-subject analysis allude to the possible effects of oral contraceptives on brain-wide connectivity, which should be evaluated in a cohort to resolve the extent to which these effects generalize across the population and the possible impact of a year-long period between conditions.


Supplementary Figures
Supplementary Figure S1: Graph theory static connectivity and Leading Eigenvector Dynamics Analysis (LEiDA) and diametrical clustering measurement of dynamic connectivity (A) System segregation and modularity describe the functional subdivision of the brain into discrete subnetworks, while characteristic path length describes the capacity for information flow across nodes (brain areas). In weighted networks edges are represented as continuous values, each one defining the edge strength; whereas in binarised networks, edges below a threshold are pruned, retaining only edges above a threshold. In LEiDA, (B) the Hilbert transform is used to separate amplitude a(t) and phase ( ) information from each regional rs-fMRI time series, for each scan session. (C) For each time-point, t, a coherence map is estimated and decomposed using the eigenvalue decomposition, retaining the leading eigenvector for further analysis. (D) The set of leading eigenvectors within and across scan sessions is clustered into k predefined clusters using diametrical clustering, e.g., three clusters c1-c3 shown here. (E) The original rs-fMRI time-series is partitioned according to the estimated brain states; fractional occurrence, i.e., the fraction of time points during which the brain is in a given brain state, can be calculated. For LEiDA, we estimated the autocorrelation across all 209 estimated brain states (for ∈ {2, … ,20}) and report the average and standard deviation of the coefficients (E). Figure S3: Modularity and bCPL estimates across a range of sparsity thresholds.

Supplementary
Both modularity and bCPL are sensitive to the threshold r2z-score as a graph pruning parameter.
Here, paired t-test estimates for modularity (A) and corresponding p-values (B) indicate a statistically significant difference between OC and naturally cycling scans across the range of thresholds. Note that the significance threshold of 0.05 is not shown in figures B and D, as this would make the points appear as straight lines.
Supplementary Figure S4: Highlighted centroids across k including structural similarity (A) and between-k Pearson correlation (B).

fMRI Preprocessing
Each of the acquired 10-minute rs-fMRI scan sessions were preprocessed separately in SPM12  (Yeo et al., 2011) were used to extract denoised regional time series for further analysis.

Hilbert transform
Regional instantaneous BOLD-phases were estimated using the Hilbert transform where * represents the convolution operator, to create the complex analytic signal ( ) = ( ) + ℎ ( ). The analytic signal encompasses both signal amplitude and oscillatory information by projecting it onto counterclockwise motion, phase-shifted by 4 . The instantaneous phases may then be found as the point-wise angle to the real axis ( ) = arctan ( ℎ ( ) ( ) ), and represent the oscillatory dynamics of the BOLD signal bounded by ∈ [− ; ). Similarly, the instantaneous amplitude of the BOLD signal may be found as the modulus of the analytic signal, however, in the attempt to uncover phase coherence networks, amplitude information is disregarded here.

System Segregation
System segregation (SS) measures the relative strength of within-network connectivity compared to between-network connectivity (Chan et al., 2014), and is defined as where ̅ (within-network) represents the average r2z score for all edges that connect two nodes from the same network, and ̅ (between-network) represents the average r2z score for all edges that connect two nodes from distinct networks. Conceptually, higher system segregation implies greater within-relative to between-network connectivity strength.

Characteristic Path Length
Characteristic Path Length (CPL), is the average shortest path length between all pairs of nodes; it is sometimes referred to as "average path length" and is the inverse of "global efficiency" (Rubinov and Sporns, 2009). CPL is a global measure of the functional integration of a system (Watts and Strogatz, 1998). Low values indicate short paths between nodes, putatively representing a highly integrated system. CPL can be described using either a weighted or binarised connectivity matrix (wCPL or bCPL, respectively). The CPL for node i is defined as: where ̅ represents the average path length from node i to all other nodes j (n = 400, being the total number of nodes) and represents the shortest path length between nodes i and j. CPL is then the average across n nodes. The shortest path between any given pair of nodes is the path that minimises the total length between a given pair of nodes, which may be direct from node i to node j, or via other nodes, in which case the path length is the sum of the lengths between each node. For wCPL individual edge lengths are defined as 1-r where r is the connectivity strength expressed as the Pearson's correlation coefficient (Rubinov and Sporns, 2009). bCPL was calculated as wCPL except using the binarised connectivity matrix for each scan session; individual edge lengths between nodes were either 1 for super-threshold edges or infinite for sub-threshold edges.
The graph-theoretical interpretation of negatively weighted connectivity is not well established (Rubinov and Sporns, 2011) and thus connectivity values are set to zero as recommended in (Rubinov and Sporns, 2011), resulting in a path length of 1.

Modularity
Modularity (Q) is a measure of the number of suprathreshold, within-network connections relative to the number of suprathreshold, between-network connections (Cohen and D'Esposito, 2016;Sporns and Betzel, 2016) and is defined as where represents the fraction of within-network connections that are suprathreshold for network i, represents the fraction of between-network connections where one node belongs to network i that are suprathreshold, while m represents the total number of networks (i.e., seven). As such, modularity is a measure of the degree to which a system is segregated into functionally distinct networks as opposed to homogenous connectivity across the system. Higher values indicate greater segregation into distinct networks.

Leading Eigenvector Dynamics Analysis
We estimate dynamic connectivity structures using LEiDA (Cabral et al., 2017), followed by diametrical clustering as in (Olsen et al., 2021)(see Figure S1B-E). The brain was parcellated into cortical regions from the Schaefer-100 atlas (Schaefer et al., 2018). Regional instantaneous BOLDphases were estimated using the Hilbert transform (see section 2.3), which separates instantaneous phases representing the oscillatory dynamics of the BOLD signal from its amplitude.
Interregional synchrony is summarised for each time point , for each scan session, in the diagonally symmetrical phase coherence map with elements , , = cos( , − , ) for regions and . As per LEiDA, the first eigenvector, 1, , following an eigenvalue decomposition of the brain-wide coherence map is retained, thereby capturing the dominant instantaneous connectivity pattern.

Diametrical clustering
Eigenvectors are ambiguous in both their norm and sign and are therefore well-clustered with respect to an antipodally symmetric unit hypersphere (Watson, 1965;Sra and Karp, 2013;Olsen et al., 2021).
Diametrical clustering serves this purpose, (Dhillon et al., 2003;Sra and Karp, 2013) and is essentially a modified k-means procedure, where distance , is measured for every time point t using the squared Pearson correlation, e.g., where is the centroid of cluster . We group all leading eigenvectors concatenated across scan sessions into k clusters (i.e., brain states). Diametrical clustering is initialised using an implementation of k-means++ (Arthur and Vassilvitskii, 2007). The optimal k is not well-defined; therefore, we ran our model for k ranging from 2 to 20, in line with previous studies examining brain dynamics using LEiDA (Cabral et al., 2017;Figueroa et al., 2019;Kringelbach et al., 2020;Olsen et al., 2021). For each k, the best of 5 trained models, in terms of average distance of points to the nearest centroid, was selected. We define brain state fractional occurrence as the fraction of volumes within a single rs-fMRI scan session assigned to a specific brain state. This produces = 60 fractional occurrence estimates per brain state.

Atlas considerations
In this study, we investigated graph-theoretical measures of brain activity and dynamic functional connectivity using the Schaefer atlases parcellated into 400 and 100 regions, respectively (Schaefer et al., 2018). The original article suggested the use of either the 400 or the 1000-region parcellations for the analysis of functional data. However, clustering methods, including diametrical clustering, suffer from the "curse of dimensionality", hence a lower number of regions is preferable. We have previously analysed a similar-sized dataset using a 90-dimensional atlas and observed reasonable convergence properties (Olsen et al., 2021). Regions in the Schaefer atlas are labelled according to a canonical functional network as defined in (Yeo et al., 2011), thereby enabling the interpretation of dFC-results in terms of network implication. As such, the Schaefer atlas is preferable over anatomical atlases where regions are not easily assigned to groups. In contrast, region-specific alterations are easily overlooked when working with network atlases, and similarly, in this study, we did not focus on anatomically localised subnetworks.