## Abstract

A primary goal of many neuroimaging studies that use magnetic resonance imaging (MRI) is to deduce the structure-function relationships in the human brain using data from the three major neuro-MRI modalities: high-resolution anatomical, diffusion tensor imaging, and functional MRI. To date, the general procedure for analyzing these data is to combine the results derived independently from each of these modalities.

In this article, we develop a new theoretical and computational approach for combining these different MRI modalities into a powerful and versatile framework that combines our recently developed methods for morphological shape analysis and segmentation, simultaneous local diffusion estimation and global tractography, and nonlinear and nongaussian spatial-temporal activation pattern classification and ranking, as well as our fast and accurate approach for nonlinear registration between modalities. This joint analysis method is capable of extracting new levels of information that is not achievable from any of those single modalities alone. A theoretical probabilistic framework based on a reformulation of prior information and available interdependencies between modalities through a joint coupling matrix and an efficient computational implementation allows construction of quantitative functional, structural, and effective brain connectivity modes and parcellation.

This new method provides an overall increase of resolution, accuracy, level of detail, and information content and has the potential to be instrumental in the clinical adaptation of neuro-MRI modalities, which, when jointly analyzed, provide a more comprehensive view of a subject’s structure-function relations, while the current standard, wherein single-modality methods are analyzed separately, leaves a critical gap in an integrated view of a subject’s neuorphysiological state. As one example of this increased sensitivity, we demonstrate that the jointly estimated structural and functional dependencies of mode power follow the same power law decay with the same exponent.

## 1 Introduction

One of the great scientific questions that modern imaging systems are increasingly capable of addressing is the structure-function relationships in the human brain. Magnetic resonance imaging MRI, in particular, provides a unique methodology for noninvasively acquiring high-spatial-resolution anatomical data (HRA) that are capable of characterizing morphological variations in exquisite detail, diffusion tensor imaging (DTI) that enables the construction of neural pathways, and functional MRI (FMRI), and, in particular, resting-state functional MRI (rsFMRI), which can detect spatial-temporal changes in brain activity.

The general idea behind characterizing structure-function relationships is to relate quantitative measures of structural features of a system with the spatiotemporal patterns of parameters changes. For example, in neuroimaging using MRI (neuro-MRI), the goal is to characterize the relationship of brain morphology, neural pathways, and brain activity by combining morphological information from HRA data, structural connectivity from DTI data, and functional connectivity from rsFMRI.

The most common approach to the analysis of neuro-MRI data is to estimate morphology, connectivity, and spatiotemporal variations separately first, then combine them to characterize system “modes.” This is, for the most part, simply a practical matter: different methods are used for each analysis and so are applied separately. For example, in studies of structure-function relationships in the human brain. Deco et al. (Deco, Jirsa, McIntosh, Sporns, & Kötter, 2009; Deco, Jirsa, & McIntosh, 2011; Deco et al., 2014) use an empirical structural connectivity (SC) matrix derived from diffusion spectrum imaging (DSI, an extension of DTI; Wedeen et al., 2005) and tractography (Hagmann et al., 2008) and an empirical functional connectivity (FC) matrix based on Honey, Kotter, Breakspear, and Sporns (2007) and Honey et al. (2009). Honey et al. (2009) construct FC maps based on Pearson correlation between BOLD time series (a correlation analysis.) The correlation analysis is also used in Dosenbach et al. (2007). Cammoun et al. (2012) use a similar analysis with DSI at multiple scales, and Zalesky, Fornito, Cocchi, Gollo, and Breakspear (2014) use a time-resolved analysis with simple correlation analysis with sliding windows.

However, the great potential of these methods is compromised because of three significant limitations in current state-of-the-art analysis approaches. First, they are, for the most part, not quantitative and, as a consequence, often are not sufficiently accurate: (1) quantitation of morphology from HRA is typically limited to surfaces and is highly error prone; (2) existing DTI analysis methods are incapable of accurately characterizing voxels containing multiple fibers, leading to highly inaccurate estimation of both local anisotropy and global tractography; and (3) resting-state rsFMRI analysis methods are essentially ad hoc and have no well-characterized definition of a “brain mode,” much less a computational approach amenable to the ranking of modes and the assessment of error estimates. Second, the data from each method are analyzed independent of one another, ignoring the inherent constraints that couple anatomy, connectivity, and function. Third, these methods (DTI and FMRI in particular) are poorly posed inverse problems with a wide range of possible “solutions” to any single estimate, some of which are nonphysical. Current analysis methods rarely incorporate constraints that require solutions to be physically realizable.

In this article, we address all three of these deficiencies by developing a general theoretical and computational framework for the analysis of neuro-MRI data based on the reformulation of the problem in terms of our recently developed entropy field decomposition (EFD) method (Frank & Galinsky, 2016b), which combines information field theory (IFT) (Ensslin, Frommert, & Kitaura, 2009), with prior information encoded on spatial-temporal patterns using our recently introduced method of entropy spectrum pathways (ESP) (Frank & Galinsky, 2014). The EFD method is a probabilistic approach that provides a general framework for the integration of HRA, DTI, and FMRI analysis, including results provided by numerical simulations (Baxter & Frank, 2013; Balls & Frank, 2009) to constrain results to those that are physically realistic. We will employ our recently developed methods for each of the three modalities: (1) the spherical wave decomposition (SWD) method for shape analysis and segmentation (Galinsky & Frank, 2014); (2) the geometrical optics-based ESP method for simultaneous estimation of local diffusion and global tractography (GO-ESP) (Galinsky & Frank, 2015); and (3) the characterization and estimation of space-time rsFMRI brain activation patterns using EFD (Frank & Galinsky, 2016b). The spatial collocation of all of these different modalities is provided by using a novel nonlinear registration method based on a general symplectomorphic formulation that constructs flexible grid mapping between modalities regularized by ESP coupling (SYMREG-ESP) (Galinsky & Frank, 2016b). We then demonstrate how these three analysis methods can be fused together and be mutually informative within the overarching theory of EFD, allowing us to construct more quantitative estimates of brain structural and functional modes.

With the constantly increasing quality of data from neuroimaging methods using MRI (neuro-MRI), such as high-resolution anatomical imaging (HRA), diffusion tensor imaging (DTI), and resting-state functional MRI (rsFMRI), there is a significant and growing interest in the development of frameworks that can perform multimodal data fusion and analysis of these different modalities. Several recently developed methods are based on some modification or ad hoc augmentation of an independent component analysis (ICA) approach (Bell & Sejnowski, 1995; Hyvarinen & Oja, 2000), including Groves, Beckmann, Smith, and Woolrich (2011), Douaud et al. (2014) and Mueller et al. (2013), or the joint/parallel ICA approaches (see the review by Calhoun, Liu, & Adali, 2009). However the performance and applicability of ICA even for a single modality (rsFMRI) analysis is predicated on a number of assumptions, such as sparsity and independence, that are questionable and have been debated (Daubechies et al., 2009; Calhoun et al., 2013). More practically, and important, it has been shown that the ICA approach can fail in the simplest of idealized simulations (Calhoun et al., 2013; Frank & Galinsky, 2016b). This motivated our development of the EFD approach, which was demonstrably more accurate than ICA in its application to rsFMRI data (Frank & Galinsky, 2016b).

In our original work (Frank & Galinsky, 2016b) we emphasized that the EFD formulation was very general, being based on a Bayesian formulation of field theory that facilitated the incorporation of whatever prior information might be available for the problem at hand. The novel feature of the EFD approach is its use of the ESP method to derive the optimal (in the maximum entropy sense) parameter configurations when the space-time coupling structure of the data is used as that prior information. For multimodality data, the EFD method can be generalized by incorporating the coupling structures not only within the data but between modalities. Here we demonstrate this capability using the three major neuro-MRI modalities: HRA, DTI, and rsFMRI. Our results demonstrate an increased resolution and specificity of both structural and functional networks in the brain. This analysis allows us to uncover a scale-free organization of brain structural and functional networks that follow virtually identical power laws.

## 2 Theory

The EFD method developed for detection and analyses of nonlinear and nongaussian modes and patterns present in heterogeneous data sets is well suited to the problem of multimodal connectivity estimation. The cornerstone of EFD estimation procedure is the construction of local connectivity (including spatial connectivity, temporal connectivity, or both) present in the analyzed data set. The extension of the EFD approach to the joint processing of multimodal data sets requires a reformulation that combines and fuses these local connectivity patterns into a common coupling matrix. The quantitative EFD procedure can then be applied to this joint intermodality coupling matrix, thus allowing the detection and ranking of the connectivity present in those individual data sets, particularly emphasizing data sets with significant intermodal similarities.

### 2.1 Overview of EFD

### 2.2 The Coupling Matrix for Multiple Modalities

Constructing either linear or nonlinear solutions for the unknown signal requires an explicit generation of either the coupling matrix itself or the coupling coefficients . Using single-modality data sets as a source of coupling allows the generation of powerful and effective single-modality data analysis approaches; the coupling matrix obtained from the diffusion-weighted imaging data produced an accurate and effective way for simultaneous estimation of local diffusion and global tractography with the geometrical optics–based GO-ESP method. Similarly, the coupling matrix based on the correlation properties of time series from rsFMRI data sets resulted in accurate and effective characterization and estimation of space-time rsFMRI brain activation patterns.

This construction of the intermodality coupling involves several simple assumptions. Two neighboring voxels and are considered to be highly coupled if (1) their structural image intensities are similar, (2) their diffusion-weighted spin density functions (equation 8 in Galinsky & Frank, 2015) point to each other, and (3) their resting-state FMRI signals are maximally correlated across all frequencies. The mathematical details that incorporate and formalize these assumptions into coupling matrix expressions are given in appendix B.

### 2.3 The Structural, Functional, and Effective Connectivity Eigenmodes

The intermodality coupling matrix facilitates the development of a unified description of brain connectivity that includes both the structural dependencies (derived from long-range DWI connectivity and HRA morphology) and the functional correlations (based on activation modes buried within the functional MRI data).

One particularly informative way in which the intermodality coupling matrix can be used is to generate the effective streamlines that characterize both structural and functional dependencies present in the data by combining the correlation function for the functional signal with the structural information from the spin density function, expressing the joint functional-structural dependencies through effective transition probabilities (see equation A.5). This is the same procedure that was used to generate fiber tracts from DTI data in the original GO-ESP method and subsequently used to perform functional tractography using the EFD method, and now extended and generalized to a multimodal coupling matrix. The resulting streamlining therefore represent high-probability pathways in the joint structural-functional space.

## 3 Implementation

- •
The processing work flow starts with HRA data set. First, it extracts the brain, doing skull stripping, and then generates the SWD for the brain volume to be used in volume segmentation and registration later.

- •
For the DWI data set the GO-ESP procedure, which generates a new data volume of the equilibrium probabilities (EP) data, is used.

- •
For the rsFMRI data set, the EFD procedure is used to generate a new data volume of the mean field data (i.e., a volume that contains averages of all activations), removing both temporal and spatial outliers.

- •
The DWI equilibrium probabilities and the rsFMRI mean field volumes are then registered with the HRA brain volume using the symplectomorphic registration SYMREG-ESP, providing forward and inverse maps.

- •
The symplectomorphic maps, together with the HRA SWD volume, are then included in GO-ESP global tractography to generate high-resolution, very-high-density fiber tracts with the possibility of more than twofold improvement of the DWI spatial resolution and potentially resolving crossing fibers at below 10 angular resolution reported in Galinsky and Frank (2016a).

- •
The high-resolution, very-high-density fiber tracts distribution is used to generate ESP long-range structural connectivity eigenmodes.

- •
The symplectomorphic maps, together with the HRA SWD volume and the long-range structural connectivity eigenmodes, are used in rsFMRI EFD processing, providing greater improvement in resolved details of rsFMRI activation modes than is possible without the information provided by the other modalities by generating high-resolution effective connectivity eigenmodes. A fusion of structural and functional coupling provides an important advantage of those high-resolution effective connectivity eigenmodes over single-modality FMRI activation modes. The effective eigenmodes include both WM and GM regions in a way consistent with both structural and functional dependencies.

Additional modalities can easily be incorporated into this work flow.

## 4 Results

### 4.1 Data

To test our approach, we used one of the publicly accessible data sets available from the open source Human Connectome Project (Van Essen et al., 2012, 2013; Sotiropoulos et al., 2013), as well as several diffusion imaging data sets collected locally at the UCSD Center for Functional MRI (CFMRI).

The HCP data set was collected on the customized Siemens 3T Connectom scanner, a modified 3T Skyra system (MAGNETOM Skyra Siemens Healthcare), housed at the MGH/HST Athinoula A. Martinos Center for Biomedical Imaging (see Setsompop et al., 2013, for details of the scanner design and implementation). A 64-channel, tight-fitting brain array coil (Keil et al., 2013) was used for data acquisition. The data set contains 96 slices of 140 140 matrix (1.5 mm linear voxel size) at four levels of diffusion sensitizations (b-values = 1000, 3000, 5000, and 10,000 s/mm^{2}) distributed over 552 total q-vectors.

The CFMRI data were acquired with a 3T GE Discovery MR750 whole body system. The anatomical T1 volumes have 168 256 256 voxel size with 1.2 0.9375 0.9375 mm^{3} resolution.

A multiband DTI EPI acquisition (Setsompop et al., 2011) developed at the CFMRI employed three simultaneous slice excitations to acquire data with three diffusion sensitizations (at b-values b = 1000/2000/3000 s/mm^{2}) for 30, 45, and 65 different diffusion gradients (respectively) uniformly distributed over a unit sphere. Several baseline (b = 0) images were also recorded. The data were reconstructed offline using the CFMRI’s multiband reconstruction routines. The DWI data sets have 100 100 72 matrix size with 2 2 2 mm^{3} resolution. Whole brain BOLD resting-state data were acquired over 30 axial slices using an echo planar imaging (EPI) sequence (flip angle = , slice thickness = 4 mm, slice gap = 1 mm, FOV = 24 cm, TE = 30 ms, TR = 1.8 s, matrix size = ). (Further details are available in Wong, Olafsson, Tal, & Liu, 2013.) All data were preprocessed using the standard preprocessing analysis pathway at the CFMRI (as described in Wong et al., 2013). Nuisance terms were removed from the resting-state BOLD time series through multiple linear regression. These nuisance regressors included: (1) linear and quadratic trends, (2) six motion parameters estimated during image coregistration and their first derivatives, (3) RETROICOR (second-order Fourier series) (Glover, Li, & Ress, 2000) and RVHRCOR (Chang & Glover, 2009) physiological noise terms calculated from the cardiac and respiratory signals, and (4) the mean BOLD signals calculated from white matter and CSF regions and their first respective derivatives, where these regions were defined using partial volume thresholds of 0.99 for each tissue type and morphological erosion of two voxels in each direction to minimize partial voluming with gray matter.

### 4.2 Joint DWI and HRA Tractography

The data for every subject with T1 HRA, DWI and rsFMRI data sets were analyzed twice, first using independent processing of each modality and then the joint modality approach developed in this article. The addition of HRA coupling into the diffusion estimation and tractography process not only provides a significant boost to the overall level of detail that can be captured but, more important, shows features that are not apparent when the modalities are analyzed independently. Figure 2 shows details of the equilibrium probability map obtained from a single-modality DWI data set using several different views in panel a and compares them with the same views of the map produced from the joint HRA + DWI coupling matrix in panel b. The jointly estimated maps clearly show improvements in image quality: the twofold increase in linear image resolution provides more small-scale detail and greater noise suppression. The localization of white matter fiber bundles is clearly enhanced significantly in the jointly estimated maps.

These joint estimation improvements can also be seen in tractography results, shown in Figure 3. The unified DWI-HRA estimation procedure is again able to detect additional details in overall fiber tract distribution compared to the DWI procedure performed independently, as well as providing increased resolution of the tracts.

This improvement is even more pronounced if one constructs a volume generated from the tract counts in every voxel produced by GO-ESP whole brain tractography (shown in Figure 4). Both the HCP and CFMRI sets of tracts were produced using 5 million seeds. The joint DWI + HRA processing clearly provides significantly enhanced level of details, allowing reconstruction of subtle white and gray matter high-resolution spatial variations, while at the same time significantly increasing the directional information available from diffusion data.

### 4.3 Joint Structural, Functional, and Effective Eigenmodes

Using joint connectivity matrices allows us to generate eigenmodes that can be instrumental in finding and understanding the most important structural and functional dependency patterns present in the underlying data and can help to uncover and quantitatively characterize an effective information exchange that happens on different scales of the brain and modulates overall brain performance.

Figure 5 shows several structural connectivity eigenmodes that were produced by running DWI + HRA GO-ESP tractography and using streamlines to generate the eigenmode (or streamline) coupling matrix that includes connections between distant voxels (see equation 2.8). Because the streamline coupling matrix contains the number of connections between every pair of voxels located both inside the white or gray matter of the brain and on the cortical surface, it provides volumetric probability weighted connectivity that is significantly more informative than binary on-off connectivity typically generated for cortical regions. Furthermore, because the seeding for GO-ESP tractography uses equilibrium and transitional probabilities to sample spatial and angular distributions of tracs, it ensures consistency with both the diffusion and anatomical data volumes that were acquired. The eigenmodes of the streamline coupling matrix can be viewed as interconnected regions of a brain that are only loosely connected with other regions. The structure of these eigenmodes includes multiple modality constraints and, hence, can provide a more comprehensive description and ranking of the connectivity patterns reflecting the input from multiple modalities.

Selected functional nonlinear EFD eigenmodes obtained using FMRI + HRA connectivity matrix are shown in Figure 6b. Those modes are similar to EFD modes generated solely from rsFMRI data shown in Figure 6a and reported in Frank and Galinsky (2016b) but have more details and show better overall spatial resolution due to the addition of high-resolution anatomical details in FMRI EFD coupling.

Combining the DWI, FMRI, and HRA data sets through the joint coupling matrix, equation 2.4, produces effective connectivity eigenmodes (see Figure 6c that spread through the entire brain volume, covering both gray and white matter areas, hence complementing cortical activation regions available from FMRI with white matter structural pathways available from diffusion data.

We emphasize that in contrast to functional modes estimated from single-modality rsFMRI processing where the signal reflects BOLD activity happening mostly in the gray matter regions, the effective connectivity modes generated by multiple modalities incorporate not only gray matter regions but may also include white matter. This does not mean that our joint modality processing detects BOLD activity in the white matter. The extension of the modes into the WM is a manifestation of the fact that both functional and structural connectivity are related to a common origin—the brain neural structure—and often functional and structural dependencies may be deduced (although with limited degree of accuracy and reliability) one from the other (Skudlarski et al., 2008; Honey et al., 2009).

The common underlying nature of both structural and functional modes can be clearly seen when plotting mode power (or eigenvalue) as a function of mode number, as shown in Figure 7. The plot shown includes 200 structural eigenmodes and 180 functional EFD modes obtained from DWI + HRA and FMRI + HRA coupling matrices. Both structural and functional dependencies show power law decay with the same exponent that is close to 1.2. This type of power law dependence can typically be associated with a manifestation of scale–free organization of brain networks (Beggs & Timme, 2012).

The effective connectivity eigenmodes provide a simple and efficient way to construct a whole brain volumetric parcellation by assigning each voxel to a mode that has the largest contribution to this voxel. Two parcellation examples are shown in Figures 8 and 9. Figure 8 uses for analysis locally (UCSD CFMRI) acquired data sets, and Figure 9 is based on the publicly available HCP data volumes. Both figures were produced by generating high-resolution volumes (with the same resolutions as the correspondent HRA volumes: 168 256 256 for Figure 8 and 176 256 176 for Figure 9) that contain in each voxel the number of streamlines going through this voxel. The high-resolution volumes were used to extract the shapes shown in figures using different thresholds. The parcellation is constructed from the high-resolution effective connectivity eigenmodes.

## 5 Conclusion

In this article we presented a new theory for the joint estimation of the structural-functional brain modes that we have applied to the three primary neuro-MRI modalities: high-resolution anatomical (HRA) data, diffusion tensor imaging (DTI) data, and resting-state functional MRI (rsFMRI) data. The approach is rooted in our previous theoretical formulations of the SWD, GO-ESP, and EFD methods, which use our general theory of entropy spectrum pathways (ESP) (Frank & Galinsky, 2014), whose central concept is that relationships between local and global aspects of a lattice (e.g., a 3D image or a 4D time-dependent image) are characterized by the relationship, or coupling, between local elements (i.e., voxels). The SWD is a straightforward characterization of the intensity variations in HRA image volume in terms of a spherical wave decomposition (Galinsky & Frank, 2014). The GO-ESP technique uses the SWD to characterize the diffusion signal in each voxel from DWI acquisitions wherein a region of -space is sampled (Galinsky & Frank, 2015). These local (i.e., voxel) diffusion profiles are then used to construct couple matrices that, through ESP theory, provide optimal spatial pathways that can be used to drive the geometric optics–based tractography. In the EFD method for FMRI, spatial-temporal coupling matrices from the 4D data are used to generate optimal space-time paths from which brain networks are constructed (Frank & Galinsky, 2016a). In this article, we have developed a general theory of structural-functional brain modes that integrates these methods through a generalization of the coupling structure of the data to incorporate the different modalities and then modifying the individual algorithms to incorporate the coupling between modalities.

Our computational implementation demonstrates a significant improvement of mode reconstruction accuracy and new levels of detail. The fiber pathway reconstruction of joint modality processing offers a more than twofold improvement of the DWI spatial resolution and the ability to resolve crossing fibers at or below angular resolution (Galinsky & Frank, 2016a). The resting-state brain network modes, eigenmodes, and tractography generated provide more than four times the improvement in spatial resolution. Moreover, this joint estimation procedure allowed us to demonstrate the novel finding of a scale-free organization of brain structural and functional networks that follow virtually identical power laws. Further exploration of this finding may help to uncover important structure-function relations in human brain development.

In addition to providing an overall increase in resolution, accuracy, level of detail, and information content that can provide valuable insight into structure-function relations of interest in basic neuroimaging studies, this method also holds the potential to enhance the clinical utility of these individual neuro-MRI methods. The current procedure of analyzing these data separately poses a significant challenge for practicing clinicians faced with the task of visually integrating the information from these complex image acquisitions. This new method provides a robust and accurate method for integrating these data to provide a more accurate assessment of brain structure and function that is nonetheless amenable to clear visualization and statistical assessment via nonlinear registration critical to the development of viable clinical protocols.

### Appendix A: EFD Summary

When the source term , the linear information propagator , and the nonlinear interaction terms are all known, or at least some more or less accurate approximations can be used for their description, the IFT approach provides an effective, powerful, and mathematically elegant way to find an unknown signal either by using the classical solution at the minimum of Hamiltonian () or with the help of summation methods (e.g., with the help of Feynman diagrams; Feynman, 1949; Ensslin et al., 2009).

### Appendix B: Coupling Matrices

To investigate long-range structural, functional, and effective connection patterns that span both gray and white matter regions, we generate the connectivity matrix using structural, functional, and effective streamlines produced by GO-ESP tractography. For structural tractography, we employ the nearest-neighbor connectivity matrix generated from DWI and HRA data sets and follow the same formalism as was used for single-modality DWI tractography.

## Acknowledgments

L.R.F and V.L.G. were supported by NSF grants DBI-1143389, DBI-1147260, EF-0850369, and PHY-1201238 and NIH grant R01 MH096100. We thank Alec Wong and Tom Liu at the UCSD CFMRI and Scott Sorg at the VA San Diego Health Care System for providing the diffusion-weighted and the resting-state imaging data. Data were provided in part by the Human Connectome Project, WU-Minn Consortium (principal investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH institutes and centers that support the NIH Blueprint for Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at Washington University. Data collection and sharing for this project was provided by the Human Connectome Project (HCP; principal investigators: Bruce Rosen, Arthur W. Toga, and Van J. Weeden). HCP funding was provided by the National Institute of Dental and Craniofacial Research, the National Institute of Mental Health, and the National Institute of Neurological Disorders and Stroke. HCP data are disseminated by the Laboratory of Neuro Imaging at the University of Southern California.