Parcellation of the human amygdala using recurrence quantification analysis

Several previous attempts have been made to divide the human amygdala into smaller subregions based on the unique functional properties of the subregions. Although these attempts have provided valuable insight into the functional heterogeneity in this structure, the possibility that spatial patterns of functional characteristics can quickly change over time has been neglected in previous studies. In the present study, we explicitly account for the dynamic nature of amygdala activity. Our goal was not only to develop another parcellation method but also to augment existing methods with novel information about amygdala subdivisions. We performed state-specific amygdala parcellation using resting-state fMRI (rsfMRI) data and recurrence quantification analysis (RQA). RsfMRI data from 102 subjects were acquired with a 3T Trio Siemens scanner. We analyzed values of several RQA measures across all voxels in the amygdala and found two amygdala subdivisions, the ventrolateral (VL) and dorsomedial (DM) subdivisions, that differ with respect to one of the RQA measures, Shannon’s entropy of diagonal lines. Compared to the DM subdivision, the VL subdivision can be characterized by a higher value of entropy. The results suggest that VL activity is determined and influenced by more brain structures than is DM activity. To assess the biological validity of the obtained subdivisions, we compared them with histological atlases and currently available parcellations based on structural connectivity patterns (Anatomy Probability Maps) and cytoarchitectonic features (SPM Anatomy toolbox). Moreover, we examined their cortical and subcortical functional connectivity. The obtained results are similar to those previously reported on parcellation performed on the basis of structural connectivity patterns. Functional connectivity analysis revealed that the VL subdivision has strong connections to several cortical areas, whereas the DM subdivision is mainly connected to subcortical regions. This finding suggests that the VL subdivision corresponds to the basolateral subdivision of the amygdala (BLA), while the DM subdivision has some characteristics typical of the centromedial amygdala (CMA). The similarity in functional connectivity patterns between the VL subdivision and BLA, as well as between the DM subdivision and CMA, confirm the utility of our parcellation method. Overall, the study shows that parcellation based on BOLD signal dynamics is a powerful tool for identifying distinct functional systems within the amygdala. This tool might be useful for future research on functional brain organization. Highlights A new method for parcellation of the human amygdala was developed The ventrolateral and dorsomedial subdivisions of the amygdala were revealed The two subdivisions correspond to the anatomically defined regions of the amygdala The two subdivisions differ with respect to values of entropy A new parcellation method provides novel information about amygdala subdivisions

ABSTRACT Several previous attempts have been made to divide the human amygdala into smaller subregions based on the unique functional properties of the subregions. Although these attempts have provided valuable insight into the functional heterogeneity in this structure, the possibility that spatial patterns of functional characteristics can quickly change over time has been neglected in previous studies. In the present study, we explicitly account for the dynamic nature of amygdala activity. Our goal was not only to develop another parcellation method but also to augment existing methods with novel information about amygdala subdivisions. We performed state-specific amygdala parcellation using resting-state fMRI (rsfMRI) data and recurrence quantification analysis (RQA). RsfMRI data from 102 subjects were acquired with a 3T Trio Siemens scanner. We analyzed values of several RQA measures across all voxels in the amygdala and found two amygdala subdivisions, the ventrolateral (VL) and dorsomedial (DM) subdivisions, that differ with respect to one of the RQA measures, Shannon's entropy of diagonal lines. Compared to the DM subdivision, the VL subdivision can be characterized by a higher value of entropy. The results suggest that VL activity is determined and influenced by more brain structures than is DM activity. To assess the biological validity of the obtained subdivisions, we compared them with histological atlases and currently available parcellations based on structural connectivity patterns (Anatomy Probability Maps) and cytoarchitectonic features (SPM Anatomy toolbox). Moreover, we examined their cortical and subcortical functional connectivity. The obtained results are similar to those previously reported on parcellation performed on the basis of structural connectivity patterns. Functional connectivity analysis revealed that the VL subdivision has strong connections to several cortical areas, whereas the DM subdivision is mainly connected to subcortical regions. This finding suggests that the VL subdivision corresponds to the basolateral subdivision of the amygdala (BLA), while the DM subdivision has some characteristics typical of the centromedial amygdala (CMA). The similarity in functional connectivity patterns between the VL subdivision and BLA, as well as between the DM subdivision and CMA, confirm the utility of our parcellation method. Overall, the study shows that parcellation based on BOLD signal dynamics is a powerful tool for identifying distinct functional systems within the amygdala. This tool might be useful for future research on functional brain organization.

INTRODUCTION
The amygdala plays a crucial role in several emotion-related functions, such as reinforcement learning (Murray, 2007;Wassum and Izquierdo, 2015), social behavior (Adolphs, 2010), and emotion regulation (Blair et al., 2007), in both animals and humans (LeDoux, 2007;Janak & Tye, 2015). Moreover, dysfunction of the amygdala is associated with a variety of psychiatric disorders, such as autism spectrum disorders (ASDs), anxiety disorders, and depression (Phillips et al., 2003;Phelps & LeDoux, 2005). It has been shown that the anatomic basis of the diverse functions of the amygdala is its complex internal structure.
According to animal studies, the amygdala comprises a set of small nuclei located very close to each other (Swanson & Petrovich, 1998). It is very likely that distinct groups of nuclei, comprising subdivisions of the amygdala, process different kinds of information. Thus, identifying the role of different subdivisions of the amygdala in humans might help in understanding both healthy emotional functioning and mechanisms of psychiatric disorders.
However, distinct subdivisions of the human amygdala have not been clearly identified. The roles of the amygdala subdivisions are difficult to examine in humans mainly because of the lack of a generally accepted method of parcellating this structure (for a review see Kolada et al. 2017, pp. 123-124). In effect, several research groups have attempted to develop a common method that divides the human amygdala into smaller subparts on the basis of their unique structural or functional properties.
Structural parcellations are typically based on the idea that a given subdivision should exhibit distinct cytoarchitectonics, topography and structural connectivity. Cytoarchitecturebased parcellations involve differentiating subparts of the brain structure on the basis of the quantity, density and localization of different types of neuronal cells. To the best of our knowledge, only one study has used a cytoarchitectonic approach to delineate small subdivisions of the amygdala ). Amunts and colleagues performed postmortem analysis of brain tissue and defined borders of three subareas of the amygdala: the laterobasal, centromedial and superficial areas. The results of this analysis are available in the SPM Anatomy toolbox for future MRI-based research (Eickhoff, et al., 2005). Topographybased parcellations rely on the manual delineation of subdivisions of brain structures on the basis of the locations of other brain areas, blood vessels, and cerebrospinal fluid spaces located nearby. At least two previous studies have used this type of parcellation scheme for the amygdala (Entis et al.,2012;Tyszka & Pauli, 2016). However, topography-based methods are researcher-dependent and may lead to inconsistent results. Finally, parcellations based on structural connectivity are grounded in the assumption that different subdivisions of the brain areas correspond to different parts of brain circuits; thus, they might be characterized by different patterns of connectivity. In human research, it is possible to conduct structural connectivity analysis using diffusion tensor imaging (DTI). Regarding the amygdala, there are several ways to apply structural connectivity analysis (Bach et al.,2011;Solano-Castiella et al., 2010;Saygin et al., 2011). As connectivity-based methods largely rely on the way of estimating the connections between structures that has been selected, they yield distinct numbers and sizes of obtained parcels.
Functional parcellations are focused on the functional specialization of brain areas or the functional connectivity between them. If a given brain area is responsible for a specific function, it is possible to prepare a fMRI localization task to selectively activate this area. Although there is no specific task for the selective activation of distinct amygdala subdivisions, it is possible to perform meta-analysis of data derived from a large number of experiments, in which even low levels of activation of the amygdala have been detected. There are at least two works that have used this method for amygdala parcellation (Bzdok et al., 2013;Yang et al., 2016).
Unfortunately, parcellations on the basis of functional specialization (especially meta-analysis) require a very large number of experiments. An alternative approach to the functional parcellation of the amygdala is parcellation based on its functional connectivity with other brain regions. Using resting-state fMRI (rsfMRI), which relies on the analysis of spontaneous fluctuations in BOLD signals in the absence of an explicit task, it is possible to identify functional connectivity between brain areas (i.e., functional interactions between them; Fox and Raichle 2007). To date, functional connectivity analysis has demonstrated that rsfMRI might be a valuable tool for dividing the amygdala into two (Mishra et al., 2014) or three parts (Bickart et al., 2014;Caparelli et al., 2017).
Although the above methods for functional parcellation have provided valuable insight into the functional heterogeneity in the human amygdala, they are encumbered by a common disadvantage: they have neglected the possibility that the spatial patterns of functional characteristics may quickly change over time. In recent years, however, there has been a growing number of neuroimaging studies that take into account the dynamics of changes in brain activity (Izhikevich et al., 2007). In fMRI studies, dynamics can be assessed by the sliding window approach in connectivity analysis (Hutchison, et al., 2013) or modeling with differential equations (Deco et al.,2008). This approach seems to be fruitful in the field of neuroimaging, and recently, it was used for parcellation of the human thalamus (van Oort et al., 2018;Kumar et al., 2017).
In the present study, we explicitly account for the dynamic nature of amygdala activity. We derive state-specific amygdala parcellations using recurrence quantification analysis (RQA), a nonlinear method derived from dynamic systems theory. Patterns of recurrence in a signal (understood as revisiting a previously visited state; Marwan & Webber, 2015) reflect the properties of underlying dynamics of a system, which generate the signal. RQA can be performed with short time series and allows us to assess multidimensional relations in signal features reflecting different aspects of recurrence. This method was successfully used for detecting task-related activity and for estimating functional connectivity patterns on the basis of fMRI data (Bianciardi, et al., 2008;Lombardi et al., 2017).
As animal and human research indicate that the amygdala consists of some subdivisions that differ in functional properties, we can expect that the dynamics of functional MRI signals differ by subdivision. It is possible to verify this idea by assessing rsfMRI data. During the resting state, different brain subsystems reveal relatively stable patterns of synchronous activity (De Luca et al., 2006;Damoiseaux et al., 2006). Thus, it might be possible to compute features of the dynamics of the BOLD signal from the resting state in all voxels within the amygdala and delineate borders between functionally distinct subdivisions. In the present work, we used rsfMRI data and RQA to estimate functionally distinct subdivisions of the amygdala based on the temporal characteristics of the signals. Analyses were performed on three datasets. To assess the biological validity of the obtained subdivisions, we compared these subdivisions with histological atlases and currently available parcellations (i.e., the anatomy probability map function based on structural connectivity and the SPM Anatomy toolbox). Moreover, we examined the cortical and subcortical projections of the subdivisions using accessible restingstate data. Our main goal was not only to develop another parcellation method but also to augment existing methods with novel information about the subdivisions. The results of the analysis on how patterns of activations within subdivisions change over time may reveal the dynamical organization of the functional subsystems that generate the signals.

Subjects.
In total, 102 subjects participated volunteered in the study and gave their informed consent on the documents prepared according to guidelines from the local ethics authorities from the University of Social Sciences and Humanities (SWPS University) in Warsaw and the University of Warsaw following guidelines of the Code of Ethics of the World Medical Association (Declaration of Helsinki). Individuals were recruited via an internet-based survey posted on social media. To participate in the study, the candidates had be aged 21 -35 years and right-handed; those who had metal medical devices located in the body, claustrophobia, psychiatric or neurological impairments and tattoos or permanent make-up were excluded from the study.
Data: All datasets were acquired using a Siemens Trio 3T scanner equipped with a 32-channel head coil. The scan parameters were as follows: Finally, the data were spatially smoothed with an 8 mm full width at half maximum Gaussian kernel and denoised with a bandpass filter between 0.01 and 0.1 Hz to minimize the effects of physiological noise.
Data analysis -pipeline. Figure 1 presents the general schema of the analysis. First, the parameters for the RQA were estimated to create the recurrence plots (RPs), and the RQA measures were computed on the signal from each voxel from all subjects included in the first dataset. Voxels within the amygdala were then grouped into clusters based on the RQA measures values (orange frames on Figure 1), and the solution with the best internal validity coefficients was selected. To confirm the results, the analysis was run on another dataset (second dataset) using the same parameters and previously selected RQA measures (blue frames). Finally, the solutions from both datasets were compared with those of existing parcellations, and we investigated the functional connectivity patterns of the obtained subdivisions (green frames) with the third dataset to validate the anatomical reliability of the obtained solution and obtain additional insight into the functional organization of the human amygdalae. In the next subsections of this paper, there are detailed descriptions of the stages. Data analysis -masks. Voxels representing the amygdala were extracted from the structural and functional data by the anatomy probability map function based on the structural connectivity (APMC) mask prepared by Bach and collaborators (2011). The outline of this mask was drawn manually for each individual and then parcellated into two subdivisions on the basis of the structural connectivity patterns derived by diffusion weighted imaging (DWI). This mask was chosen because it allows the results obtained on the basis of functional features to be compared directly with those obtained by parcellation performed with structural data. Since the structural representation of the brain areas is not completely congruent with the functional representation, a comparison of these two perspectives might enrich our analysis.  Data analysis -Clustering. We designed six matrices, one for each measure (example shown in Figure 3). The rows represent the subjects, while the voxels in the MNI space correspond to consecutive columns. We assumed that the distribution of the RQA measures in the population could be a good indicator of these dynamics. As the dynamics of the signal in different subdivisions should be similar in the whole population, it might be possible to classify voxels into clusters representing functionally separate subdivisions on the basis of at least one RQA measure. To determine which measure is the best for this purpose, we tested all of them and selected the most appropriate one. For classification purposes, we used hierarchical clustering, which divides the dataset into smaller subparts on the basis of similarity. The similarity between vectors was assessed with ward metrics. Additionally, we imposed spatial constraints: two points could connect only if they were immediate neighbors to each other (voxels located next to each other), and the number of clusters was set to two. Data Analysis -Connectivity Analysis. To determine whether the resulting subdivisions of the amygdala corresponded to anatomical groups of nuclei, which can be characterized by distinct connectivity patterns, we performed ROI-to-ROI functional connectivity analysis on the third dataset using the common part of the obtained subdivisions in the first and second datasets as masks of the different amygdala subparts. Our analysis was conducted with the CONN toolbox, and we used the Harvard-Oxford Atlas for the cortical and subcortical regions.  Table 1   The distribution of the entropy values in the amygdala voxels was left skewed ( Figure   5). This pattern suggests that the values of entropy were higher in the VL part of the amygdala than in DM part of the amygdala (Figure 6).    Figure 7. All computations were performed using the SciPy toolbox from Python 3.7.

Repeated Analysis on the Second Dataset
To verify the stability of the solution, the whole procedure was carried out on the second dataset using the same RQA parameters, and Shannon's entropy of diagonal lines was used to differentiate the amygdala subdivisions. The results of this confirmatory analysis are presented in Figure 8. Similar to the previously obtained results, the VL part of the amygdala (left: 128 voxels, right: 101 voxels) had a larger volume than did the DM part (left: 60 voxels, right: 51 voxels). The distribution of the entropy values in the amygdala voxels in the second dataset ( Figure 9) had a similar tendency to be left-skewed, but in general, the values were smaller than those in the first dataset. For that reason, some voxels with negative values of entropy might have been classified by the algorithm as being in the VL subdivision ( Figure 10). Thus, the spatial distributions of the entropy values for the two datasets seemed to be reasonably similar.
The mean square error of the entropy values was smaller than 5%.   Table 3. Mean values of entropy in the amygdala subdivisions obtained for the second dataset. Fig. 11. Mean values of entropy with standard errors for the amygdala subdivisions obtained for the second dataset.

Comparisons of the Obtained Solutions
Spatial Representation. Here, we present a quantitative comparison of the parcellation solution for the first dataset, the parcellation solution for the second dataset and those of existing parcellations: structural connectivity-based parcellation (AMPC mask) by Bach et al. (2011) and the cytoarchitectonic atlas of the amygdala by Amunts et al. (2005) generated with the Anatomy toolbox (AT mask).
Additionally, we show the similarities between first and second dataset results (voxels that were classified to the same subdivisions in both analyses, see Figures 12 and 13).  Table 4). Comparison with existing parcellations -external validation. In addition to the visual comparison, we compared the two solutions using external validation measures, taking the solution obtained from the second dataset as the ground truth (Table 5). Regardless of the measured selected, the similarity between the datasets was approximately 70%.  Table 5. Comparison of the parcellations for the first and second datasets in terms of the external validation measures (the parcellation for the second dataset served as the ground truth for the computation of the external validation measures for the first dataset).
Moreover, we compared our parcellations with existing parcellations. As the masks used in the analysis differed in terms of the spatial distribution, the voxels that were localized in the amygdala in one mask could be assigned to another brain area (for example, in the comparison between the APMC mask and AT mask). For that reason, when we compared our results to those of existing parcellations, we needed to take into account this variability. We defined a limited comparison as one that compares only voxels that are included as parts of the amygdala in both masks, and a full comparison as one that takes into account all the voxels of the masks compared (in Figures 14 and 15, the regions are shaded in green and blue). We also provided a full comparison, taking into account all areas present in both masks (the whole areas in Figures   14 and 15). Coordinates: x = -21, y = -4, and z = -20; views: sagittal, coronal, and axial. Fig. 15. Comparison of the parcellations obtained for the first and second datasets with the parcellation obtained on the basis of cytoarchitecture by Katrin Amunts and colleagues (2005). The blue and green regions were classified in the same way in both parcellations; the rest of the colors denote that the solution was classified differently in the two datasets or that the voxel cannot be considered a part of the amygdala in one of the masks. Coordinates: x = -21, y = -4, and z = -20; views: sagittal, coronal, and axial.
The results of both the limited and full comparisons between the parcellation of the first dataset and the parcellation derived on the basis of structural connectivity patterns were the same and revealed a large degree of similarity between the two parcellations. The Fowlkes-Mallows score reached the level of approximately 75% with respect to the APMC mask, while the value of this indicator was much smaller than that obtained with respect to the SPM Anatomy toolbox mask. There are two reasons for this finding. First, the masks compared differ in terms of the spatial location of the voxels assigned to the amygdala -the mask from the Anatomy toolbox is significantly larger than the mask proposed by Bach et al. (2011). Second, parcellation on the basis of cytoarchitectonic features was performed on the three subparts, while our results indicated two functional subsystems within the amygdala. The values of the indices are presented in Table 6. We also calculated validation measures for the parcellation solution for the second dataset. As the proportions of sizes of the subdivisions were similar between the results of the second dataset and the SPM Anatomy toolbox, the external validation measures showed better results for this parcellations than for the comparison with the APMC parcellation (for example, value of Fowlkes-Mallows score for the limited comparison with the SPM Anatomy toolbox, Table 7). On the other hand, this comparison in the full scenario is similar to the results obtained on the first dataset. However, the Fowlkes-Mallows score was higher for the parcellation with first dataset than for the APMC parcellation and the parcellation with the second dataset. This finding proves that the results from the first dataset are more similar to the parcellation based on structural connectivity, while the results from the second dataset are more spatially congruent with the parcellation based on cytoarchitectonic features. Connectivity results of the third dataset. We showed that the obtained results are similar to those of previously existing parcellations. Until now, we labeled the obtained subdivisions as ventrolateral and dorsomedial subdivisions. To be able to align these names with anatomical nomenclature, we decided to perform a functional connectivity analysis on yet another dataset (the third dataset) with masks obtained in our analysis to identify the patterns of connections with different brain areas. We use the common areas between the ventrolateral and dorsomedial parts (as shown in Figure 12) and treated them as seeds in our analysis. To identify the distinct connectivity patterns characteristic to a specific amygdala subdivision, we compared all connections for the VL and DM subdivisions and determined which ones were had the strongest connections with a particular subpart.
Then, we repeated this analysis, but the main seed was in the mask of the dorsomedial amygdala and presented brain areas with significantly higher correlations with the spontaneous activity of this part than with the resting-state activity of the ventrolateral part. Thus, the dorsomedial part of the amygdala had stronger correlations with the subcortical structures such as the putamen, nucleus accumbens, thalamus and some cortical regions, such as the frontal medial cortex and postcentral gyrus. All of the correlations are presented in Table 9. In figure 16, the brain areas with stronger connections with the ventrolateral part are presented in red, and the brain areas with stronger connections with the dorsomedial subdivision are shown in blue.   Table 9. Functional connectivity of the dorsomedial part based on the common part between the two datasets, obtained for the third dataset.

DISCUSSION
Our study introduced a novel approach to parcellate the human amygdala on the basis of the dynamic characteristics of the BOLD signals in the amygdala in voxels, as assessed by recurrence quantification analysis (RQA). We clearly partitioned the amygdala into two subdivisions: the ventrolateral (VL) and dorsomedial (DM) subdivisions. These names were given because of their location in the brain.
In our analysis, we assessed the system's trajectory in the phase space and, subsequently, created a recurrence plot (Eckmann et al., 1987) for each voxel. Then, on the basis of the RPs, we computed several known measures for recurrence quantification analysis (RQA). From the set of these measures, Shannon's entropy of the distribution of diagonal lines appears to be the most suitable for the differentiation of amygdala subdivisions on the basis of neuronal dynamics. It measures how unpredictable the trajectory of the analyzed system is. The larger the value of the entropy is, the more uniform the distribution of the lengths of the diagonal lines, and the more variable the lengths, which means that it is more difficult to predict how long the system will repeat a previously traveled sequence of states. Thus, Shannon's entropy of diagonal lines focuses on the stability and regularity of a signal. For that reason, in our analysis, entropy was assessed for different features of data, whereas in previous investigations in the field of fMRI neuroimaging, the measure of entropy was directly applied to the BOLD time series (Wang et al., 2014;Mikoláš et al., 2012).
Regarding the masks we obtained, compared with the DM part, the VL part can be characterized by higher values of entropy. This finding may suggest that the VL subdivision is influenced by more brain systems. We tested this idea with the functional connectivity analysis, where masks of the VL and DM subdivisions were treated as seeds. VL activity was found to be correlated extensively with the cortical, occipital, parietal, temporal and frontal regions, whereas DM activity predicted activity primarily in the subcortical structures, such as the striatum, thalamus, and hippocampus. These two major amygdala complexes, with differential connectivity patterns, were consistent with connectivity patterns of the basolateral amygdala (BLA) and centromedial amygdala (CMA) (Roy et al., 2008). In contrast, the CMA is essential for controlling the expression of emotional responses through projections to subcortical structures, including the thalamus, hypothalamus, striatum, and brainstem (LeDoux 2000; Davis & Whalen 2001). The fact that the VL seems to be similar to the BLA shows that the less stable dynamics (higher entropy of the length of repeated trajectories in the phase space) of the brain signal might be related to the contribution of evolutionarily younger brain subsystems associated with emotion expressions rather than the DM or the CMA.
The similarity in functional connectivity patterns between the VL and BLA, as well as between the DM and CMA, confirm the utility of our new parcellation approach. Moreover, our method yielded parcellations similar to previously described parcellations performed on the basis of structural connectivity patterns (Bach et al., 2011) and cytoarchitectonic characteristics . It is a well-known fact that the spatial representation of brain areas highly correlates with but also partly differs from the results of structural and functional analysis (Honey et al., 2008). For that reason, the comparison of our results with those of parcellation on the basis of structural connectivity revealed high similarity but not identical solutions.
Moreover, Bach and colleagues (2011) performed parcellation based on the structural connectivity patterns acquired via DWI and two regions with which connections differentiating the subdivisions of the amygdala were used. Nevertheless, the connectivity patterns of the masks obtained in our study correspond well with the brain circuits and networks of the BLA and CMA.
On the other hand, our results were less similar to those of parcellation based on cytoarchitectonic features than those of MRI-based research. This finding might be caused by the fact that the amygdala was partitioned into three parts in the work by Amunts et al. (2005), while we divided the amygdala into two subareas. Moreover, the mask generated with the SPM Anatomy toolbox was significantly larger than the masks we used for the purpose of our analysis. Thus, these two masks overlapped only partially, and these two parcellations were less similar to that in our results than was the APMC parcellation.

CONCLUSIONS
In summary, we present a new parcellation approach for the human amygdala. We used RQA, which revealed the dynamic characteristics of the functional subsystems within the amygdala. This method also allowed us to further explore the nature of the BOLD signals in other brain areas. Although few related attempts have been made in the field of fMRI research