Identifying the bioimaging features of Alzheimer's disease based on pupillary light response-driven brain-wide fMRI in awake mice Xiaochen Liu , David Hike , Sangcheon Choi, Weitao Man, Chongzhao Ran, Xiaoqing Alice Zhou , Yuanyuan Jiang & Xin Yu Pupil dynamics has emerged as a critical non-invasive indicator of brain state changes. In particular, pupillary-light-responses (PLR) in Alzheimer's disease (AD) patients show potential as biomarkers for brain degeneration. To investigate AD-specific PLR and its underlying neuromodulatory sources, we combine high-resolution awake mouse fMRI with real-time pupillometry to map brain-wide event-related correlation patterns based on illuminationdriven pupil constriction (P c ) and post-illumination pupil dilation recovery (amplitude, P d , and time, T). The P c -driven differential analysis reveals altered visual signal processing and reduced thalamocortical activation in AD mice in comparison with wild-type (WT) control mice. In contrast, the post-illumination pupil dilation recovery-based fMRI highlights multiple brain areas associated with AD brain degeneration, including the cingulate cortex, hippocampus, septal area of the basal forebrain, medial raphe nucleus, and pontine reticular nuclei (PRN). Additionally, the brain-wide functional connectivity analysis highlights the most significant changes in PRN of AD mice, which serves as the major subcortical relay nuclei underlying oculomotor function. This work integrates non-invasive pupil-fMRI measurements in preclinical models to identify pupillary biomarkers based on brain-wide functional changes, including neuromodulatory dysfunction coupled with AD brain degeneration.

Pupil diameter changes reflect cognitive processing [1][2][3][4] , presenting a unique, non-invasive index of brain state fluctuation in both normal and degenerative conditions, e.g., Alzheimer's Disease (AD) [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23] . Although several studies associate pupillary responses with the risk for AD, it remains highly speculative given the lack of understanding of its mechanistic linkage [24][25][26][27][28] . Multiple neuronal sources drive the spontaneous pupil dynamics, as well as pupillary light responses (PLR), through autonomic neuromodulation to control pupil dilation and constriction [9][10][11][29][30][31][32][33][34][35][36] . Besides evoked visual pathways, widespread neuromodulatory pathways, including those involving acetylcholine [37][38][39] , norepinephrine 39,40 , and serotonin 41 , are recruited to regulate PLR. In addition, several task-related studies have reported that PLR is mediated by higher-order brain function with either direct cortical microstimulation 42 or loss-of-function by brain lesion 43 . Given that ADrelated cognitive decline is linked to impaired neuromodulatory circuits, particularly cholinergic deficits [44][45][46][47] , an intriguing hypothesis proposes that altered PLR could serve as an indirect indicator of brainwide dysfunction resulting from neuromodulatory impairment in AD.

Combining fMRI with real-time pupillometry provides a noninvasive mapping approach that bridges altered pupillary responses with corresponding brain-wide functional changes in AD brains. Existing studies have primarily focused on analyzing spontaneous pupil diameter changes with resting-state (rs) fMRI signals. Human fMRI studies with pupillometry have revealed that pupil dynamic changes are associated with the activation of subcortical brain nuclei along the ascending arousal network, as well as the salience network including the cingulate cortex and insula 32,35,[48][49][50][51][52][53][54][55] . Elman et al. have linked the locus coeruleus (LC)-related functional network with the task-related pupil responses, indicating potential LC malfunctionbased pupil dynamic biomarkers for AD risk 22 . Whereas, the noradrenergic regulation of pupil-fMRI interactions, e.g., the noradrenergic cell group 5 (A5) nuclei at the brainstem, was detectable only at a certain brain state of rat brains 56 . By measuring the axonal activity of neuromodulatory neurons, it has been shown that pupillary responses to locomotion are more closely coupled with sustained cholinergic activity, but less with noradrenergic phasic activity 57 . This observation is further confirmed in the AD transgenic mice, which exhibit pupil dilation coupled with locomotion but less associated with the diminished astrocytic activity controlled by noradrenergic inputs 58 . Besides cholinergic and noradrenergic systems, other neuromodulatory pathways are also reported to regulate pupil dynamics at resting state 41 . A recent work demonstrates the highly varied spectral distribution of spontaneous pupil dynamics according to different activation patterns of neuromodulatory systems 57 . These studies have illustrated the challenges of dissecting the neuromodulatory contribution to spontaneous pupil dynamics based on rs-fMRI correlation analysis, which would be insufficient when bridging the AD-specific pupil dynamics changes with brain-wide functional alteration due to neuromodulatory dysfunction in AD brains.

Instead of relying on spontaneous pupil dynamics, light-triggered PLR responses can be controlled with fixed interstimulus durations, allowing for phase-locked analysis of fMRI and pupillary changes. In this study, an approach combining awake mouse fMRI with simultaneous real-time pupillometry was developed to elucidate altered brain function and potential neuromodulatory dysfunction underlying AD-specific PLR. Using an implantable radiofrequency (RF) surface coil, which also served as a headpost for head-fixation during scanning, robust high-resolution brain-wide BOLD functional maps (100 × 100 × 200 µm resolution) were acquired from awake mice undergoing visual stimulation. Epoch-specific PLR responses were measured during awake mouse fMRI, including the illuminationinduced pupil constriction amplitude (P c ) and post-illumination pupil dilation recovery time (T) and amplitude (P d ). Using these time-varied PLR dynamic features, the P c -based fMRI differential maps between AD and wild-type (WT) mice were presented to identify altered brain activation in the visual pathways. In contrast, post-illumination pupil dilation features (T, P d )-based fMRI differential maps highlighted several cortical (e.g., cingulate and retrosplenial cortex (RSP)) and subcortical regions (e.g., septal areas in the basal forebrain, hippocampus, and pontine reticular nucleus (PRN)) involved in the altered neuromodulation underlying the AD-specific PLR. Meanwhile, the region of interest (ROI)-based correlation analysis also identified altered functional connectivity among these PLR-related neuromodulatory nuclei, highlighting the most significant connectivity changes originating from the PRN between AD and WT mice. These results establish a pupil-fMRI mapping scheme for identifying AD-specific PLR dynamic features, revealing altered correlation patterns in functional nuclei potentially affected during brain degeneration.

The real-time pupillometry was simultaneously performed inside the 14 Tesla scanner with high-resolution echo-planar-imaging (EPI) scans of awake mice (Fig. 1A). The awake mouse real-time pupillometry setup includes i. MR-compatible miniaturized complementary metal-oxidesemiconductor (CMOS) sensor with a focus adjustable lens to maintain at least a 2 cm distance from the mouse to avoid B0 field distortion, ii. A plane mirror to reflect the pupil toward the camera mounted on the animal holder, iii. The optical fibers to deliver the switchable blue/ green light pulses (530 nm and 490 nm, flashed at 5 Hz and 5.1 Hz, respectively) for visual stimulation and the red light (660 nm) to illuminate the eye for pupillometry, as mice are 1000 times less sensitive to red light with wavelengths over 660 nm 59 . Mice were implanted with an RF surface coil that served as the headpost for head-fixation while inside the animal holder during scanning (the detailed setup was described in our previous work 60 ). Figure 1B showed a representative trace of the PLR from an awake mouse, measured using DeepLabCut 61 . To improve the accuracy and efficiency of PLR measurements, we also implemented a grayscale-based pupil size measurement to verify different PLR dynamic features across different pupillometry recording trials (Supplementary Fig. 1). The mean PLRs of two groups of mice (AD vs. WT) were plotted in Fig. 1C, showing the dynamic PLR features during constriction (P c ) and dilation recovery phases (P d , T). The recovery time (T) was quantified based on exponential fitting of the pupillary dilation curve from each stimulation on/off epoch. Figure 1D showed the quantitative comparison of mean PLR features, presenting a significantly shorter T in AD mice compared to WT mice, but no significant difference in P c and P d between the two groups. These results provided robust PLR dynamic measurements of head-fixed awake mice during fMRI scanning. However, the altered PLR dynamic features detected in AD mice still need to be defined as potential biomarkers of AD. Instead, combining PLR measurements with simultaneous fMRI to identify the brain-wide functional changes underlying AD-specific PLR features may provide a promising approach to identifying bioimaging markers of AD.

High-resolution fMRI was applied to acquire EPI images of awake mice with 100x100x200µm resolution using a horizontal 14 T MRI scanner (Supplementary Fig. 2). During the awake mouse fMRI scanning, motion artifacts caused by the voluntary movement of the mice led to distorted images (Supplementary Fig. 2E), which can be removed from the 3D time series for regression analysis (Supplementary Fig. 2F, Movie 1). Figure 2 showed the visual stimulation-evoked brain-wide activation in the visual cortex (VC), superior colliculus (SC), lateral geniculate nucleus (LGN), retrosplenial cortex (RSP), and anterior cingulate cortex (ACA) in both WT and AD mice. Robust positive BOLD signals were detected in both groups from the time course data of specific ROIs segmented based on the Allen mouse brain atlas 62 (Fig. 3), demonstrating the reliable visual stimulation-evoked fMRI signals in awake WT and AD mice.

To investigate the difference of functional maps between the two groups, we first integrated the P c -based amplitude modulated (AM) hemodynamic response function (HRF) as the regressor across different trials. Since the varied pupil constriction at each stimulation on/ off epoch altered retinal illumination, specific visual pathway activation would be better quantified by integrating the P c -based AM regression scheme. Quantitative group analysis was performed to produce voxel-wise differential maps between WT and AD mice. Figure 4 showed significantly higher BOLD responses in WT mice compared to AD mice in the primary and secondary visual cortices and LGN, but little difference was detected in the SC and cingulate cortex. In particular, well-separated higher visual cortical areas, including the adjacent RSP, were highlighted in the differential maps (Fig. 4B). These results indicate that the higher-order thalamocortical processing of visual signals is impaired in AD mice based on the P c -fMRI AM regression analysis.

Brain-wide functional changes associated with post-illumination pupil dilation were investigated using the pupil-fMRI mapping scheme. The pupil dilation features P d and T were quantified for each stimulation on/off epoch with awake mouse fMRI, enabling a PLR-based fMRI regression analysis. Figure 5 showed the P d -based fMRI correlation patterns in both WT and AD mice. In both groups, a negative correlation was observed in subcortical midbrain regions. In WT mice, this negative correlation was located in the serotonergic raphe nuclei, while in AD mice it was predominantly detected in sparse pontine areas adjacent to the LC, laterodorsal tegmental area (LDT), and the mesencephalic pedunculopontine tegmental nuclei. Also, P d -based negative correlation with the cingulate cortex, hippocampus, and lateral septal (LS) area were detected in AD mice. Group differential analysis revealed the most salient differences in the MnR, the cingulate/retrosplenial cortex, the hippocampus (dentate gyrus (DG) and CA1), and the LS of the basal forebrain regions.

Figure 6 showed the T-based fMRI correlation patterns, highlighting the negative correlations in the PRN, the cingulate cortex, and the hippocampus in WT mice. Furthermore, negative correlations were observed in the midbrain reticular formation and central thalamic nuclei in WT mice. In contrast, AD mice showed positive correlations across the hippocampus and cingulate cortex, and the hypothalamus. The statistical differential maps between WT and AD mice highlighted the significantly different brain regions located at the PRN, the hippocampus, the cingulate/retrosplenial cortex, and central thalamic nuclei, including the paraventricular (PVP) and dorsal thalamus.

To further identify the neuromodulatory nuclei responsible for the altered pupil-fMRI coupling in AD brains, we also performed ROIbased correlation network analysis given the reiterated visual stimulation paradigms. Figure 7A illustrated the spatial correlation network of functional connectivity among the highlighted nuclei from PLRbased differential maps between WT and AD mice. By comparing the ROI-based correlation between WT and AD mice, significantly altered connectivity among paired functional nuclei was highlighted in the spatial correlation network maps, showing increased connectivity between hippocampus-PRN, LS-PRN, and decreased connectivity between LS-VC/SC, ACA-VC/SC in AD mice. Figure 7B showed the corresponding correlation matrices, summarizing the paired functional nuclei with significantly different connectivity between the two groups. It should be noted that we also conducted permutation controls using randomly selected brain-wide regions of interest (ROIs) for differential analysis. This approach revealed a significantly reduced number of different connections compared to those identified in the nuclei selected based on the PLR-based fMRI analysis (Fig. 7C). Meanwhile, brain-wide correlation analysis of ROIs, independent of the PLR-based fMRI mapping results, was performed to identify the ROIs with the most significantly altered connectivity in AD mice. The top 20 paired ROIs with significantly altered connectivity between WT and AD mice highlighted that over 40% of the connectivity changes originated from PRNc/PRNr in AD mice (Fig. 7D), as identified by the PLR-based fMRI differential analysis. These results further confirm that the PLRbased fMRI differential analysis identified a critical set of neuromodulatory nuclei most responsible for AD-specific PLR.

The present study combined awake mouse fMRI with real-time pupillometry to identify brain-wide functional changes associated with ADspecific PLR dynamics. This awake mouse pupil-fMRI mapping scheme provides a proof-of-concept method for identifying non-invasive PLR linked to altered brain function, particularly within neuromodulatory pathways, as potential biomarkers of AD brain degeneration.

PLR in AD patients has been reported to show altered pupil constriction phases [24][25][26] . These studies have suggested the potential of pupil dynamics as biomarkers of AD degeneration. However, a few studies have failed to reproduce altered PLR features when comparing AD patients with control groups 27,28 . Numerous efforts have been made to link pupil dynamics with various behavioral variables such as the subject's arousal fluctuation [5][6][7] , designated tasks 5,[8][9][10][11][12] , or specific cognitive processes [13][14][15][16][17][18][19][20][21] . The richness of behavioral correlates of pupil dynamics, especially in AD patients with cognitive impairments, could lead to confounding PLR results given the complex controlling mechanisms. In the 5xFAD mouse model, altered PLR features, such as the significantly shortened post-illumination dilation recovery time (T) (Fig. 1), could serve as alternative behavioral readouts reflecting brain functional changes associated with cognitive impairments, which have been extensively studied using behavioral tests in this transgenic mouse line [63][64][65][66][67] . Besides the well-documented behavioral impairments, detailed studies have also identified various neuroanatomical abnormalities in 5xFAD mice [68][69][70][71] . Our work addresses the gap in bridging PLR features with brain-wide functional changes using fMRI in 5xFAD mice, which could directly contribute to the reported cognitive impairments in this mouse model.

One intriguing technological development of this work is the high-resolution awake mouse fMRI at 14 T. Hike et al. developed a novel implantable RF coil, which also serves as the headpost for headfixation during awake mouse fMRI scanning 60 . In contrast to the inductively coupled wireless RF coil 72,73 , this new design has enabled ultra-high-resolution EPI-based fMRI of awake mice with 100 × 100 × 200 um resolution, allowing the characterization of activated functional nuclei in mice with only a few hundred-micron size. Both AD and WT mice have gone through the acclimatization procedure with fMRI as described previously 60,74 , enabling the reliable detection of visual stimulation-evoked fMRI maps (Figs. 2 and 3). It should be noted that the functional maps showed brain-wide activation not only in the visual pathway but also included the thalamic areas and ACA in both groups. This demonstrated extended high-level cortical function involved in the cognitive processing of awake mice. This advanced awake mouse fMRI allowed the brain-wide functional characterization of AD-specific PLR features.

We have performed PLR-based fMRI correlation analysis based on the stimulation on/off epoch across different trials. In contrast to the direct comparison of evoked functional maps between WT and AD groups, the P c -based AM regression analysis enabled the characterization of evoked brain function based on illumination-driven retinal activity. The BOLD signal difference detected along the central visual pathways indicated significantly reduced visual cortical activity from primary and associated cortices in AD mice, as well as in the LGN, but no signal difference was detected in the central thalamic areas and cingulate cortex (Fig. 3A). This altered visual pathway activity in AD mice could be directly linked with impaired retina function due to vascular and neuronal degenerative effects [75][76][77][78] . The differential maps highlighted V1 and associated visual areas with LGN, indicating impaired high-order visual signal processing. Previous studies have reported altered laminar distribution of neurofibrillary tangles (NFTs) in the visual cortex 79 and decreased long synaptic projection from area V2 to V1 in AD patients 80 . The P c -based fMRI provides a bottom-up mapping scheme to characterize the malfunction of visual signal processing in AD brains.

In contrast to the P c -based fMRI correlations, the post-illumination pupil dilation reflected the top-down neuromodulatory processes. The P d =T-based fMRI revealed altered coupling with multiple functional nuclei in AD mice. First, the posterior cingulate cortex and hippocampus of AD mice showed altered coupling features with P d =T compared to WT mice. Since both the posterior cingulate cortex and hippocampus are impaired in AD due to the neurodegeneration [81][82][83][84][85] , it is plausible that the reduced neuronal activation in these areas was directly involved in the altered pupil dilation features, especially the shortened dilation recovery time (T) that presented the most significant changes in PLR measurements (Fig. 1C). This observation is consistent with reports, showing that hippocampal damage in AD patients is related with impaired oculomotor behavior and networks [86][87][88] as well as altered pupil response [89][90][91][92] . From the ROI-based correlation analysis, we observed fewer or weaker connections related to the hippocampus (CA1) and cingulate cortex in AD mice compared to WT mice. This finding is consistent with other functional connectivity studies in AD mice [93][94][95] .

Another intriguing observation is the identification of subcortical brain regions underlying the altered pupil dilation recovery in AD mice. First, the PLR (T)-based fMRI coupling in PRN was reduced in AD mice compared to WT mice. The PRN is the major relay nuclei in the brainstem controlling oculomotor function [96][97][98] . By performing the brain-wide ROI network connectivity analysis based on the reiterated visual stimulation (Fig. 7D), significant changes in PRN-based network connectivity in AD mice further confirmed the PLR (T)-based PRN coupling as a crucial estimate of the reliable pupil-fMRI biomarker for AD brains. For the P d -based fMRI correlation, negative correlations at median raphe nuclei (MnR) were observed in WT mice but were less pronounced in AD brains, indicating impaired serotonergic regulation of pupil dilation. The pupil dilation and constriction can be controlled by serotonergic agonists and antagonists, respectively [99][100][101] , and significant neuronal reduction is also reported in the MnR of AD patients' brains 102 . In contrast, the LS showed more negative correlations in AD than WT mice in the P d -based fMRI maps. The septal area of the basal forebrain is composed of glutamatergic, GABAergic, and cholinergic neurons and is reciprocally connected to multiple subcortical regions 103 , including the medial septal area (MS) which sends cholinergic projections throughout the cortex 104 . Cholinergic hypofunction is a hallmark of AD 24,105 , characterized by impaired cholinergic projections from septal areas to widespread brain regions, as observed in various AD transgenic models [106][107][108][109][110][111] , including 5xFAD mice 112 . Since cholinergic function directly influences pupil dilation 57 , the increased negative correlation in LS suggests that LS-mediated septal cholinergic function in AD brains contributes more significantly to regulate postillumination pupillary dilation compared to healthy brains. Our work points to a future research direction on investigating the causal relationship between cholinergic deficits and altered PLR dynamics in AD. In particular, the multi-modal fMRI platform, which combines simultaneous fMRI and fiber photometry to record acetylcholine signals, can be integrated with chemogenetic tools to regulate the cholinergic neuronal activity in the septal areas of awake AD mice.

A counterintuitive observation in this work is the primary negative correlation patterns from multiple nuclei for the P d -based fMRI maps. The pupil-based fMRI negative correlation has been previously reported in anesthetized rodent brains 56 . And brain-wide negative BOLD signal was also coupled with the vigilant brain state changes based on eye open/close arousal indices 113 . The P d -based negative BOLD fMRI correlation could be mediated by the neurovascular coupling event underlying arousal state changes in awake mice involving multiple neuromodulatory pathways. By incorporating the restingstate spatiotemporal pupil-fMRI correlation patterns, principal component analysis (PCA) has been implemented to highlight different neuromodulatory pathways based on clustered pupil dynamic spectra 114 . To further elucidate the neuronal basis of the negative coupled pupil-fMRI relationship, cellular-specific fiber photometry Ca 2+ recordings 56 need to be implemented with awake mouse fMRI and real-time pupillometry. In addition, the brain state-dependent spontaneous pupil dynamics and locomotion-coupled pupil dynamics of awake mice could provide further insights into the coupled neuromodulatory dysfunction in AD mice.

In contrast to existing awake mouse fMRI studies 115 , we implemented a 5-week training scheme that has reported significantly reduced signs of stress-related behaviors, e.g., eye movements and spontaneous pupil dynamic changes 60 . Specifically, eye movements of head-fixed mice are much less frequent than in free-moving mice 116 MRI acoustic training 60 . In addition, stress-related pupillary dilation predominates in the spontaneous pupil dynamics during the early phase of training 74 but diminishes as acclimation progressed and leads to low-frequency fluctuations corresponding to brain state changes of normal behaving mice 60 . This work reported the highly comparable illuminance-induced pupil constriction between WT and AD mice, presenting similar acclimated states of both groups. And reduced corticosterone level of head-fixed mice has been reported with acclimation training in awake mouse fMRI studies 117 . However, corticosterone levels in head-fixed mice might still be higher than those in freemoving mice 117 , with variations depending on the duration of training and the length of the scanning session. Thus, the impact of stressrelated confounding factors on the fMRI-pupil relationship of awake AD mice requires further investigation in future studies.

In conclusion, our study demonstrates that the epoch-specific post-illumination PLR features, i.e., the pupil dilation recovery time (T) and amplitude (P d ), correlate with specific neuromodulatory nuclei that exhibit impaired functionality in AD brains. These nuclei also show dramatically altered network connectivity during visual stimulation in AD mice (Fig. 7), suggesting that AD-specific PLR features are linked to underlying changes in neuromodulatory function.

To further explore the relationship between PLR-related neuromodulatory dysfunction and cognitive impairments, future research will need to connect these altered functional maps with behavioral changes in AD mice across different stages of neurodegeneration. Eventually, this pupil-fMRI mapping could offer non-invasive measurements for identifying pupillary biomarkers of brain degeneration in AD patients.

A total of nine female 5xFAD transgenic AD mice and thirteen matched genetic background WT (C57BL/6 J) mice (6 females and 7 males) at 6-12 months old were used, similar to previous studies 60 . The age distribution of the mice during scanning is shown in Supplementary Fig. 3. For real-time pupillometry with fMRI, only nine (4 females and 5 males) of the thirteen WT mice were used because of the poor quality of the pupil recordings from inside the 14 T MRI scanner. For the AD mice, only females were used since female 5xFAD mice display more pronounced molecular changes, including cortical and hippocampal Aβ accumulation and other inflammatory gene induction, than male 5xFAD mice [64][65][66]118,119 . Also, female 5xFAD mice exhibit significant Alzheimer's-related behavioral and physiological changes by the age of 6 months [63][64][65][66][67] . Mice were group-housed (2-4/cage) under a 12-h light/ dark cycle with food and water ad libitum. All animal procedures were conducted in accordance with protocols approved by the Massachusetts General Hospital (MGH) Institutional Animal Care and Use Committee (IACUC), and animals were cared for according to the requirements of the National Research Council's Guide for the Care and Use of Laboratory Animals.

All mice underwent surgery under anesthesia (1-2% isoflurane with 1 L/min of medical air and 0.2 L/min additional O2 flow) to implant the RF coil to the skull 60 . Briefly, anesthetized mouse heads were stabilized in a stereotaxic stage, and the skull was exposed after the aseptic treatment. Before the coil was implanted, both 0.3% H 2 O 2 and PBS was applied to clean the skull. It is critical to wait until the skull is fully dried before the coil implantation. The coil was secured with cyanoacrylate glue and dental cement to prevent movement and ensure stability. Analgesic (Buprenorphine, 0.1-0.3 mg/kg) was administrated before the surgical procedure and was also injected periodically (once per 12-24 h) for 3 days after the surgery. Mice were allowed at least a week of recovery in their home cages for neck strengthening to resume normal head movements before training. Animals were acclimated to the MRI environment through a modified training program spanning five consecutive weeks. The training duration increased progressively increasing from 10 min to 60 min, transitioning from head-fixed mock MRI sessions to in-bore training with actual MRI scanning 60,74 . The stress levels of head-fixed mice during the training procedure were previously reported based on the monitored eye movements and pupil dynamics 60 . The stressrelated confounding issue is discussed in the Discussion section.

MRI data was acquired with the 14 T horizontal MRI scanner (Magnex Sci, UK) located at the Athinoula A. Martinos Center for Biomedical Imaging in Boston, MA. The magnet is equipped with a Bruker Avance Neo Console (Bruker-Biospin, Billerica, MA) and is operated using ParaVision 360 V. 3.3. A microimaging gradient system (Resonance Research, Inc., Billerica, MA) provides a peak gradient strength of 0 10 20 30 40 0 200 400 600 PRN_R PRN_L MnR_R MnR_L SC VC_R VC_L RSP CA1_R CA1_L ACA LS_R LS_L PRNc_R PRNc_L MnR_R MnR_L SC VC_R VC_L RSP CA1_R CA1_L ACA LS_R LS_L PRNc_R PRNc_L MnR_R MnR_L SC VC_R VC_L RSP CA1_R CA1_L ACA LS_R LS_L PRNc_R PRNc_L MnR_R MnR_L SC VC_R VC_L RSP CA1_R CA1_L ACA LS_R LS_L PRN_R PRN_L MnR_R MnR_L SC VC_R VC_L RSP CA1_R CA1_L ACA LS_R LS_L PRN_R PRN_L MnR_R MnR_L SC VC_R VC_L RSP CA1_R CA1_L ACA LS_R LS_L ACA LS SC RSP VC CA1 PRN MnR WT AD WT-AD 8% 2% 5% 8% 8% 5% 2% 2% 25% 20% 5% 5% 2% 2% DR IV MRN NLLh OT PAG PBG PPN PRNc PRNr PSV RR SCdg SCiw Number of significant different connections Number of repeats 1.2 T/m within a 60-mm diameter. The MRI data for anatomical registration was obtained with a multi-slice T1-weighted 2D gradient echo fast low angle shot (FLASH) sequence (TE: 3 ms, TR: 475 ms, flip angle: 30°, 4 averages for an approximate acquisition time of 4.5 min, 100 × 100 × 200 µm resolution). The awake mouse fMRI images were obtained with a multi-slice 2D EPI (TE: 6.2 ms, TR: 1 s with two segments, BW: 278 kHz; Matrix size: 144 × 96 × 36, in-plane resolution: 100 × 100 µm, thickness: 200 µm). The 6.2 ms TE was chosen based on the measured T2* values of the cortex (~7-8 ms) of awake mice at 14 T. A total of 205 repetitions are acquired per trial for a total of 6 min 50 s. For each mouse, three trials of visual stimulation were conducted per scanning session. Within each trial, a block-design paradigm consisting of 10 stimulation on/off epochs was applied. For each epoch, the visual stimulation lasted for 8 s, followed by 32 s of inter-stimulus intervals. The optical stimulation during the "8 s on" period included two LED lights flashing at specific frequencies: 490 nm at 5.1 Hz and 530 nm at 5 Hz, with light pulse duration at 20 ms.

The pupil recordings were obtained using a 3D-printed animal cradle with mounts for a camera, a mirror, and three fiber optic cables used as light sources as shown in Fig. 1A. The pupil dynamics were recorded through the mirror. The red light (660 nm) was used for video illumination since mice are 1000 times less sensitive to red light with wavelengths over 600 nm 59 , while the green and blue lights (530 nm and 490 nm) were used for visual stimulation.

To correlate the whole-brain fMRI signals with specific pupillary response features of normal and AD mice, the detection of pupil dynamics is first needed. Figure 1B shows an example of pupil dynamic recording. As shown in Fig. 1B, the edge of the mouse pupil was detected and tracked using DeepLabCut 61 . To further confirm the PLR results, an alternative pupil size detection method based on ocular and global grayscales was proposed. By comparing the average grayscale of the eye with the average grayscale of the whole image, the size of the pupil could be well characterized based on the signal intensity (Supplementary Fig. 1). Eventually, the normalized pupil dynamics data P can be calculated as:

Where g eye is the grayscale of the eye area and g all is the grayscale of the whole image.

To further map the correlation between pupil and fMRI, we extracted different features of the PLR in AD and WT mice with MATLAB R2023b (MathWorks, Natick, MA). The PLR time courses were divided into two components (constriction and dilation) for phasedependent analysis of PLR. The constriction component represents the immediate reaction of the pupil to illumination, primarily regulated by the parasympathetic nervous system. To quantify pupil constriction, we used the difference of the trial-specific z-scored pupil dynamic data before and after constriction, denoted as P c . For each epoch, the pre-constriction pupil size is calculated by averaging the pupil dynamic recording 1 s before the stimulation. The constricted pupil size is calculated by averaging the pupil recording below 75% of the epoch-specific minimum z-score. The dilation component, mainly controlled by the sympathetic nervous system, represents automatic pupil modulation and can provide insights into associated brain activity. Regarding the dilation component, we used an exponential function to accurately model the changes in pupil size.

Where PðtÞ is the pupil dilation defined from above 75% of the epochspecific minimum z-score after stimulation to 75% of the epochspecific maximum z-score. Through this modeling approach, we obtained two parameters: P d , indicating the difference in pupil size before and after dilation; T, characterizing the shape of the pupil dilation curve and the potential full dilation time.

To generate the visual-stimulation BOLD map, the fMRI images were processed using Analysis of Functional Neuroimages (AFNI_23.0.02) software 120,121 . A detailed awake mouse fMRI imaging processing procedure has been described previously 60 . First, the "to3d" function was used to transform the fMRI dataset into AFNI BRIK format. For group analysis, the images were registered to a mouse brain template. The mouse brain template acquired from the Australian Mouse Brain Mapping Consortium (AMBMC) 122 was resampled at 100 µm isotropic resolution. For each scanning session, the anatomical FLASH image was registered to the AMBMC template using the "3dAllineate" function, creating a transformation matrix (TM1). Meanwhile, averaged EPI data was registered to the anatomical image, creating a second transformation matrix (TM2). By integrating the two transformation matrices, EPI time courses from multiple trials per scanning session will be registered to the brain template for concatenation. After the registration, several pre-processing procedures were implemented before the general linear regression using commands in "afni_proc.py".

The time series for the scanning session was first concatenated using the "3dTcat" command. The "despike" command was used to truncate spikes in each voxel's time series before each EPI time point was registered to the atlas via a 6-degree transformation with the "volreg" command. The "mask" command was then applied to separate the brain from the background and other tissues. The "blur" command was used to smooth the transformed data by applying a Gaussian filter to average signal over nearby voxels before the "scale" command was applied to ensure compatibility and consistency across the analysis pipeline. In addition, the motion parameters were characterized and regressed out in the final general linear regression. The "censor_motion" parameter was set at 0.5 and the "censor_outlier" parameter was set at 0.02 based on the acquired awake mouse fMRI datasets (e.g., Supplementary Fig. 2F). The "3dDeconvolve" function was used to create the final statistical β coefficient maps with a "BLOCK" function.

To characterize PLR feature-driven fMRI maps, we applied the stimulation on/off epoch-specific PLR features to modulate the regression function. To improve the stability of regression, we concatenated the three trials from each scanning session and then obtained the PLR feature-driven fMRI maps. For group analysis, voxel-wise Student's t-tests were performed to identify the differential PLR-based fMRI correlation maps between AD and WT mice. The permutation-based false discovery rate (FDR) was calculated to verify the statistical results. We compared the observed statistic to the null distribution by shuffling the labels to estimate the permutation-based p-value.

To analyze the constriction component, the pupil constriction parameter P c was applied to the amplitude-modulated pupil-fMRI correlation analysis. The P c was used as the auxiliary behavioral information to modulate the amplitude of the hemodynamic response function. The AM1 function (Supplementary Fig. 4A) of AFNI was applied with a regressor calculated as the following:

Where P K k = 1 hðt À τ k Þ is the ideal HRF, P c ðkÞ is the value of the P c at time k. K is the time length of the fMRI data.

For the dilation component, pupil dilation parameters P d and T were used as the auxiliary behavioral information which was expected to be modulated by activation in specific brain areas. Those brain areas were expected to be associated with automatic pupil modulation and to exhibit variation between AD and WT groups. The AM2 function (Supplementary Fig. 4B) of AFNI was applied to incorporate the ideal HRF, as well as two regressors representing the epoch-specific modulation for pupil dilation parameters (P d and T), which were calculated as the following:

Where aðkÞ is the value of the pupil dilation parameter at time k. a is the average pupil dilation parameter.

For the statistical differential maps, significant voxels were presented based on the p-values ranging from less than 0.04 to 0.001, according to different comparison schemes. The minimal cluster size was set at 180-200 voxels according to different comparison schemes.

Two-sample t-tests on the three pupil-fMRI correlation maps (P c , P d , and T) were also performed to compare male and female WT mice, aiming to characterize the gender-related differences. The threshold for p-value and cluster size were set the same as in the WT-AD analysis. No significant difference was observed for PLR-based fMRI differential maps between male and female WT mice in the ROIs highlighted from the PLR-based fMRI differential maps between AD and WT mice (Supplementary Fig. 5).

For the correlation network analysis, we defined all the ROIs based on the Allen mouse brain atlas given the brain-wide delineation of the functional nuclei. This Allen mouse brain atlas was registered to the previously chosen AMBMC brain template. For each ROI, we averaged the time course from the preprocessed fMRI datasets based on the 21 sessions of 9 WT mice and 21 sessions of 9 AD mice. Pearson correlation analysis was performed to calculate the correlation coefficients among paired ROIs that were identified based on the PLRfeature-driven BOLD differential maps between WT and AD mice (a total of 13 ROIs were selected). A one-sample Student's t-test was performed to identify the significantly correlated paired ROIs in two groups of animals. A two-sample Student's t-test was performed to compare the correlation difference of the paired ROIs between the two groups. To verify the differential correlation networks of PLR-specific ROIs, permutation control tests were performed by randomly selecting the same number of ROIs (excluding those pupil feature-related ROIs) 5000 times. For each permutation test, the number of randomly paired ROIs with significantly different correlations between the two groups was calculated. Gaussian fitting was applied to verify the normal distribution for the permutation tests. Alternatively, Pearson correlation analysis was also performed across all the ROIs based on the Allen mouse brain atlas (a total of 1086 ROIs), and the two-sample Student's t-tests were performed to compare correlation coefficients of paired ROIs between the two groups. Among the paired ROIs showing significant correlation differences between the two groups (p < 0.01), the top 20 ROIs with the largest coefficient difference were chosen to plot the pie chart, which shows the individual ROIs' contribution to the altered correlation network between WT and AD mice. Supplementary Fig. 6 also showed the pie charts from the top 10 to the top 200 verified paired ROIs.