Functional ultrasound imaging of intrinsic connectivity in the living rat brain with high spatiotemporal resolution

Long-range coherences in spontaneous brain activity reflect functional connectivity. Here we propose a novel, highly resolved connectivity mapping approach, using ultrafast functional ultrasound (fUS), which enables imaging of cerebral microvascular haemodynamics deep in the anaesthetized rodent brain, through a large thinned-skull cranial window, with pixel dimensions of 100 μm × 100 μm in-plane. The millisecond-range temporal resolution allows unambiguous cancellation of low-frequency cardio-respiratory noise. Both seed-based and singular value decomposition analysis of spatial coherences in the low-frequency (<0.1 Hz) spontaneous fUS signal fluctuations reproducibly report, at different coronal planes, overlapping high-contrast, intrinsic functional connectivity patterns. These patterns are similar to major functional networks described in humans by resting-state fMRI, such as the lateral task-dependent network putatively anticorrelated with the midline default-mode network. These results introduce fUS as a powerful novel neuroimaging method, which could be extended to portable systems for three-dimensional functional connectivity imaging in awake and freely moving rodents.

T he brain dynamically integrates and coordinates responses to internal and external stimuli across multiple spatiotemporal scales through large-scale functional networks. Assessment of its functional connectivity (FC), through the measurement of regionally correlated, spontaneous, lowfrequency (0.01-0.1 Hz) fluctuations in blood oxygen leveldependent (BOLD) signals with functional magnetic resonance imaging (fMRI), particularly during resting-state/task-free periods (resting-state fMRI or rsfMRI) 1 , has greatly advanced our understanding of the functional organization of the human brain 2 . Intrinsic connectivity networks, such as the default-mode network 3,4 , ventral and dorsal attention networks 5,6 , and salience network 7 , have been intensely studied in both basic and clinical cognitive neuroscience fields. These correlated resting-state BOLD fluctuations appear to be a fundamental property of the brain because they are present during sleep 8 and even during general anaesthesia 9 . Indeed, these fluctuations are temporally coherent among brain areas that are structurally connected and functionally related 1,10,11 . Most, if not all, neurological and psychiatric diseases involve the disruption of large-scale functional and/or structural properties in the brain 11,12 . These include major pathologies such as schizophrenia, depression and Alzheimer's disease. Consequently, investigation of the FC under well-controlled and experimentally accessible conditions is of major scientific interest, because it may lead to better diagnostic or prognostic indicators, and more targeted and controlled drug treatments. Therefore, the development of the corresponding translational rodent models is also of great interest. Indeed, recent reports show the presence of prominent intrinsic networks in the monkey 9 and rat brain 13 .
To date, only rsfMRI has been able to image intrinsic brain networks with the appropriate spatial resolution and coverage. However, the low-frequency oscillations measured in rsfMRI studies can be contaminated by higher frequency (41 Hz) physiological noise, such as the cardiac cycle and respiratory motion 14,15 . Therefore, the unambiguous and general exclusion of physiological noise requires novel techniques that are able to temporally resolve signals above 1 Hz. Magnetic resonance imaging (MRI) has additional technical and logistic challenges, such as the necessity for high magnetic fields, problems with motion artefacts, electromagnetic compatibility issues, high costs and lack of portability. These concerns currently hinder the broad dissemination of FC research in translational and pre-clinical research settings. Therefore, the validation of complementary or alternative methods for in vivo imaging of intrinsic FC is an important scientific objective. Recently, very high frame rate ultrasound imaging (410,000 frames per second) was shown to enable high-resolution and high-sensitivity power Doppler imaging of cerebral blood flow (CBF) 16 , leading to functional ultrasound (fUS) imaging of task-evoked changes in cortical activity 17 . However, this study used highly invasive craniotomy, which is not optimal for the imaging of resting-state functional networks. Indeed, craniotomy was shown to induce important behavioural, morphological, biochemical and vascular changes in the rat brain [18][19][20] . In consequence, currently we do not know whether fUS is capable of mapping intrinsic FC patterns and we do not know the relative performance of this method as compared with the standard fMRI-based approach.
Here we develop a large thinned-skull cranial window preparation and an ultrasonic sequence for ultrafast imaging of intrinsic FC of the living adult rat brain. This approach enables us to identify highly contrasted intrinsic FC patterns. Our approach represents a simple, portable, cost-and space-effective, highly resolved, and motion and pulsatility artefact-free method, capable of imaging FC deep in the rodent brain.

Results
Identification of putatively interconnected seed areas. We have previously shown that fUS allows measurement of functional haemodynamic changes induced by sensory stimulation 17 . To test whether fUS is capable of detecting correlated patterns, we aimed to measure spontaneous haemodynamic fluctuations between contralateral homologous cortical areas that share similar function and are massively interconnected by axonal projections through midline commissural structures such as the corpus callosum 21 . To achieve this, we first identified functionally correlated contralateral cortical areas activated by electrical stimulation of the right or left sciatic nerve. Preliminary mapping of the rat brain vasculature (Fig. 1a-d) allowed recognition of the blood vessel architecture of large-scale brain regions on the coronal level where the fUS imaging was performed. Figure 1e,f show the vascularization of the rat brain detected with Ultrafast Doppler imaging, over which we superimposed the evoked fUS signal, measured using the fUS sequence number 1. We have denoted regions that correspond to the rat brain atlas of Paxinos 22 . Accordingly, five consecutive electrical stimulations induced a robust and reproducible increase in haemodynamic response in the contralateral primary sensory cortex, hindlimb part (S1HL) and primary motor cortex (M1) regions (Fig. 1g).
fUS detection of FC. Next, we asked whether functionally and anatomically connected regions show correlated changes in spontaneous fUS signal fluctuations. Once the bilateral S1HL and M1 regions were functionally identified (Fig. 2a,b), the spontaneous fUS signal of this coronal slice was measured, without stimulation, for 10 min. Figure 2c compares the temporal variation of the spontaneous power Doppler signals of the two 'seed regions', respectively, identified based on the previously evoked blood flow response of the right (blue delimitation) and left (red delimitation) sciatic nerve stimulations (Fig. 2b). Strikingly, the fUS signals from these two contralateral, but functionally similar regions were highly correlated (correlation coefficient r40.8). The correlation function (blue curve in Fig. 2f) also indicated that the two signals were in phase, as the maximum of the correlation function occured at a zero lag time point. Notably, not all cortical regions were similarly correlated. Figure 2d compares the signal from the left seed (red delimitation in Fig. 2b) to a signal extracted from the left anterior secondary cingular cortex (green delimitation in Fig. 2b). These signals were less correlated than the previously analysed contralateral regions (correlation coefficient r ¼ 0.25). The low-amplitude fluctuations of both correlation functions in Fig. 2f display no secondary maximum, indicating that the fluctuations of the different time signals were random. Analysis of the power spectral density of the signal of the left seed (red delimitation in Fig. 2b) revealed that its temporal variations were of low frequency (90% of the total power spectral density is below 0.15 Hz, Fig. 2e). Mapping the cortical regions where the signal had an elevated correlation coefficient with the left or right S1HL and M1 seed regions enabled labelling of the bilateral, motor and sensory cortical areas (Fig. 2g,h) with high inter-animal reproducibility (N ¼ 6, Fig. 3). The mean Pearson correlation coefficient for both right and left 'seed-based' FC maps was 0.65±0.04 (Po0.001, N ¼ 6).
In conclusion, our results demonstrated that the fUS-measured, spontaneous low-frequency fluctuations of the CBF were highly correlated in the sensorimotor system. The frequency range and frequency distribution of the correlated spontaneous CBF fluctuations are in accordance with values reported by previous BOLD-based studies 1,10,23,24 .
Seed-based determination of FC at different coronal planes. In addition to the coherent activity of the sensorimotor system reported above, BOLD-based studies have reported patterns of intrinsic FC in many other neuroanatomical systems, including visual 25 , auditory 26 , hippocampal 26 , dorsal attention 25 and ventral attention systems 5,6,25 . To investigate whether coherent fUS measurements would reveal similar spatial patterns, we assayed the FC in three different coronal planes, based on calculation of the mean Pearson correlation factor in the fUS signal between the anatomically defined regions of interest (ROI). These ROIs were obtained through overlaying the power Doppler vascular map with the spatial referential frame of the Paxinos atlas 22 . The results are presented either as correlation matrices (Fig. 4a,b) or as the projection of these correlation factors on the schematic atlas view (Fig. 4c). The correlation coefficient matrices were highly reproducible between individual animals. The mean Pearson correlation coefficients of the FC matrix intercorrelation coefficient, using only the non-diagonal values, were 0.85 ± 0.03 (Po0.001, N ¼ 6) at Bregma À 0.6 mm, 0.86 ± 0.04 (Po0.001, N ¼ 4) at Bregma þ 0.84 mm and 0.86±0.04 (Po0.001, N ¼ 5) at Bregma À 2.16 mm. The mean correlation matrix of each coronal plane is displayed in Figs 4a and 5c,d. All of the matrix coefficients that were 40.2 were reproducible (Po0.05). Only the reproducible coefficient values (Po0.05) are displayed in these figures.
Furthermore, the brain areas of the 'sensory-motor restingstate network' were overall highly connected to each other (for instance, the S1HL with the motor cortex, Fig. 4a), as previously described 27 . However, the cingulate cortex (both primary and secondary areas) was not correlated with other cortical structures, except the neighbouring primary motor cortex (Fig. 4a), as previously described 28 . By contrast, cerebral structures such as the CPu, septum, thalamus and hippocampus showed no correlation with the 'sensory-motor resting-state network' (Figs 4a,c and 5c) and little or no correlation with each other. These observations are consistent with the previous description of their respective implications in segregated pallidumlike (septum), thalamic (thalamus) and retrohippocampal (hippocampus) resting-state networks 27 .
In addition to reproducing the spatial patterns previously reported by fMRI BOLD imaging, the high sensitivity and high spatial resolution (100 mm Â 100 mm in-plane pixel size) of fUS imaging allowed the identification of additional specific connectivity patterns between parts of the above brain areas, such as the preferential connection of the cingulate cortex with the primary, but not secondary, motor cortex and a strong connection between different parts of the primary sensory cortex, such as the S1HL and S1FL (Fig. 4a,c). In conclusion, fUS imaging was able to detect FC patterns of distinct neuroanatomical systems with 100 mm spatial resolution.
Filtering of high-frequency physiological noise. The cerebral blood volume fluctuates with the cardiac cycle (around 5 Hz), which could lead to biased measurements of the low-frequency components, a well-known concern for fMRI BOLD imaging 29,30 . Therefore, to precisely evaluate the components of the spontaneous fUS signal fluctuation patterns, we continuously measured the fUS signal in the neocortex for 500 s at our standard 500 Hz sample rate (fUS sequence number 2). This high frame rate, which was chosen to adequately sample the Doppler ultrasound signal without aliasing 31 , oversamples the fUS variations resulting from both intrinsic FC and cardiac pulsatility. Indeed, as shown by the red curve in Fig. 7b,c, which was computed from the area boxed in red in Fig. 7a, both low-and high-frequency fluctuations were observed. The power spectral density revealed two peaks: one low frequency (o1 Hz) and a second one at B7 Hz (Fig. 7d,e), which corresponded to the pulsatility of the CBF. The low-frequency components were grouped around 0.1 Hz (Fig. 7e), the frequency range of the intrinsic FC 1 . This demonstrates that the blood pulsatility signal is not negligible compared with the intrinsic FC signal.
As the pulsatility was correctly sampled with fUS sequence number 2, accurate filtering of the pulsatility could be performed.  1 | Functional identification of putatively interconnected seed regions using task-evoked haemodynamic changes. Anatomical organization of the rat brain at Bregma À 0.6 mm (a) and mapping of the cerebral blood vessels at the same level using either staining with DiI (b) or India ink (c,d). (c,d) Black staining: India ink; blue counter staining: toluidine blue. Between the two techniques for vasculature staining, DiI is more sensitive. However, India ink more clearly shows the typical curved vascular staining observed in the choroid plexus, which is prominently displayed in the fUS image. (d) High-power magnification of c. Functional ultrasound (fUS) imaging evoked in the left (e) or right (f) somatosensory cortex, hindlimb part (S1HL) using electrical stimulation (5 Hz, 0.2 mA, 100 ms width for 10 s, separated by 20 s) of the right (e) or left (f) sciatic nerve applied either on the right (e) or left (f) side, respectively. (g) Time course changes in the evoked fUS response in the left S1HL following stimulation of the contralateral sciatic nerve. Black bar: duration of the stimulation. Electrical stimulation of the left hindpaw induced a reproducibly increased fUS signal in the contralateral (right) S1HL. Scale bars, 2.3 mm (a,b), 1.5 mm (c), 375 mm (d), 1 mm (c-f).
Indeed, low-pass filtering (butter third order, 0.5 Hz frequency cutoff) resulted in the complete removal of the high-frequency fluctuations (Fig. 7f), without changing the low-frequency components (Fig. 7g), generating the blue curve in Fig. 7b,c.
Because of the practical limitations in computing bandwidth (see Methods), the continuous imaging sequence number 2 is currently not adapted for real-time analysis of large brain areas. However, the cardiac cycle was already fully sampled with the intermittent (2 s) low-bandwidth sampling sequence number 1 (500 Hz frame rate with a sampling time of 400 ms, corresponding to two cardiac cycles). Therefore, the mean value of this signal acted as an efficient low-pass filter that removed most of the pulsatility signal. Accordingly, the results obtained by using the sequence number 2 were fully consistent with the results given by sequence number 1 (Fig. 7h). The Pearson correlation coefficient computed between the two FC matrices using only the non- Previously determined seed regions, that is, the right (blue) and left (red) S1HL þ M1. (c) Spontaneous variations in the power Doppler signal in the seed regions marked by the evoked blood flow response of the right (blue curve) and left (red curve) S1HL show high temporal correlations. The curves shown are typical in terms of fluctuation amplitude. (d) Example of temporal variations in signals that are weakly correlated in the S1HL þ M1 and the ipsilateral secondary cingulate cortex. (e) Frequency distribution of the power spectral density in the seed regions (S1HL þ M1, in red in b). (f) Temporal correlation function between the signals obtained at the seed region (left S1HL þ M1, red delimitation in b) and the contralateral S1HL þ M1 region (blue curve, corresponding to the zone delimited in blue in b) or between the same seed region and the ipsilateral secondary cingulate cortex (green curve, corresponding to the zone delimited in green in b). Averaged spatial pattern of cortical regions temporally correlated to the right (g) or left (h) S1HL þ M1 seed regions (N ¼ 6). Scale bars, 1.1 mm (b), 1 mm (g,h). ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms6023 diagonal coefficients of the two different acquisition types was r ¼ 0.97, establishing that the intermittent sequential acquisition sequence number 2 was not biased with pulsatility artefacts.
In conclusion, both of the fUS acquisition sequences efficiently filtered out the cardiac pulsatility noise. Consequently, all coherent spatial patterns of intrinsic connectivity were based on signals of low-frequency changes specific to the modulation of CBF.
Data-based identification of anticorrelated FC patterns. Interestingly, BOLD-based connectivity studies have suggested that regions with apparently opposing functionality display temporally anticorrelated spontaneous signals, both in human and rodent brains 28,32,33 . Specifically, externally focused or task-related networks are anticorrelated with the intrinsic default-mode network 28,32 . In BOLD-based studies, pre-processing removal of both neuronal and non-neuronal global noise is often required to visualize these patterns, but the potential introduction of artefactual coherence patterns is a subject of ongoing debate 30,34-36 . As fUS-based measurements of intrinsic connectivity are free of high-frequency physiological noise (see  ARTICLE above), we investigated in detail the connectivity patterns in rodent brains, without global noise extraction, to reveal these complex intersystem relationships. For this purpose, we used singular value decomposition (SVD) as a linear decomposition operator, which provided powerful data decomposition without any spatial or temporal information a priori. After SVD-based noise removal (see Methods), we identified the global spatial modes (GSMs) of the SVD decomposition, which represent the mean (that is, global) FC for the entire population of N animals at each of the three investigated coronal planes ( Fig. 8).
At least five or six GSMs were above the noise level at each coronal plane (at Bregma þ 0.84, þ 0.6 and À 2.16 mm), reproducibly revealing highly contrasted FC networks (Fig. 9a,b). The first GSM, corresponding to the most prominent connectivity pattern at each coronal plane, was dominated by cortical regions including the sensory, motor, frontal, cingulate and retrosplenial cortices (Fig. 9a), indicating strong, direct bilateral communication across the cortical ribbon 24 . Notably, this pattern did not include the secondary anterior cingulate cortex, a region that showed a strong correlation with the primary anterior cingulate cortex using seed-based analysis (Fig. 4a). Interestingly, the second most prominent connectivity pattern (GSM 2) at each coronal plane identified the midline-structure anterior cingulate and retrosplenial cortices, and the dorsal hippocampus at Bregma À 2.16 mm. These structures were negatively correlated with the more lateral parts of the cortical ribbon of GSM 1 on all three coronal planes, with the exception of the most medial part of the secondary motor cortex (M2). The third GSM at the two most rostral planes revealed a positive correlation of the contralateral anterior cingulate cortices and the dorsal part of the CPu, and these networks were also negatively correlated with the more lateral motor and somatosensory cortices (Fig. 9a). Finally, the fifth GSM, at Bregma þ 0.84 mm, and the fourth GSM at Bregma þ 0.6 and À 2.16 mm, showed a strong anticorrelation of the midline anterior cingulate and retrosplenial cortices with the more lateral motor cortices (M1 and M2), as well as a strong bilateral connection of the subcortical regions, such as the CPu at Bregma 0.84 and À 0.6 mm (Fig. 9a), and the dorsal hippocampus and thalamus at Bregma À 2.16 mm (Fig. 9a).
In conclusion, SVD analysis of fUS connectivity correlations separates, at high spatial resolution, several overlapping and previously unreported FC patterns that correspond to anatomically well-defined structural connectivity networks.

Discussion
Coordinated spontaneous brain activity, indicating functional brain connectivity 1 , has received considerable attention in recent years. It has been investigated with multiple modalities, CPu Septu m S1FL predominantly based on fMRI 1,37 , but also with other functional imaging techniques such as optical imaging 38 , positron emission tomography 3 , magnetoencephalography and electroencephalography 39,40 . The data obtained through these different modalities, each with its own specificity, sensitivity and spatiotemporal resolution, represent different types of physiological signals and are affected by different artefacts 41 . In this study, by using a large thinned-skull cranial window preparation and a novel ultrasonic sequence for ultrafast imaging of intrinsic FC of the living adult rat brain, we showed that fUS is highly efficient for mapping FC in a rodent model. This introduces a new neuroimaging modality for FC research with several novel key features. First, we demonstrated that we were able to detect, both using seed-based and SVD analysis, strong and highly contrasting spatial coherence signals in low-frequency (o0.1 Hz) spontaneous fUS signal fluctuations. The resulting intrinsic FC patterns were similar to known major functional networks described in humans by using resting-state fMRI, such as the lateral taskdependent network putatively anticorrelated with the midline default-mode network. This type of anticorrelated activity is an important feature of dynamic brain organization 32 , putatively important for efficient cognitive function 42,43 , and has clinical predictive power in the treatment of depression 44 . However, preprocessing removal of global noise in fMRI-based studies potentially introduces artefactual coherence patterns and is a subject of ongoing debate 30,[34][35][36]41,45,46 . Here the high temporal resolution of the fUS imaging coupled with the SVD analysis enabled us to show, without preprocessing removal of neuronal or non-neuronal global noise, the prominent presence of highly contrasting and anticorrelated connectivity patterns, which are similar to functional networks previously described in humans. Current developments are now underway to apply fUS to transcranial imaging in mice, through increasing the ultrasonic frequency range, to exploit the genetic tools available in this model organism.
Second, the reproducibility and robustness of this ultrasound method was demonstrated by comparing the seed-or SVD-based FC maps obtained from different animals. Notably, the mean Pearson intercorrelation coefficient of the seed-based FC matrix using only the non-diagonal values was at least 0.85±0.03 (at Bregma À 0.6 mm) and was typically statistically significant (Po0.001) on all coronal planes. Similarly, for the three coronal planes investigated in this study, SVD analysis reproducibly identified at least five distinct mean GSMs (N ¼ 6) above the noise level. Third, fUS imaging is able to spatially resolve changes in cerebral haemodynamics at 100 mm in-plane, which is higher than the typical 400 mm in-plane voxel dimension in recent fMRI studies of the rat brain 13,28 . This fourfold increase in spatial resolution considerably multiplies (4 2 ¼ 16-fold per slice) the amount of potentially available connectivity data. In addition to a more precise and detailed mapping of the FC networks, this may provide better imaging sensitivity at relatively small voxel volumes because of the elevated possibility of spatial filtering, a method that most investigators apply to improve the functional signal-to-noise ratio (SNR) of BOLD measurements 47 . In the near future, the extension of fUS to three-dimensional dynamic imaging using two-dimensional matrix array technology could lead to a 4 3 ¼ 64 times increase per volume in the number of available FC data compared with fMRI.
Fourth, the ultrafast 2 ms time resolution of fUS enables unambiguous discrimination between the resting-state activity and pulsatility artefacts before any SVD or independent component analysis (ICA) processing. This is important because the issue of cardiac and respiratory motion remains a potential concern in fMRI acquisitions 30 . Indeed, cardiac and respiratory cycles occur at around 6 and 1 Hz, respectively. Consequently, these signals become strongly aliased at typical repetition times (B1.7-5 s) in small animal fMRI, causing significant noise at frequencies typical of resting-state networks. Although De Luca et al. 29 showed that probabilistic ICA can partly solve this issue, and Chang and Glover 46 proposed a model-based noise-removal approach in human brain imaging, the usability of these approaches has not yet been demonstrated in rodent fMRI, where the temporal constraints for cardiac and respiratory motion are much higher.
How closely does the fUS signal reflect neuronal activity, and how does this compare with the sensitivity and fidelity of the BOLD signal, the gold standard method for FC studies? In the current models of functional brain imaging, neural activity produces complex local changes in the CBF, the cerebral metabolic rate of oxygen (CMRO 2 ) and cerebral blood volume following task-activation and during resting state (recent review by Liu 48 ). In the BOLD signal model, neural activation-induced fractional BOLD signal change is related to underlying changes in the CBF and CMRO 2 (refs 49,50). Because of the limited diffusibility of oxygen from the blood to the brain, the activationinduced increases in CBF largely overcompensate that of CMRO 2 . The difference between these two quantities yields a relatively moderate (compared with the CBF change in amplitude) reduction in the deoxyhaemoglobin content of the brain microvasculature, elevating the local magnetic resonance signal 49,50 . Consequently, the major parameter that determines the change in the BOLD signal, similarly to fUS, is the change in CBF. However, the amplitude of the BOLD signal (and consequently its SNR) is reduced by the parallel increase in CMRO 2 . Therefore, the direct measurement of local changes in CBF using fUS is a valid readout of localized neuronal activity, which is also a potentially more sensitive method than BOLD. In addition, by directly measuring the changes in the microvascular haemodynamics, fUS measurements are also free of another potential fMRI confound-blood CO 2 concentration, which can change with the respiration rate, leading to BOLD signal changes that are unrelated to neural activity 30 . Importantly, previous flow-sensitive MRI imaging in human subjects revealed coherent fluctuations in the CBF with similar spatial patterns to those seen with BOLD, both after task activation and during the resting state [51][52][53][54][55] . Based on these premises, the CBF was already proposed as a biologically more specific correlate of neural activity compared with the BOLD contrast image [56][57][58] , because of the more direct link between CBF and neuronal activity, as discussed above. Indeed, change in the CBF has the potential for a more accurate estimation of the location, magnitude and longitudinal variation of neural function, being more linearly related to changes in neural activity than BOLD 59 and having decreased intersubject and intersession variability, which permit a smaller sample size for a given effect size 56,57,60 . However, current MRI-based measurements of CBF have significantly lower sensitivity and temporal resolution than BOLD 58 . The fUS approach presented in our study overcomes these limitations and provides a novel perspective in FC imaging through direct, sensitive and highly resolved measurements of a biologically well-defined physiological quantity, CBF.
In present, constraint to three-dimensional (two dimensions in space plus time) imaging is a limitation of this study compared with fMRI. The extension to four dimensions will require the use of matrix arrays coupled to ultrasonic platforms with a great number of channels embedded and adequate graphic processing unit computational power. Nevertheless, such exponential developments are envisioned in the near future. Of particular interest, the number of achievable raw data sets used in fUS should reach a 64-fold increase (typically 100 mm versus 400 mm) in space and 1,000-fold increase in time (typically 1 ms versus 1 s) compared with fMRI. This will lead to a fourfold increase in the available data compared with current state-of-the-art imaging modalities. The impact of this increase on the SNR and robustness of data processing algorithms for FC studies should be carefully studied in future works.
In comparison with fUS, optical techniques such as optical intrinsic signal imaging can reach higher spatial (B1 mm) and comparable temporal resolutions (1 ms) than fUS by using a similar thinned-skull approach on rats, but it remains limited to the cortex surface due to the elevated scattering of light in living tissue. To increase the imaging depth, two strategies were developed. First, approaches based on fNIRS (functional nearinfrared spectroscopy) or diffuse optical imaging can be used to This matrix is highly similar to the connectivity matrix obtained using the intermittent frame rate, low-bandwidth sampling sequence number 1, shown previously in Fig. 4a. The mean Pearson correlation coefficient of the FC matrix intercorrelation coefficient, using only the non-diagonal values, was 0.80 ± 0.05, Po0.001, n ¼ 5.
recover deeper brain images, however with significantly lower spatial resolution (41 cm) 61 . Second, new multi-wave approaches such as photo-acoustics propose to solve this problem by taking benefit of the interaction between acoustic and optical waves, getting optical contrast and acoustical resolution in deep tissues. However, these approaches require the complex use of ultrasonic scanners coupled with laser sources 62 .
Compared with optical approaches, the major advantage of fUS consists in providing whole-brain functional imaging with a very good spatial and temporal resolution. Moreover, the fact that fUS relies solely on the use of ultrasonic waves at ultrafast frame rates enables to design very small imaging probes that can be implanted in the skull, paving the way to functional brain imaging of awake and freely moving animals. Finally, ultrafast ultrasound imaging (the core of resting-state fUS) can provide several other valuable information such as vascular resistivity 63 , tissue strain imaging and tissue pulsatility.
The thinned-skull imaging window size is also a putative limitation that impeded, in the present study, the imaging of laterally located areas, such as the amygdala, rostral structures such as the olfactory bulb, and caudal structures such as the brain stem. However, as deeper lying brain structures, such as the ventral parts of the CPu (see Fig. 8a), were well resolved, the experimental protocols may be easily adapted for fUS imaging of almost all lateral, rostral and caudal brain areas.
The prospect of resting-state fUS are exciting because the current development of light, portable ultrasonic probes will soon allow FC mapping in awake and freely moving animals, without motion artefacts that may perturb fMRI studies, even in human subjects 30 . Although recent studies have shown that fMRI-based FC imaging is possible in restrained, awake rats after training 24,64 , important experimental protocols, such as the task-related deactivation protocol that is cardinal for the definition of the default-mode network 65 , have not yet been realized because of the requirement for tight physical restraint during BOLD imaging.
Finally, the translation of resting-state fUS imaging to clinical applications is very promising. Although passing through the skull bone with ultrasonic beams remains a challenge, two clinical applications of major interest are not limited by this restriction Thresholding algorithm ARTICLE and are currently under development. First, non-invasive fUS imaging can be performed straightforwardly in newborns through the fontanel or the temporal window. Second, fUS imaging can be performed peroperatively during surgery. These clinical applications are of major interest, as the use of fMRI is extremely difficult in these practical configurations. Consequently, restingstate fUS imaging may become complementary to resting-state fMRI. In particular, the ability of fUS to perform functional imaging of newborn brain activity with a portable bedside technology is a major advantage compared with fMRI. Therefore, recent ultrafast Doppler mapping of cerebral vasculature in neonates paves the way to resting-state fUS in paediatrics 63 . After more than 25 years of functional neuroimaging, fMRI, positron emission tomography and electroencephalography/ magnetoencephalography are well-established techniques. Our study demonstrates that fUS is a novel, additional neuroimaging modality, which enables FC mapping in rodents with unprecedented resolution and sensitivity. Preparation of the large thinned-skull imaging window. To perform ultrasound imaging through the skull of adult rats, the skull was thinned to 50-100 mm over an area of B1.2 cm Â 0.9 cm at 1-7 days before imaging. This window dimension allows to use the central 128 elements of the transducer (0.08 mm per element Â 128 ¼ 10.24 cm large). Under anaesthesia (intraperitoneal (IP) injection of medetomidine (Domitor, 0.3 mg kg À 1 ) and ketamine (Imalgène, 40 mg kg À 1 )), the head of the animal was placed in a stereotaxic frame and the three layers of bone were consecutively removed by drilling (Foredom, USA) at low speed using a micro drill steel burr (Fine Science Tools, catalogue number 19007-07). During the thinning procedure, which typically lasted 90 min per rat for a skilled experimenter, the skull was frequently cooled with saline and an airstream, as suggested previously 66 , resulting in a lack of heating, swelling or oedema of the cerebral cortex. The thinned window was protected by a small (1 cm Â 1 cm) plastic cover and the skin was sutured using 5.0 non-absorbable Ethicon thread. Preliminary experiments showed that this method enabled good quality ultrasound imaging for as long as 1 week after skull thinning. The optimal results were obtained 24 h to 3 days after the preparation, as the bone tends to re-grow and the size of the ultrasound transparent window reduces with time, especially on the lateral sides.

Methods
Animal preparation for ultrasound imaging. Animals were anaesthetized using an initial IP injection of medetomidine (Domitor, 0.3 mg kg À 1 ) and ketamine (Imalgène, 40 mg kg À 1 ), followed by hourly IP injections of medetomidine (0.1 mg kg À 1 h À 1 ) and ketamine (12.5 mg kg À 1 h À 1 ). Imaging sessions lasted 4-6 h. This dose was chosen according to previously published data on the observation of resting-state FC in rats 67,68 . In addition, we performed a preliminary experiment where we compared the use of medetomidine alone versus a combination of medetomidine and ketamine. We measured stable heart and respiratory rates, similar stimulation-evoked haemodynamic responses and intrinsic connectivity patterns with both methods. However, to minimize the possible pain and discomfort that the animal may be subjected to during the electrical stimulations of the sciatic nerve, we chose a mix of medetomidine and ketamine because of the analgesic properties of the ketamine. We observed that the animals needed to be under this anaesthesia for 2 h before the acquisition of intrinsic connectivity patterns to have stable and reproducible results. For fUS imaging, the animals were placed in a stereotaxic frame (Stoelting; Chicago, IL, USA). Body temperature was maintained constant using a heating pad (Gaymar Industries, New York, NY, USA). For animals that received electrical stimulations of the sciatic nerve, a 1-cm incision was made in the skin and muscle of both hindpaws, above the femur. The sciatic nerve was exposed and isolated. A hook stimulating electrode was gently inserted to stimulate the sciatic nerve. The nerve was allowed to rest for 10 min before the beginning of the stimulation. Trains of five C fibre stimulations (5 Hz, 0.2 mA and 100 ms width for 10 s), separated by 20 s, were performed on each sciatic nerve, with coupled fUS imaging.
Staining of the cerebral vascular architecture. Ten naive Sprague-Dawley rats were anaesthetized by IP injection of sodium pentobarbital (80 mg kg À 1 ). Once the animal was deeply anaesthetized, a thoracotomy was quickly performed and an incision was made in the right atrium. For DiI staining (N ¼ 5 rats), 2 ml of saline, followed by 15 ml of DiI 69 and then 10 ml of 4% paraformaldehyde were perfused intracardially at a rate of 7 ml min À 1 . The brain was removed and fixed for 2 days in 4% paraformaldehyde. Alternatively, (for N ¼ 5 rats), 5 ml of India ink was injected intracardially in a 1-min bolus. The brain was then immediately removed and fixed for 1 month in 4% paraformaldehyde. The brains were cryoprotected in 30% sucrose for two days. They were then frozen in cooled isopentane and sectioned (50 mm sections) using a cryostat (Microm Microtech France, Rhône, France). Free-floating sections were collected in 0.02 M phosphate buffer. Half of the sections were immediately mounted on Superfrost slides and mosaic images were acquired using a Â 5 objective with an Axio Imager M1 microscope (Zeiss, Jena, Germany). The other half of the sections were mounted on Superfrost slides, air dried and counter-stained with 5% toluidine blue for 2 min. Finally, after two washes in distilled water, sections were dehydrated for 2 Â 5 min in 50% ethanol, 2 Â 5 min in 70% ethanol, 2 Â 5 min in absolute ethanol and 2 Â 5 min in xylene. Once dehydrated, coverslips were mounted with DPX (Sigma Aldrich, St Louis, MO, USA).
Imaging the rat brain with fUS. Two hours after the induction of anaesthesia, the sutures and protective plastic cover were removed. The thinned skull was rinsed with sterile saline and 1 cm 3 of ultrasound coupling gel was placed on the window. The linear ultrasound probe (15 MHz central frequency, 160 elements; Vermon, Tours, France) was positioned directly above the cranial window surface. We used the central 128 elements of the transducer, which was connected to an ultrafast ultrasound scanner (Aixplorer, SuperSonic Imagine, Aix-en-Provence, France). The software-based architecture of the scanner enabled Matlab (MathWorks, Natick, Massachusetts, USA) programming of custom transmit/receive ultrasound sequences. At the beginning of each session, an anteroposterior Doppler scan was performed to visualize and localize the shape of the blood vessels and the corresponding anteroposterior coordinates. Acquisitions were performed at three anteroposterior coordinates: Bregma þ 0.84, À 0.6 and À 2.16 mm. Note that ' þ ' (positive) signs mean 'in front of' (or 'rostral to') the Bregma, while ' À ' (negative) signs mean 'behind' (or 'caudal to') the Bregma, corresponding to the rat brain atlas of Paxinos 22 .
Ultrasound sequences. The concept of ultrafast Doppler relies on compounded plane-wave transmissions 70 . The brain was insonified with a succession of ultrasound plane waves. The backscattered echoes were recorded and beamformed to produce an echographic image for each transmission. The frame rate of ultrafast ultrasound can reach more than 10 kHz. However, a 500-Hz frame rate allows correct sampling of the ultrasound signals backscattered by the red blood cells without aliasing, in the rat brain, as demonstrated inref. 31. To increase the SNR of each echographic image taken at 500 Hz, we compounded the echographic images by transmitting several tilted plane waves and added their backscattered echoes. This compounded sequence resulted in enhanced echographic images, increasing the sensitivity of the Doppler measurement. In this study, the ultrasound sequence consisted of transmitting five different tilted plane waves ( À 4°, À 2°, 0°, þ 2°and þ 4°tilted angle) with a 2,500-Hz pulse repetition frequency (PRF). The backscattered echoes were added to produce enhanced echographic images at a 500-Hz frame rate.
The current technological limit of the ultrafast imaging system is the high rate of recorded raw data (2 GB s À 1 ), which cannot be stored in the computer memory and has to be treated in real time. This treatment consists of a beamforming step (image formation), which is time consuming. Therefore, we decided to perform two sequences: Sequence number 1 consisted of insonifying the rat brain for 400 ms (2 cardiac cycles, 5 angles and a 2,500-Hz PRF), every 2 s, with sufficient spare time between acquisitions for in-depth beamforming of all the raw data. This allowed us to study the FC of both superficial and deep brain structures for acquisition periods lasting several minutes.
Sequence number 2 extended the temporal coverage of fUS acquisition by trading off spatial coverage. It consisted of in-depth insonification of the whole brain continuously for minutes at a 2,500-Hz PRF, but only beamforming in real time the dorsal part of the rat brain (the neocortex) because limited computer memory (6 hyperthreaded core processing unit and 24 GB random access memory). Indeed, limiting the region of interest to the neocortex resulted in 3 Â less raw data and led to viable real-time processing with our current research platform.
Power Doppler data treatment. The backscattered signals from the rat brain were composed of tissue and blood signals. As the blood moves faster than tissue, its signal frequency is higher and can be extracted by time filtering the data with a high-pass filter called clutter filtering (fourth order Butterworth filter with a 75-Hz frequency cutoff). For sequence number 1, as the acquisition was performed by blocks of 400 ms of insonification, the time fluctuations of the blood signal of each block were incoherently averaged to obtain a power Doppler image proportional to the cerebral blood volume. The sequence number 2 data treatment consisted of simply removing the tissue signal with the same high-pass filter.
Building activation maps. Maps of 'activated' pixels were built using the Pearson correlation coefficient r between the local power Doppler temporal signal computed from each spatial pixel of the fUS acquisition and the temporal binary pattern induced by the electrical stimulations of the sciatic nerve. Activations were considered significant for a correlation r42s, where s is the spatial s.d. of the correlation map. The time course for a given region was calculated by averaging the power Doppler signal over time for all pixels in the activated region (r42s).The amplitude of the power Doppler was represented as the percentage of change relative to the baseline in the activated region±s.d.
Building FC maps. We employed the seed-voxel approach, which used the time course of the haemodynamic signal extracted from a ROI (the seed region), and determined the temporal correlation between its signal and the time course from all other brain voxels 10 . This created a whole-brain voxel-wise FC map of covariance with the seed region. This is a hypothesis-driven method, because it gives direct answers to specific hypotheses about the FC of the seed region. Such 'seed-based' FC maps correspond to maps of the Pearson's product-moment correlation coefficient r between the time signal a(t), which is the mean signal of the seed region and the power Doppler signal s d (t) for each pixel, defined as: where Bdesignates the mean value of the variable over time. To superimpose 'seed-based' FC maps on the Doppler maps, we displayed only the pixels of the 'seed-based' FC maps in which coefficient r42s, where s is the spatial s.d. of the correlation map. The reproducibility of the 'seed-based' FC maps between the different animals was determined without thresholding.
To increase the SNR, the time signal of each spatial pixel was filtered with a band-pass filter with cutoff frequencies of 0.05 and 0.2 Hz. The filtering of frequencies below 0.05 Hz was chosen to correct for a drift in the baseline, which could bias the correlation value. Filtering of frequencies over 0.2 Hz was chosen to remove noise, because we assumed that the frequencies of brain FC were below 0.2 Hz.
SVD of the fUS data. Alternatively, data-driven analysis, such as SVD or ICA, can also be used to analyse whole-brain FC patterns without the need for a priori information. Here we built FC maps based on SVD analysis. The FC maps corresponded to the spatial singular vector of the spatiotemporal matrix containing the power Doppler signal time profiles of all brain pixels (see Methods). SVD processing decomposed the fUS data into separate space and time variables, without the need for computational optimization or a priori information. We searched for the GSMs of the SVD that represented the FC on N acquisitions of different animals for each coronal plane. One fUS acquisition could be interpreted as a three-dimensional matrix: a two-dimensional spatial map of vascularization was acquired at multiple time points. We introduced fUS(x, z, t) as the spatiotemporal matrix of one fUS acquisition, where x and z are the spatial variables that describe the lateral distance and depth of one power Doppler image, respectively, and t indicates the time at which the power Doppler image was acquired.
The SVD may be thought of as decomposing a matrix into a weighted, ordered sum of the separate matrices. Consequently, the acquisition matrix fUS(x, z, t) was written after SVD as: fUS x; z; t ð Þ¼ X r l r u r ðx; zÞv r ðtÞ ð 2Þ NATURE COMMUNICATIONS | DOI: 10.1038/ncomms6023 ARTICLE Using SVD, fUS(x, z, t) was decomposed into different orthogonal, spatiotemporal representations, each expressed by u r (x, z) and v r (t): one spatial image, modulated by a temporal signal, v r (t). In other words, all pixels of the image u r (x, z) behaved with the same fluctuating time signal v r (t). The main advantage of this orthogonal representation of the matrix fUS(x, z, t) was that each spatiotemporal representation was weighted by l r ; this parameter allows the ranking of these spatiotemporal representations according to their energy weight in the matrix fUS(x, z, t). A high value of l r means that the associated spatiotemporal representation u r (x, z) and v r (t) are highly representative of the spatiotemporal fluctuations of the matrix, that is, they describe the FC of the brain, whereas a low value of l r means that the associated spatiotemporal representation u r (x, z) and v r (t) may represent only noise. Using the SVD on fUS data permits extraction of the main characteristics of the FC by choosing the spatiotemporal representation with the highest associated l r .
The main remaining problem with SVD is that a clear threshold on the physically relevant values of l r is difficult to define with an unbiased approach. For this reason, the method proposed in this article aimed to extract functional brain connectivity patterns directly from the entire population of N rats. Using the first SVD of each individual fUS acquisition as a noise filter, we suppressed the noise space associated with the lowest singular values from the spatiotemporal representations of the fUS data of each individual animal. Then, for each coronal plane we integrated the fUS acquisition data of N rats and applied a second SVD to this global matrix, representing the combined noise-filtered data of N rats to find the GSM. The last processing step consisted of establishing a clear threshold to take advantage of the several N fUS acquisitions. The specific steps of the algorithm are detailed in the Methods and are illustrated in Fig. 8.
Reshaping the fUS acquisitions. For each coronal plane, one fUS acquisition was performed on N different rats. One fUS acquisition output is a three-dimensional matrix of power Doppler images of size N X (samples on the lateral direction) and N z (samples in the depth direction), and acquired at N T different times. Figure 8a shows the resulting three-dimensional fUS matrix. We start by changing every three-dimensional acquisition (depth, lateral direction and time) into a twodimensional matrix by reshaping the space dimensions (depth and lateral direction) on one side (N S ¼ N X Â N Z samples) and the time dimension on the other side (N T samples).
For each fUS acquisition (rat number i), we obtain a two-dimensional matrix A i , iA{1,2,y,N} for which each column vector represents a Doppler image. The reshaped matrix for the first rat A 1 is displayed in Fig. 8b.
Centring and normalizing the data. The SVD is based on an energy criterion. Therefore, to take in account all of the scales of spatiotemporal fluctuations, we centre the time profile (mean value of the temporal variation centred to zero) and normalize each time profile of the fUS acquisition. Therefore, for each matrix A i , iA{1,2,y,N}, each row vector is centred and normalized. For the sake of clarity, we will still call this matrix A i . SVD of A i An SVD for each A i is performed and this matrix can be written as: where T stands for the matrix transposition, U i the matrix whose column vectors are the spatial singular vectors of the matrix A i , V i the matrix whose column vectors are the temporal singular vectors of the matrix A i and S i a diagonal matrix of coefficients containing the singular value of the matrix A i . (Note that the number of non-zero diagonal coefficients is the exact rank of the matrix A i ). In fact, the r th column vector of U i and V i can be identified as an image and a time signal, respectively, defining a spatiotemporal representation of the matrix A i with the r th diagonal coefficient of S i , the associated weight.
The results for the first acquisition U 1 , S 1 and V 1 of the SVD for the first rat (matrix A 1 ) are shown in Fig. 9c.
Noise removal. The number of singular vectors where a non-zero singular value is found is equal to the smallest dimension of the matrix A i . In our case, as the number of spatial samples (N S B10,000) exceeds the number of temporal samples (N T ¼ 300), we obtain N T singular vectors. Most of these vectors contain only noise and are irrelevant for determining the GSM of the brain FC. Therefore, we keep only the N R singular vectors that correspond to the highest singular values, and we increase the SNR of these matrices by removing the last singular vectors. The enhanced matrix computed from A i can be expressed as: where S SNR i is a diagonal matrix of the same size as S i containing only the first N R singular value of the matrix A i .
Extracting spatial modes of FC from fUS acquisitions. If we look at the spatiotemporal characteristics of the FC, we notice that the temporal behaviour of the FC signal is random in different animals; only the spatial characteristics will be conserved. Therefore, to unify the N acquisitions from each fUS acquisition, we only keep the highest spatial singular vectors with their weight (singular value), making only the first N R vectors of the matrices U i and S SNR i relevant. From this set of relevant singular vectors, we can compute the matrix A final i , whose size is N S Â N R , for each fUS acquisition. These matrix column vectors are the N R first column vectors of U i multiplied by their singular value contained in the diagonal elements of S SNR i . Unifying the information of the N fUS acquisitions. To unify the information over all of the N fUS acquisitions, we create a matrix A tot , which is the concatenation of the A final i matrices computed from the fUS acquisitions on the i ¼ 1 to N rodents. A tot can be expressed as: A tot ¼ A final 1 ; A final 2 ; :::; A final N Â Ã . Figure 9d shows this matrix A tot .
Computing the GSM of the FC. The GSM of the FC is found using SVD of A tot . To compute such a process, the line vectors of this matrix must be first centred and normalized as in step 2. For the sake of clarity, we will still call this matrix A tot . We can decompose the matrix A tot as: The GSM can be computed within a matrix M. These matrix column vectors are the N R first column vectors of U tot multiplied by their singular value contained in the diagonal elements of S tot . We only keep the first N R column vector of U tot , because we can consider that there are less relevant global vectors than N R . The GSM can be computed by reshaping the two-dimensional matrix M into a three-dimensional matrix. We performed the opposite process of step 1, except that the third dimension of the reshaped matrix will be the rank of the GSM. The remaining problem is that M contains more column vectors than the number of GSMs representing the FC in the total population of N animals. The column vector of M can contain specific FC in one animal but only biological noise in other animals, and the last column vectors of M can contain only acquisition noise. Therefore, we must find an efficient way to threshold which vector column of M is representative for the FC in the total population of N animals and can be called the GSM (global meaning relevant for all animals). To measure whether one vector column of M is a GSM, we studied its strength in the FC of each animal by projecting it onto each matrix A final i . If this vector is strongly present in each A final i matrix, it was considered a GSM.
Projecting the column vector of M on each fUS acquisition. To explain why we project the column vector of M, we take the example of the k th column vector, which can be expressed as a space vector m k 2 R Ns , as this vector has N S coefficients. We shall also express the matrix describing the spatial FC of the i th rat A final i as vectors: the column r of this matrix can be expressed as a i r 2 R Ns with rA{1,2,...,N R } , these N R vectors are orthogonal due to the SVD properties.
The main question is how can we measure the strength of the vector m k among the vectors a i r , rA{1,2,...,N R }? We can project the vector m k using the scalar product and introduce the projection of m k on the vectors a i r , rA{1,2,...,N R }: m i k . m i k ¼ P NR r¼1 a i;k r a i r with a i;k r ¼ a i r j m k is the scalar product between the vector a i r and the vector m k .
Next, we can measure the collinearity of m k and its projection m i k using a normalized scalar product: If a space vector m k is perfectly represented in the FC of the i th rat, that is, one a i r is collinear to m k so m i k ¼ a i;k r a i r $ m k , the associated c i k will be equal to 1. If m k is noise, the value of the associated c i k will be lower. We will still have a residual colinearity, because at step seven the m k are computed from an SVD of the matrix A tot containing the a i r , rA{1,2,...,N R }. As we aim at finding the GSM on N rats, we will study the average of c i k on N animals: c k can be expressed as: In Fig. 9e, blue dots show the c k for k ¼ 1 to N R for the coronal level at Bregma À 0.6 mm (in this study we chose N R ¼ 150). We notice that the first six coefficients c k have a high value (40.6), which means their associated column vector m k is strongly present in the matrix A final i . The c k for k varying from 7 to 35 have a low value that fluctuates. This means that their associated m k represents some functional characteristic of a few rats over the total N rats. As they are not present in all rats, we consider that they are highlighting the biological noise and we do not keep them as relevant singular vectors. Therefore, the last c k (k435) are due to the acquisition noise.
Thresholding the c k distribution. To accurately threshold which lowest value of c k is representative of the global FC over the total N rats and is not due to any kind of noise, we can apply the same process (the height first steps) to a set of N threedimensional random matrices (N X Â N X Â N Z ) containing only Gaussian random noise. The obtained mean coefficients (c noise k ) are displayed by red dots in Fig. 8e.
We clearly see that it matches the last c k (k435) of the fUS experimental acquisitions, because these last are due to the acquisition noise, which can be considered as random. We can select the GSM, which has a c k 4c noise k þ 2Dc noise k , where Dc noise k is the s.d. of c noise k over N rats: This curve is displayed in Fig. 8e with the red line. At the end, we select the first six column vectors of M to be the GSM of the FC over N rats for the coronal slice at Bregma À 0.6 mm. These first six GSM are displayed in Fig. 9f. The same process is applied to the other two coronal slices.