Changes in resting-state functional connectivity in neuropsychiatric lupus: A dynamic approach based on recurrence quantification analysis

There is growing interest in dynamic approaches to functional brain connectivity (FC), and their potential applications in understanding atypical brain function. In this study, we assess the relative sensitivity of cross recurrence quantification analysis CRQA) to identify aberrant FC in patients with neuropsychiatric systemic lupus erythematosus (NPSLE) in comparison with conventional static and dynamic bivariate FC measures, as well as univariate (nodal) RQA. This technique was applied to resting-state fMRI data obtained from 45 NPSLE patients and 35 healthy volunteers (HC). Cross recurrence plots were computed for all pairs of 16 frontoparietal brain regions known to be critically involved in visuomotor control and suspected to show hemodynamic disturbance in NPSLE. Multivariate group comparisons revealed that the combination of six CRQA measures differentiated the two groups with large effect sizes (.214 > η2 > .061) in 40 out of the 120 region pairs. The majority of brain regions forming these pairs also showed group differences on nodal RQA indices (.146 > η2 > .09) Overall, larger values were found in NPSLE patients vs. HC with the exception of FC formed by the paracentral lobule. Determinism within five pairs of right-hemisphere sensorimotor regions (paracentral lobule, primary somatosensory, primary motor, and supplementary motor areas), correlated positively with visuomotor performance among NPSLE patients (p < .001). By comparison, group differences on static FC displayed large effect sizes in only 4 of the 120 region pairs (.126 > η2 > .061), none of which correlated significantly with visuomotor performance. Indices derived from dynamic, temporal-based FC analyses displayed large effect sizes in 11/120 region pairs (.11 > η2 > .063). These findings further support the importance of featurebased dFC in advancing current knowledge on correlates of cognitive dysfunction in a clinically challenging disorder, such as NPSLE.


I. INTRODUCTION
It is widely acknowledged that brain function is highly dynamic. This property greatly complicates efforts to assess functional associations in the neurophysiological activity displayed concurrently by two or more brain regions, which is essential in generating all psychological phenomena. Such associations are considered as evidence of functional connectivity (FC) between regions. In studies utilizing functional Magnetic Resonance Imaging (fMRI) to measure temporal variations in hemodynamic activity, FC is more commonly estimated through static measures, such as correlation coefficients between two time series, at a time. This approach, however, does not take into account a well-known characteristic of brain signals, namely non-stationarity. Moreover, it is not designed to capture deterministic behavior in the form of high degrees of predictability of the hemodynamic signal recorded in one brain region from the signal recorded in other regions. Accordingly, there has been a growing recognition of the need to model brain activity through dynamic approaches [1].
Dynamic Functional Connectivity (dFC). The fluctuation of temporal dynamics is typically estimated using a slidingwindow approach. Analogous to static FC, Pearson correlations are computed between pairs of regional time series for shorter, overlapping portions of the time series. Many dynamic functional connectivity approaches are based on the study of several network states and their fluctuation throughout the duration of the scan. For instance, [2]- [4] utilize an unsupervised machine learning method (k-means clustering) for the identification of repeating connectivity patterns. A related approach, where functionally more well-defined temporal states are identified and examined, can be found in [5]- [7]. The connectome of each time-block (∼53 sec) is characterized as predominately integrated or segregated.
Functional modules are groups of nodes that exhibit high connectivity with the module's nodes and weak or negative with nodes from other modules. Integrated networks exhibit strong global connections mediated through hubs, whilst segregated networks have more pronounced strongly connected communities that are weakly connected with each other. The state identification process is achieved by evaluating a combination of network characteristics via their joint histogram, namely participation coefficient and within-module degree zscore. These measures quantify each node's inter-and intramodular connectivity, accordingly. This version of node degree measures the connections of a particular region with other regions belonging to the same module. On the other hand, the participation coefficient quantifies a node's connectivity to nodes of other modules. For this reason, networks with high mean participation coefficient are often characterized as integrated.
Feature-based dFC. Although brain signals are not stationary, they may display the phenomenon of recurrence, which is typical in dynamic systems. A fundamental property of a recurrent system is that similar situations often evolve in a similar manner. Recently, recurrent patterns of task-related hemodynamic activity have been modeled through recurrence quantification analysis (RQA), which was shown to be quite robust to variations in hemodynamic responses across regions, tasks and participants [8]. RQA relies on the construction of recurrence plots (RP) and computes measures of the degree to which activity in a given brain region displays similar patterns (feature states) over time. This method can be extended to allow quantification of recurrent signal patterns displayed by two concurrent signals, by computing cross recurrence plots (CRP) [9].
The basic idea of CRPs is to study the similarity of two different processes, whose phase space trajectories belong in the same phase space. The CRP was introduced to exploit and analyze the dependencies and interrelations between two different systems, by comparing their states. Suppose that we aim to reconstruct a one-dimensional time series which comes from a higher-dimensional dynamics, of dimension d. The main goal is to start embedding time-delayed copies of the one-dimensional time series, of delay τ , thus creating a phase space which theoretically can reach infinity. This phase space's dimension is based on the embedding dimension parameter m, which essentially determines how many time-delayed copies we need so as to reconstruct the original time series dynamics. For each new necessary copy, we examine its location in the phase space in relation to the previous one. Accordingly, the m parameter cannot exceed the value of 2d + 1 and helps the embedding procedure to converge to the minimum possible dimension of the phase space, in which all the copies can exist. For two time series the above procedure has to be repeated.
Recurrent signaling may be even more prominent when the brain is at rest, that is, without the constraints imposed by specific (cognitive or motor) task demands. Measurements of such activity can be obtained using EEG signals (ref. [10] for a recent application of cross-recurrence analysis) and resting-state fMRI (rs-fMRI), which purportedly demonstrate self-organizing (i.e., at least partially deterministic) behavior. In that case, systematic alterations in the degree of signal recurrence may serve as a more sensitive index of aberrant brain function.
A disorder that is characterized by a variety of neuropsychiatric symptoms in the absence of remarkable lesions in the brain is neuropsychiatric systemic lupus erythematosus (NPSLE). Recently, our group has identified regional hubs that optimally differentiate NPSLE patients from healthy controls based on graph theoretic metrics derived from static connectivity indices [11]. We have also found evidence of impaired perfusion in the normal-appearing white matter of the semioval center in NPSLE patients utilizing dynamic susceptibilitycontrast-enhanced T2 * -weighted perfusion MRI [12]. Major fiber bundles transverse this area, including the superior longitudinal fasciculus, which interconnects superior frontal and parietal regions in both hemispheres. Given that these regions play an important role within the circuit that supports visuomotor function [13]- [15], it is likely that individual differences on functional connectivity among NPSLE patients would be significantly associated with visuomotor task performance.
In this work, we examine the suitability of cross recurrence quantification analysis (CRQA) to capture recurrent dynamic characteristics of rs-fMRI time series, which could serve as complementary indices of functional connectivity. Specifically, we assess the relative sensitivity of this method in identifying aberrant FC in a group of NPSLE patients as compared to a gender and age-matched group of healthy volunteers. In view of the potential vulnerability of brain regions normally involved in visuomotor coordination and control to NPSLE-related pathology, we adopted a region of interest approach focusing on aberrant dynamic connectivity between 16 frontoparietal regions (eight in each hemisphere). The majority of the selected regions are considered as parts of the sensorimotor resting-state network (primary sensory and motor cortex, dorsal premotor and supplementary motor cortices, superior parietal and adjacent paracentral lobule) and the Default Mode Network (angular gyrus, precuneus; ref. [16]). We further hypothesized that CRQA-based features, for those pairs of regions that differentiate the two groups, will significantly correlate with performance on a task of visuomotor coordination and control (Grooved Pegboard task). The sensitivity of CRQA metrics of FC is evaluated against a conventional static FC approach, estimates of nodal RQA, as well as with metrics of dFC derived from a previously established approach (ref. [5], [7], [17]).

A. Participants
The main group that provided fMRI and cognitive assessment data for the present analyses comprised 45 patients diagnosed with NPSLE, recruited from the Rheumatology outpatient clinic of the University Hospital of Heraklion by their attending physician. All patients met the revised ACR 1997 classification criteria for SLE [18]. NPSLE diagnosis and attribution was based on physician judgement, following a multidisciplinary approach and considering patient age and risk factors for NPSLE (anti-phospholipid antibodies, prior neuropsychiatric manifestation, generalized disease activity, findings of conventional MRI imaging and other diagnostic procedures). Additional inclusion criteria were: age over 18 years, ability to speak and read Greek fluently, negative history of thromboembolic cardiovascular disease or other primary CNS diseases. Comparison rs-fMRI data were available on 35 age-matched (p > .1) and gender-matched (p > .4) healthy volunteers (see Table I). The hospital review board approved the study and the procedure was thoroughly explained to all patients and volunteers who signed informed consent.

B. Neurocognitive assessment
Visuomotor coordination and efficient fine motor control was assessed using the Grooved Pegboard test (Lafayette Instruments) posing similar demands as the 9 Hole Peg test employed by [14]. The test involves placing 25 3 mm diameter pegs into evenly-spaced keyholes arranged on a 5 × 5 matrix as fast as they could. Performance data using the dominant hand is analyzed here. The time to complete the trial serves as the performance metric. This goal-directed task requires intact motor programming and fine motor control to manipulate the pegs into the correct orientation and efficient visuomotor and proprioceptive integration to place them in their respective keyholes.
C. fMRI data acquisition and preprocessing All brain MRI examinations were performed on a clinical 1.5T whole-body superconducting imaging system (Vision/Sonata, Siemens/Erlangen), equipped with highperformance gradients (Gradient strength: 40 mT/m, Slew rate: 200 mT/m/ms) and a two-element circularly polarized head array coil. The conventional magnetic resonance (MR) imaging protocol consists of a 3D T1-w MPRAGE, a T2wTSE and a Turbo fluid attenuated inversion recovery (FLAIR) sequence. The acquired images were interpreted by a senior neuroradiologist (EP), blinded to the clinical and laboratory data, who reported any incidental findings not related to Systemic Lupus Erythematosus (SLE), or findings due to focal SLE-related abnormalities, such as acute or old infarcts, haemorrhages and focal brain atrophy. Patients with any of the aforementioned MRI findings were excluded from the study.
Rs-fMRI was derived from a T2 * -weighted, fat-saturated 2D-FID-EPI sequence with repetition time (TR) 2300 ms, echo time (TE) 50 ms, and field of view (FOV) 192 × 192 × 108 (x, y, z). Voxel size was equal to 3×3×3 mm, with a whole-brain scan consisting of 36 transverse slices acquired parallel to the plane passing through the anterior and posterior commissures (AC-PC line with 3.0 mm slice thickness and no interslice gap).
Each blood oxygenation level-dependent (BOLD) time series consisted of 150 dynamic volumes (the first five were omitted from the analyses). Preprocessing steps included slicetime correction, realignment, segmentation of structural data, normalization into standard stereotactic Montreal Neurological Institute (MNI) space, and spatial smoothing using a Gaussian kernel of 8 mm full-width at half-maximum using SPM8. Since functional connectivity can be affected by head motion, we performed motion artifact detection and rejection using Artifact Detection Tools 1 . Grey matter, white matter and cerebrospinal fluid (CSF) mean signals were regressed out of all voxel time series, in order to mitigate their effects on fMRI BOLD time courses. The first five principal components of white matter and CSF regions were regressed out of the signal, along with their first order derivatives. Finally, the fMRI time series were detrended and bandpass filtered in the 0.008−0.09 Hz range to eliminate low-frequency drift and high-frequency noise. The above steps were carried via the CompCor tool [19], that is implemented within the CONN preprocessing module and executed in MATLAB.
In the current analysis, we selected a set of 16 frontoparietal brain regions (defined in the automated Anatomical Labeling (aal) atlas [20]), which are known to play an important role in motor control and visuomotor integration in the healthy brain: Primary motor cortex (precentral gyrus, MI), primary somatosensory cortex (postcentral gyrus, SI), supplementary motor cortex (SMA), dorsal aspect of the superior frontal gyrus (SFG, including the dorsal premotor area), superior parietal lobule (SPL) and adjacent paracentral lobule (PCL), and inferior parietal lobule (angular gyrus; ANG), and precuneus.

D. Cross-recurrence quantification analysis
Univariate RQA. The bivariate RQA technique (CRQA) originated from the univariate RQA approach, proposed by Eckmann et al. [21], which relies on recurrence plots (RPs) as a means of visualizing recurrences within a given dynamic system. The objective of univariate RQA is to identify and quantify similarities between pairs of state vectors, which belong to the same time series and correspond to a unique phase space of dimension m. As previously mentioned, the dimension of the phase space depends on the frequency of a measurement's appearance in a time series. This analysis considers a time series of the fMRI data matrix X = [x 1 , . . . , x N ] ∈ R K×N , where x i ∈ R K is the time series acquired from the ith brain region, and N , K denote the number of regions of interest (ROIs) and time points (i.e., samples) per ROI, respectively. Then, through univariate RQA analysis (henceforth referred to as nodal RQA) and having defined the embedding parameter m and time delay τ , we move to the phase space of the trajectory of the state vectors a i ∈ R m , i = 1, . . . , N a , where N a denotes the number of state vectors. A detailed description of the computation of the parameters m and τ is presented in the following paragraphs. The Euclidean distance between two state vectors, located on a trajectory in a phase space, is given by where, the parameter ε determines the criterion of how close two state vectors should be. If two state vectors are located in an area of radius ε then the Heaviside function Θ would return "1". Otherwise, "0". This procedure is illustrated in Fig. 1. The last significant step is to quantify the positive (i.e., "active") dots of the recurrence plot matrix R via various mathematical ways through a number of complementary nodal RQA metrics (see below). Bivariate RQA (CRQA). The objective of the bivariate extension of RQA is to find similarities between pairs of state vectors corresponding to two distinct phase space trajectories with equal embedding dimension m. As before, let the fMRI data matrix be denoted by is the time series acquired from the ith brain region, and N , K denote the number of regions of interest and time points (i.e., samples) per ROI, respectively. In the bivariate RQA extension (CRQA), we examine the trajectories corresponding to pairs of ROIs. Let y i ∈ R m , i = 1, . . . , N y , and z j ∈ R m , j = 1, . . . , N z , denote the associated state vectors of two trajectories, whose number can be, in general, different (N y = N z ). In our implementation N = 16 and K = 145.
The following equation gives the mathematical formulation for generating the key component of CRQA, namely, the cross recurrence plot, where ε is a threshold defining the radius of the area in which two state vectors are considered to be neighboring. If two state vectors are located within the area defined by this radius they are considered as neighbors and the Heaviside function Θ(·) returns 1; otherwise, it returns 0. The first step entails specifying the values of the embedding dimension m and threshold ε, by taking into account the specific characteristics of the recorded time series. We note that, for computational simplicity, the time delay τ is typically set equal to 1, without affecting the sensitivity of the CRQA.
The estimation of the two key embedding parameters, namely, the τ and m, is based on the average mutual information (AMI) criterion and the false nearest neighbor (FNN) rule, respectively [22]. The first step of the (C)RQA process is to reconstruct the dynamics of the original one-dimensional time series via time-delay embedding in a m-dimensional space.
The intuition behind the AMI criterion is to search for a time delay that yields independent time delay coordinates of the original time series against its shifted versions by τ (for various values of τ ), as measured via their mutual information I(x(t), x(t + τ )). In particular, the AMI criterion suggests to set the optimal τ as the first minimum of the AMI versus τ curve. In the CRQA case, distinct time delays, τ y and τ z , are estimated for each time series, whilst the final time delay is given by τ = min{τ y , τ z }. Nevertheless, as mentioned above, in the current implementation we set τ = 1 without affecting the performance of CRQA.
Having estimated τ , the embedding dimension m determines the appropriate number of time-delayed samples for the construction of the phase space. The FNN rule examines, for increasing values of m, whether two state vectors that are neighbors (i.e., their distance is less than a predefined threshold ε) in the current m-dimensional space are still neighbors in the next higher dimensional (m+1) space. Doing so, the optimal m is set equal to the first minimum of the FNN versus m curve. Notice that, in practice, two distinct ROIs are characterized by different embedding dimensions, m y , m z . To address this issue, the final optimal m utilized by CRQA is given by m = min{m y , m z }. Having estimated the optimal embedding dimension m, which determines the dimension of the phase space trajectory, ε is defined as a percentage c of the maximum Euclidean distance within each pair of state vectors. The percentage c is typically set in the range [0.25, 0.35].
The next step quantifies the proximity of two trajectories in phase space, according to a predefined threshold ε, and results in a binary matrix that describes the dependencies between the state vectors of two distinct dynamical systems (represented here by two distinct brain ROIs). The final step in our analysis involves the computation of measures that quantify the density and structure of the CRP. Notice that the CRQA process follows the same steps shown in the block diagram presented in Fig. 1, except for a difference in the estimation of the joint m and τ , as described above. The CRQA measures employed by our method rely on the diagonal and vertical lines in a CRP. Relevant notation in the representation of RP's, concerns the frequency distributions of the diagonal line lengths P ε k (l) = {l i = 1, . . . , N l }, where N l is the number of diagonal lines, for each diagonal parallel to the main diagonal line of identity (LOI). Each diagonal line consists of the points C i,j (ε) with j − i = k. If k = 0 the line is the LOI, if k > 0 the diagonal is located above the LOI, whilst if k < 0 it is located below the LOI. • Recurrence Rate (RR): Recurrence rate is the simplest measure that can be extracted from a CRP, and quantifies the density of cross recurrence points as follows, where M denotes the dimension of C. Notice that M is smaller than the number of samples in the time series K. • Determinism (DET): Determinism relies on a key feature of the elements comprising a recurrence plot, namely, the presence of structures. In the case of a deterministic system, these structures appear as long diagonal lines, which implies highly correlated processes. Alternatively, in weakly similar processes, as in a system displaying chaotic behavior, these binary matrices consist of very short diagonal lines or even isolated recurrence points.
DET is defined as the ratio of the recurrence points forming diagonal lines of at least length l min over all the recurrence points, The threshold l min excludes the diagonal lines which are formed by the tangential motion of the phase space trajectory. • Average diagonal line length (L): L represents the average time that two segments of the examined trajectory are close to each other, and is related to the degree of divergence of trajectory segments of length l, • Maximal diagonal line length (L max ): An extension to the interpretation of the diagonal lines via their divergence, is for how long the segments of a trajectory converge. The longer the convergence time of the trajectory segments, the longer the diagonal lines, where N l = l≥lmin P ε k (l) is the total number of diagonal lines. • Entropy (ENTR): Entropy is computed according to Shannon [25] and is based on the probability p(l) = P ε k (l) N l of finding a diagonal line of length exactly equal to l. In this context, an abundance of short diagonal lines would indicate low complexity in a pair of regions that display uncorrelated activity. Entropy is defined by, • Maximal length of vertical lines (V max ): V max considers the vertical lines in a CRP, which reflect the tendency of the corresponding trajectory to follow tangential motions in phase space. High V max values indicate that two time series are trapped in a dependency for a long time, that is, their states are highly similar. V max is calculated as follows, where N ν denotes the absolute number of vertical lines. RR and DET take values in the range of [0, 1], L, L max and ENTR can take values greater than 1, with higher values indicating greater similarity in dynamic patterns of activity between the two regions examined each time. Finally, small V max values imply rather isolated recurrence points in the vertical direction, while high V max values indicate large trapping times between the two systems.

E. Comparison Metrics
The relative sensitivity of CRQA features (RR, DET, L, L max , V max , ENTR) in differentiating NPSLE patients from healthy volunteers, as correlates of visuomotor performance, was assessed through comparisons with three sets of indices: (a) Nodal RQA measures computed for each of the 16 ROIs, (b) a conventional static FC metric (zero-order Pearson correlation between two time series computed for all 120 connections among the 16 ROIs), (c) dynamic FC based on functional identification of two temporal states (for all 120 unique connections among the 16 ROIs in each state).
Nodal RQA. From the set of 13 possible nodal RQA-based features the following three did not contain any undefined values and were used in further analyses: recurrence rate (RR), the maximal diagonal line length (L max ) and the maximal vertical line length (V max ), as these are defined in [23]. The higher rate of undefined values in computing nodal RQA metrics may have been caused by the use of a common radius ε = 0.35 of the maximum Euclidean distance between two state vectors, in order to render the results somewhat comparable to CRQA metrics, although this value may not have been appropriate for all 16 ROIs.
Temporal-based Dynamic FC. For each individual window, an assignment of nodes to functional modules was achieved using the Louvain greedy algorithm [24]. These module assignments were then used for the calculation of the nodal participation coefficient and within-module degree zscore. The joint histogram of these metrics [7] served as the set of features in a k-means clustering model (k = 2) in order to identify integrated and segregated states. The identification using clustering uses the features derived from all of each subject's time-resolved windows. The cluster with the higher mean participation coefficient is considered to represent the predominately integrated state. State graphs are calculated as the edge-wise median of each state's assigned networks. Median graphs consist of the 120 ROI-ROI connections between the above mentioned 16 ROIs. Pairwise connectivity indices (Pearson correlation coefficients) between these ROIs as they participate in each of the two representative networks, are compared with corresponding indices derived from the proposed CRQA method, as well as static FC.

III. STATISTICAL ANALYSIS
Differences between the two study groups (NPSLE vs. HC) were assessed using three sets of one-way MANOVAs performed separately for each pair of ROIs. The first set of MANOVAs were performed on ROI-to-ROI CRQA measures and included all six CRQA measures (RR, DET, L, L max , V max , ENTR) as complementary indices of connectivity (for a total of 120 tests). The second set of MANOVAs was performed on nodal RQA measures (RR, Lmax, Vmax, for a total of 16 tests). The third set of MANOVAs was performed on dFC measures (pairwise connectivity during the integrated state and corresponding index during the segregated state, for a total of 120 tests). A fourth set of comparisons entailed univariate ANOVAs on the static FC metric (total of 120 tests). Given that Bonferroni correction for multiple comparisons would result in prohibitive p-values, we relied primarily on multivariate effect sizes (η 2 ) as a measure of the magnitude of group differences.
Associations between CRQA measures (and comparison static and dynamic FC metrics) and patient performance on the Grooved Pegboard task were assessed using Pearson's correlation coefficient evaluated at p < 0.0001 uncorrected (corresponding to a critical r > 0.55).

A. Visuomotor capacity of NPSLE patients
Evaluated against age-and education-adjusted population norms, the NPSLE group displayed impaired average performance on the Grooved Pegboard test using the dominant hand (mean z = −1.51, SD = 2.03). Importantly, individual scores displayed adequate dispersion with 50% of patients performing in the significantly impaired range (z < −1.5).

B. Differences between NPSLE and HC groups on FC metrics
CRQA-based FC. MANOVA analysis revealed at least medium effect sizes (η 2 > 0.61) in 40 out of the 120 ROI pairs, with values ranging from 0.061 to 0.214. For 31 region pairs inspection of average CRQA measures indicated greater connectivity in the NPSLE as compared to the HC group (ref. Fig. 2A). The opposite trend was found in the remaining 9 pairs involving a coupling of the right paracentral lobule with 5 other parietal regions (bilateral MI, SMA, and precuneus, and the right SI) and the superior frontal gyrus bilaterally (which includes the dorsal premotor area; ref. Fig. 2C). The most frequently featured regions within the 40 pairs were four parietal lobe ROIs bilaterally: SI, SPL, angular gyrus, and PCL. Inspection of the results of univariate group comparisons revealed that group differences approached significance only for determinism within four region pairs, all including the angular gyrus, bilaterally (0.001 < p < 0.05, uncorrected; associated with large effect size ranging from 0.134 to 0.155).
Nodal RQA. Large effect sizes (ranging between 0.09 and 0.146) were found in 9/16 ROIs in the MANOVAs comparing the two clinical groups on the three RQA metrics (RR, Lmax, Vmax). In all cases NPSLE patients displayed larger values than the HC group.
Static FC. The effect size associated with one-way ANOVAs on ROI-to-ROI connectivity assessed using conventional correlation coefficients was greater than 0.61 in only 4 out of the 120 ROI pairs, all of which also displayed large effects sizes (NPSLE>HC) in the MANOVAs on CRQA measures (ref. Fig. 2B, D). However, adding sFC metrics to the 6 CRQA measures in the same MANOVA model, improved the effect size (from medium to large) for 6 out of the 120 ROI pairs.
Temporal-based Dynamic FC. Large effect sizes (ranging between 0.063 and 0.11, in all cases indicating NPSLE>HC), were found for 11/120 ROI pairs in the MANOVAs assessing group differences on connectivity metrics computed during the segregated and integrated states. Interestingly, large effect sizes in univariate ANOVAs were found for all 11 connections during the segregated state (involving connections between primary somatomotor regions and the SMA, the SMA and PCL, bilaterally, and also between the left and right angular gyri, and the left angular gyrus and precuneus) and in only 3 of the 11 connections during the integrated state (between the SMA and PCL).
Notice that, these figures were created through the BrainNet tool [26].

C. Associations between FC metrics and visuomotor capacity in NPSLE
CRQA metrics. Correlations with Grooved Pegboard performance were computed with the determinism index of FC in the 40 region pairs, where substantial group differences were found. Significant negative associations were detected in three of the ROIs pairs, where NPSLE patients displayed reduced DET values as compared to controls: PCL-SI (r = −0.626), PCL-MI (r = −0.566), and PCL-SMA (r = −0.565; all in the right hemisphere). Marginally significant associations (−0.55 < r < −0.50) were found between Grooved Pegboard performance and DET in the following ROI pairs: right PCLright SFG, right PCL-left SFG, right PCL-left precuneus, right PCL-left SMA (see lower 2 rows of Fig. 3). Additionally, significant associations with performance were found in two ROI pairs where DET values were larger among NPSLE patients as compared to the HC group: right SI-left SI (r = −0.58) and right SI-right MI (r = −0.59) (see upper row of Fig. 3). These findings suggest that higher DET was associated with faster performance on the Grooved Pegboard test with the dominant hand, reflecting higher visuomotor capacity.
Temporal-based Dynamic FC. Significant associations of performance with dFC estimates (p < .001) were in the oppo-site direction (positive) for the following ROI pairs: left MIleft Precuneus and right SMA-right Precuneus (r = 0.574 and r = 0.618, respectively, with connectivity in the segregated state), left MI-right Precuneus (r = 0.570, integrated state), left PCL-left Precuneus (r = 0.637 and r = 0.575, during the segregated and integrated states, respectively).
Nodal RQA and static FC metrics. Associations of performance with nodal RQA and sFC estimates did not approach significance anywhere (|r| < 0.32).

V. DISCUSSION
The present results highlight the potential utility of the CRQA technique as a complementary means to assess functional brain connectivity from rs-fMRI data. CRQA not only proved to be more sensitive in detecting signs of aberrant FC than static FC, but also helped to account for concurrent performance deficits in a clinically challenging patient group. Moreover, CRQA was found to be more sensitive than an established dFC approach, relying on temporal fluctuations of dynamic states, in identifying functional connections that may be reduced in NPSLE. Our results underscore the potential benefits from taking into account recurrent patterns in the signals recorded from distinct brain regions, to derive indices of cross-regional similarities in hemodynamic activity as it evolves over time. These indices, and especially determinism, may reveal functionally significant properties of self-organizing neuronal signaling, at rest, between regions normally involved in visuomotor coordination and control. Determinism (DET) characterizes the proportion of elements in a given recurrence plot that form long diagonal lines. The longer the diagonal lines the higher the determinism values. A CRP with long diagonal structures denotes trajectories with long autocorrelated times and further, time series with highly similar state vectors. Thus, even though rs-fMRI time series are not stationary, specific sets of densely interconnected brain regions appear to display similar features over time.
Our findings bear on the topic of changes in functional organization that can take place in NPSLE, even in the absence of remarkable structural changes on conventional MRI. Importantly, these changes appear to be crucial for relatively preserved motor function in this population. As a group the present cohort of NPSLE patients performed significantly below age-and education-adjusted population norms on a test of visuomotor control and coordination (Grooved Pegboard test). They also displayed reduced average functional connectivity (indexed by DET) between the paracentral lobule and several other fronto-parietal regions (SMA, MI, SI, dorsal premotor area). Importantly, those patients with more preserved connectivity performed better than their peers who displayed relatively lower connectivity. Visuomotor performance was additionally predicted by increased functional connectivity between primary sensorimotor areas (SI and MI-both relatively and in comparison with healthy controls). Thus, it appears that compensatory increases in co-recurrent patterns of BOLD fluctuations in fronto-parietal areas known to be involved in the sensorimotor and visuomotor integration and control contribute to the relative preservation of visuomotor capacity.
In light of the functional-systems approach adopted in the present study, it may not be wise to attempt direct comparisons regarding clinical group differences with the results of the only other rs-fMRI study that specifically assessed NPSLE patients [27]. In that study, network-level analyses revealed evidence of widespread hyper-as well as hypo-connectivity involving primarily the Default Mode, Sensorimotor, and Central Executive Networks. However, some notable similarities with this earlier work can be observed, considering that the regions examined in the present study are components of the sensorimotor (MI, SI, SMA, dorsal premotor, SPL, PCL) and Default Mode Networks (angular gyrus), and Executive Networks (precuneus). Specifically, the authors of [27] also found increased connectivity of the sensorimotor cortex among NPSLE patients as compared to HC. However, they reported reduced connectivity of the angular gyrus with the precuneus and other Default Mode Network (DMN) regions, whereas we observed increased connectivity of the angular gyrus primarily with regions outside this network (SI, MI, SMA, SPL).

A. Comparison with prior work
The performance of the CRQA method as an alternative to the computation of pairwise regional interrelations was compared with the most well-established way of approximating brain FC, namely Pearson coefficients between the entire time series. Compared to this method, CRQA proved to be more sensitive in detecting interrelations of potential clinical importance, given that they vary systematically among NPSLE patients as compared to the HC group.
The idea of applying the RQA theory on fMRI time series is not new, albeit limited in previous studies to the RQA technique with task-related fMRI data. In a seminal work [28], Bianciardi et al. compared the performance of the RQA method with the conventional general linear model (GLM), applied on simulated and real fMRI data, to extract regional patterns of hemodynamic activity. As RQA does not make any distribution assumptions, does not rely on signal stationarity, and is largely unaffected by variations in the shape and timing of the underlying neuronal and hemodynamic responses, it proved to be more robust in detecting recurrent patterns in the recorded time series. Moreover, an evaluation of the univariate RQA is presented in [29]. In this study, researchers aimed to examine how effective is the RQA in detecting chaotic behaviours from EEG recordings. However, on the one hand this study refers to a different modality than the rs-fMRI, and on the other hand, it examines the effectiveness of RQA from different perspectives. The current work takes RQA one step further to assess cross-regional recurrent patterns of hemodynamic fluctuations as they evolve over time in the absence of external stimulation. Importantly, we showed that the degree of cross-regional similarity in such recurrent patterns may vary systematically in disease and parallels the impact of the disease on behavioral performance.
In our results, signals recorded in ROIs where notable group differences were also detected on CRQA indices, also varied systematically on nodal-RQA measures (NPSLE>HC in MI, SI, dorsal premotor, angular gyrus, SPL and PCL). Interestingly, these regions featured in pairs displaying both increased (e.g., MI, SI, dorsal premotor, angular gyrus, SPL) and reduced CRQA-based connectivity in NPSLE vs. HC (PCL). Thus it appears that cross-recurrence and auto-recurrence hemodynamic phenomena may be dissociable. This notion is supported by the finding that only CRQA measures correlated significantly with visuomotor performance.
There were further notable dissociations between the results of CRQA and the comparison dynamic FC approach used in the present study. The latter method relies on identifying time-varying spatial patterns of connectivity, as opposed to the CRQA technique which focuses on recurrent signal similarities between two regional time series at a time. Interestingly, CRQA appeared to be more sensitive than temporal-based dFC to instances of reduced connectivity among NPSLE patients, as compared to the healthy control group. Both methods were capable of identifying instances of hyperconnectivity in the NPSLE group. Notably, relative increases in connectivity strength in these ROI pairs correlated strongly, albeit positively, with time to complete the Grooved Pegboard task (indicating slower performance). Therefore, it appears that the observed hyperconnectivity in the NPSLE group is a manifestation of compensatory neural reorganization processes involving the ipsilateral MI, SMA and PCL (all components of the somatomotor resting-state network) and the precuneus (which is part of the executive network).

B. Limitations and future work
The present work has some limitations that need to be investigated in a thorough future study. Firstly, FC was examined for only a subset of brain regions which were defined a priori based on functional criteria (i.e., their known role in supporting visuomotor coordination and control in the healthy brain) and prior evidence of disturbed cerebral perfusion in NPSLE patients [12]. Consequently, we did not assess the capacity of CRQA measures of FC to reveal the overall organization of resting-state hemodynamic activity in both patients and healthy volunteers at the level of known rs-fMRI networks and, consequently, our results are not directly comparable to earlier findings [27].
Secondly, each CRQA measure was treated as a separate, complementary index of FC, and we did not attempt to combine these measures into a single vector representing a composite measure of FC, which could also include the conventional static index (correlation coefficient). As a future work, integration of complementary CRQA measures into a composite score for each pair of brain regions would render them more amenable to network-level analyses, e.g., using graph theory.
Furthermore, our results should be compared to and integrated with indices of global connectivity (such as the Intrinsic Connectivity Contrast; ICC) to appreciate functional changes taking place in NPSLE more comprehensively. Moreover, integration of connectivity with hemodynamic measures is important to help to understand the pathophysiological changes that trigger the observed changes in FC.
Finally, associations between CRQA measures and performance were computed only among NPSLE patients. It may be the case that such associations are also evident in the healthy brain or, alternatively, that they emerge (or become accentuated) as a compensatory phenomenon when the brain networks that normally support a particular cognitive function are affected by brain pathology.

VI. CONCLUSIONS
This study focuses on the cross recurrence quantification analysis as a method for estimating similarities in hemodynamic activity between brain regions at rest. CRQA was applied to functional associations between brain regions known to be involved in visuomotor coordination and control. These associations are supported by major fiber bundles traveling through the semioval center, which appears to be underperfused in NPSLE [12]. In this context, CRQA measures, and particularly determinism, were successful in revealing a pattern of aberrant FC that was not evident in analyses relying on static FC metric (Pearson's correlation). Crucially, the observed changes in FC accounted for systematic variability in the degree of visuomotor impairment displayed by individual NPSLE patients. Overall, feature-based connectivity estimation methods such as CRQA indicated a strong potential in detecting differences between the two populations especially when compared to static FC. Comparison dynamic FC also exhibited strong differentiating power as well as correlations with visuomotor performance which static FC did not. The apparent strong-suit of CRQA-based FC compared to previously used dynamic FC, according to the results presented, is its ability of highlighting connections weaker in NPSLE patients. On the basis of this preliminary work, we surmise that CRQA measures can be integrated with other dynamic indices of FC and expanded to assess whole-brain, networklevel interdependencies in hemodynamic activity at rest.