Authors: Kei Kanari, Taihei Tsutsumi
Categories: Article, Pupillometry, Generalized additive mixed models (GAMM), Sensory processing, Autonomic nervous system, Individual differences, Temporal profile, Neuroscience, Psychology
Source: Scientific Reports
Authors: Kei Kanari, Taihei Tsutsumi
Sensory processing traits describe how individuals respond to and regulate everyday sensory input, yet objective physiological markers of these traits remain difficult to establish. The pupillary light reflex has been proposed as a noninvasive window into sensory differences, but findings have been inconsistent, possibly because complex pupil dynamics are often reduced to single summary measures such as peak constriction. Here we analyzed the pupillary light reflex in 39 young adults. Participants completed the Adolescent/Adult Sensory Profile and questionnaires assessing psychological symptoms, and pupil responses were recorded to brief full-field light flashes in four colors. We quantified proportional pupil change over 0–3 s after stimulus onset using an approach that evaluates the entire response waveform rather than a small set of scalar indices. Trait associations appeared mainly as model-based differences in global fitted waveform shape rather than mean-level shifts. Effects were most pronounced for Sensation Avoiding. Waveform-based analyses may help characterize global temporal response profiles associated with sensory traits, although the present findings do not identify reliable divergence at specific time points or underlying mechanisms.
The online version contains supplementary material available at 10.1038/s41598-026-48654-5.
Sensory processing—the reception of sensory input, its central integration, and its translation into behavioral responses—varies widely across individuals and is thought to be continuously distributed in both clinical and nonclinical populations^1^. Such individual differences have been linked to everyday participation and quality of life^2,3^. Thus, beyond describing sensory processing traits as subjective tendencies, an important goal is to identify objective physiological signatures that can provide an empirical basis for these individual differences.
Dunn’s framework^4,5^ conceptualizes sensory processing styles as four patterns arising from the intersection of two neurological threshold (low–high) and self-regulation strategy (passive–active). The Adolescent/Adult Sensory Profile (AASP)^6^ operationalizes this framework in a widely used self-report measure. Specifically, the four AASP quadrants are Low Registration (high threshold, passive response tendency), Sensation Seeking (high threshold, active response tendency), Sensory Sensitivity (low threshold, passive response tendency), and Sensation Avoiding (low threshold, active response tendency). Dunn’s framework is often expressed as four quadrant scores, but conceptually these quadrants arise from the intersection of two underlying neurological threshold and self-regulation strategy. In the present study, we therefore used theory-guided summary scores for low threshold and active strategy to examine these two dimensions directly, while retaining quadrant-based analyses as complementary tests of convergence with the original AASP structure. A key implication of Dunn’s model is that individuals may differ not only in how strongly they respond to sensory input but also in how they regulate its impact over time—that is, in post-input adjustment processes following stimulation.
The pupillary light reflex (PLR) provides a noninvasive window into such processes. The PLR is a multiphasic response initiated by retinal light input, with pupil constriction followed by recovery (redilation)^7,8^. Because a single scalar index (e.g., peak constriction) can miss time-dependent features, the full temporal profile has been proposed as a more appropriate target for analysis^9,10^. As recovery is sensitive to arousal and sustained attention and has been linked to locus coeruleus–noradrenergic activity^11,12^, time-resolved PLR waveforms may be informative for examining post-stimulus temporal differences associated with sensory traits.
Prior work has reported associations between sensory traits and PLR measures in clinical populations, but evidence has been mixed, and generalization to nonclinical samples remains uncertain. For example, Daluwatte et al.^13^ reported that greater sensory dysfunction was associated with smaller constriction amplitude in children with autism spectrum disorder, whereas this association was not observed in their typically developing group. More broadly, autism-related PLR findings have varied across studies^14^.
At least two factors may contribute to these discrepancies. First, analyses that reduce PLR responses to scalar metrics (e.g., latency and peak amplitude) may miss trait-related differences expressed in time-dependent dynamics, particularly during recovery^14,15^. Second, because sensory traits often co-occur with psychological traits (e.g., autism traits, ADHD traits, anxiety, and depression), insufficient separation and statistical control may complicate interpretation^1,16–19^. In addition, PLR waveforms are sensitive to stimulus physical properties (photometric and spectral factors), which should also be considered when interpreting individual differences.
Against this background, we examined how sensory processing traits grounded in Dunn’s model relate to PLR temporal dynamics in nonclinical young adults by treating the PLR as a waveform over 0–3 s after stimulus onset, rather than reducing it to a small set of scalar indices. This window was selected as the primary analysis interval to focus on constriction and early recovery under the present design, while later intervals were evaluated separately in sensitivity analyses. We analyzed the data using generalized additive mixed models, which flexibly represent nonlinear response trajectories and test trait-related modulation as time-varying effects^20^. To strengthen interpretability, we adjusted for co-occurring psychological traits and explicitly accounted for stimulus physical properties, including luminance and spectral characteristics.
Our hypotheses were as follows. In the low-threshold domain, Sensory Sensitivity may be associated with more sustained responses and delayed recovery, whereas Sensation Avoiding may be associated with reduced constriction amplitude and faster recovery. In the high-threshold domain, Low Registration may be associated with relatively rapid convergence toward baseline, whereas Sensation Seeking may be associated with greater response maintenance and altered recovery dynamics. We further expected that these temporal associations would remain evident after adjusting for stimulus physical properties and co-occurring psychological traits.
After preprocessing, five of the 44 participants were excluded based on the predefined inclusion criterion (≥ 6/10 valid trials in the White condition within the 0–3 s analysis window), leaving 39 participants (20 female, 19 male; M = 20.69 years, SD = 1.59) for analysis. Descriptive statistics for psychological and sensory processing traits are shown in Table 1.
The mean AQ score was 20.51 ± 6.70, comparable to the mean reported for Japanese university students (20.7 ± 6.4)^21^. The mean ASRS-v1.1 total score (18 items) was 29.38 ± 8.98, broadly similar to values reported in Japanese university samples (25.3 ± 8.8)^22^. The mean BDI-II score was 11.15 ± 8.13, close to the mean reported for Japanese adults (8.95 ± 6.49)^23^, and fell within the “Minimal” range (0–13) for the BDI-II. Trait anxiety (47.03 ± 11.73) was comparable to Japanese adult norms (men: 47.91 ± 10.40; 46.89 ± 11.28)^24^, whereas state anxiety (38.31 ± 8.85) was somewhat lower than those norms.
The four AASP quadrant scores Low Registration 34.44 ± 6.32, Sensation Seeking 40.08 ± 8.59, Sensory Sensitivity 33.85 ± 7.35, and Sensation Avoiding 33.77 ± 8.07. Compared with Japanese normative data for the 18–34-year group (30.21 ± 7.89, 38.42 ± 8.74, 33.23 ± 8.88, 32.71 ± 8.42, respectively)^25^, these values were all within the typical range. Some participants exceeded screening cutoffs on certain scales (e.g., AQ; n = 3). Because the present study focused on continuous individual differences rather than diagnostic group comparisons, these participants were retained in the analyses.
Results of sensitivity analyses varying the valid-trial thresholds (White ≥ 5/10 and ≥ 7/10) and the analysis windows (0–2, 0–3, 0–4, and 0–5 s) are provided in the Supplementary Materials (Supplementary Table S1 online), and changes in sample size are described in Supplementary Methods S3.1 online. We also examined participant inclusion/exclusion and the association between the valid-trial rate within the 0–3 s window and individual traits (rejection bias) and found no clear evidence of systematic bias (Supplementary Table S2 online).
We used a generalized additive mixed model (GAMM) to examine the association between sensory processing traits (LowTh and Active) and the PPC temporal profile over 0–3 s after stimulus onset. Population-level predicted PPC waveforms by condition are shown in Fig. 1 (point estimates and 95% confidence intervals). In this model, between-condition differences can be expressed not only as shifts in mean level but also as differences in waveform shape. As parametric fixed effects relative to the reference level (White), Red (b = 0.0221, t = 36.10, p < 0.001), Green (b = 0.0050, t = 8.15, p < 0.001), and Blue (b = 0.0117, t = 19.33, p < 0.001) showed overall shifts toward smaller constriction (i.e., more positive PPC). In addition, the time-by-condition smooth term was significant for White, Red, and Blue (ps < 0.001), indicating that waveform shape could differ by condition; in contrast, this term did not contribute for Green (p = 0.999). Main effects of age and sex were not significant (age: p = 0.628; p = 0.577) (Supplementary Table S4 online).Fig. 1Condition-wise predicted PLR waveforms from the main model (M1). Population-level predictions from M1 for each stimulus condition (White, Red, Green, and Blue) over 0–3 s, holding covariates at their means (z = 0). Shaded bands indicate pointwise 95% confidence intervals.
Main effects of sensory traits (i.e., associations with the average response level) were not significant (LowTh: b = 0.00564, t = 0.777, p = 0.437; Active: b = 0.00339, t = 0.426, p = 0.670). In contrast, nonlinear time-by-trait interactions were significant for both traits, indicating that waveform shape could vary as a function of trait levels (LowTh: F(196.60, 235.04) = 17.39, p < 0.001; Active: F(181.68, 219.83) = 32.41, p < 0.001). Under the present conditions and modeling framework, these results suggest that sensory traits may manifest more as differences in temporal patterns than as mean-level differences.
To improve transparency, Fig. 2 shows observed PPC mean traces for descriptively defined High and Low groups of LowTh and Active (median split), whereas Fig. 3 shows the corresponding population-level predicted waveforms with LowTh and Active set to ± 1 SD, and Fig. 4 shows the corresponding difference waveforms (High − Low) (see Supplementary Methods S6 online for interval definitions and interpretation). The observed mean traces in Fig. 2 are provided for descriptive visualization only; inferential conclusions remain based on the continuous-trait GAMM. For the difference waveforms (Fig. 4), the simultaneous confidence intervals included zero across the entire 0–3 s interval, and no localized time window with a reliable difference was identified. Accordingly, the significant time-by-trait smooths are interpreted here as evidence for global differences in fitted waveform shape under the present modeling framework, rather than as evidence for reliable divergence at specific time points. In this context, the smooth interaction terms indicate that the fitted response functions vary with trait level across the analyzed interval, whereas the simultaneous confidence intervals indicate that the present data do not support confident localization of that variation to any specific time window. To illustrate between-participant variability for the same descriptive grouping, participant-level PPC traces are shown in Supplementary Figure S5.Fig. 2Observed PPC mean traces for LowTh and Active groups. Observed proportional pupil change (PPC) waveforms over 0–3 s for descriptively defined High and Low groups of LowTh and Active (median split). Curves show group means across participants, averaged as equal-weight means across the four stimulus conditions (White, Red, Green, Blue). Shaded bands indicate 95% confidence intervals. These plots are provided for descriptive visualization only; inferential conclusions remain based on the continuous-trait GAMM.Fig. 3Trait-modulated PLR waveforms for Dunn’s two axes in the main model (M1). Predicted waveforms for High (+ 1 SD) and Low (− 1 SD) levels of Low threshold and Active strategy, with other predictors fixed at their means (z = 0; equal-weight mean across the four stimulus conditions). Shaded bands indicate pointwise 95% confidence intervals.Fig. 4Trait-modulated PLR difference waveforms (High − Low) for Dunn’s two axes in the main model (M1). Difference curves computed as High (+ 1 SD) minus Low (− 1 SD) predicted waveforms for Dunn’s two axes (Low threshold and Active strategy) in the main model (M1), based on the equal-weight mean across the four stimulus conditions (White, Red, Green, and Blue) with covariates held at their means (z = 0) and random effects excluded. Shaded bands indicate simultaneous 95% confidence intervals across 0–3 s (FWER-controlled). Simultaneous confidence intervals are shown to assess whether the High–Low contrast can be localized to specific time points across the 0–3 s interval; intervals overlapping zero throughout indicate that no reliable localized divergence was identified.
As an aid to interpretation, we computed the scalar summary metrics defined in Supplementary Methods S6.3 online (early/late mean differences based on a common boundary t*, as well as PPCmin, tmin, T50, and their differences). The estimated t* was 0.85 s (95% CI [0.83, 0.87]). For the High − Low difference waveform, the early and late mean differences were 0.0038 (95% CI [− 0.0296, 0.0368]) and 0.0142 (95% CI [− 0.0151, 0.0442]) for LowTh, and 0.0267 (95% CI [− 0.0128, 0.0658]) and − 0.0013 (95% CI [− 0.0348, 0.0315]) for Active. High − Low differences (Δ) in PPCmin, tmin, and T50 also had 95% CIs that included zero (see Supplementary Table S5 online for details).
We also computed functional effect-size measures summarizing the High − Low difference waveform over 0–3 s. For LowTh, the maximum absolute difference was 0.023 (95% CI [0.011, 0.066]), the RMS difference was 0.013 (95% CI [0.006, 0.041]), and |Δ|AUC was 0.034 (95% CI [0.013, 0.118]). For Active, the corresponding values were 0.035 (95% CI [0.020, 0.087]), 0.017 (95% CI [0.010, 0.045]), and 0.042 (95% CI [0.025, 0.120]) (Supplementary Table S5 online). These indices complement the significant time-by-trait smooth terms by summarizing the overall magnitude of waveform-level differences across 0–3 s, but they do not identify a specific time window of reliable divergence.
Robustness was examined using a subject-level cluster bootstrap (200 iterations). The time-by-trait interactions for LowTh and Active were consistently supported across iterations, and the shapes of the difference waveforms broadly matched the primary analysis (Supplementary Fig. S3 online; Supplementary Table S6 online).
Overall, within the scope of the present study, associations with sensory processing traits appeared more as differences in waveform shape over 0–3 s than as differences in average response level.
To examine whether the temporal modulation of PLR waveforms observed in the primary analysis could be explained by co-occurring psychological traits, we fitted a robustness model (M2) that additionally adjusted for psychological covariates (AQ, ASRS, BDI-II, and STAI).
In M2, the added psychological covariates and the main effects of sensory traits did not reach conventional significance levels (BDI-II showed a p = 0.080; other ps ≥ 0.14). In contrast, the nonlinear time-by-trait interactions remained significant after adjusting for psychological covariates (LowTh: F(196.52, 235.02) = 17.39, p < 0.001; Active: F(181.59, 219.69) = 32.43, p < 0.001) (Supplementary Table S4 online). These results suggest that temporal modulations of PLR waveforms associated with sensory processing traits may persist even when the measured psychological traits are simultaneously considered.
To test whether the temporal modulation observed in the primary analysis could be explained by physical properties associated with stimulus conditions (luminance and melanopic DER), we fitted a physical-parameter model (M3) in which categorical condition was replaced by physical indices.
In M3, the linear main effect of melanopic DER was significant, with higher melanopic DER being associated with PPC shifts toward more negative values (b = − 0.00523, t = − 19.80, p < 0.001). In contrast, the linear main effect of Lv was not significant (b = − 0.0449, t = − 0.001, p = 0.999). However, time-varying effects indicating that waveform shape could vary with the level of physical indices were significant for both Lv and melanopic DER (Lv: F(80.67, 92.80) = 185.62, p < 0.001; DER: F(67.82, 83.27) = 111.86, p < 0.001) (Supplementary Table S4 online). Note that across the four conditions, Lv and melanopic DER were not independent (correlation based on condition means, Pearson’s r = − 0.58), which may limit strict separation of their contributions to mean level via linear main effects.
Even after explicitly modeling physical indices, the main effects of sensory traits were not significant (ps > 0.40), whereas nonlinear time-by-trait interactions remained significant (LowTh: F(72.79, 102.66) = 1.98, p < 0.001; Active: F(77.10, 98.94) = 5.56, p < 0.001). Thus, within the present conditions and modeling framework, LowTh- and Active-related temporal modulation is unlikely to be explained solely by between-condition differences in physical properties (luminance and melanopic DER) (Supplementary Fig. S4A-B online).
As supplementary analyses, we examined the association between PLR temporal dynamics and each AASP quadrant score (Low Registration, Sensation Seeking, Sensory Sensitivity, and Sensation Avoiding), adjusting for psychological covariates.
Across all four quadrants, main effects (differences in average response level) were not significant, whereas tensor-product time-by-quadrant interactions were significant for all quadrants (Low Registration: F(184.98, 199.92) = 40.36, p < 0.001; Sensation Seeking: F(161.31, 184.46) = 8.83, p < 0.001; Sensory Sensitivity: F(154.31, 181.72) = 9.46, p < 0.001; Sensation Avoiding: F(157.52, 183.06) = 14.06, p < 0.001) (Supplementary Table S7 online). As in the primary analysis (M1), these results indicate that quadrant effects may manifest not as mean-level differences in PPC but as modulations of waveform shape over the 0–3 s post-stimulus interval.
Figure 5 shows predicted waveforms with each quadrant score set to ± 1 SD, and Fig. 6 shows the corresponding difference waveforms (High − Low) (see Supplementary Methods S6 online for interval definitions and interpretation). For all quadrants, the 95% simultaneous confidence intervals for the difference waveforms included zero across the entire 0–3 s interval, and no localized significant time windows were identified. These quadrant-level results are therefore interpreted as global differences in fitted waveform shape under the model, rather than as evidence for reliable divergence during a specific sub-interval of the response. That is, the interaction terms capture variation in fitted waveform shape across the interval, whereas the simultaneous confidence intervals do not identify a specific portion of the waveform in which that variation can be localized with confidence.Fig. 5Predicted PLR waveforms for AASP quadrant models (S1–S4). Population-level predicted waveforms for High (+ 1 SD) and Low (− 1 SD) levels of each AASP quadrant score (S1–S4), with covariates held at their means (z = 0) and random effects excluded. Waveforms are shown as the equal-weight mean across the four stimulus conditions (White, Red, Green, and Blue). Shaded bands indicate pointwise 95% confidence intervals.Fig. 6Difference waveforms (High − Low) for AASP quadrant models (S1–S4). Difference curves computed as the High (+ 1 SD) minus Low (− 1 SD) predicted waveforms for each AASP quadrant model, based on the equal-weight mean across the four stimulus conditions (White, Red, Green, and Blue), with covariates held at their means (z = 0) and random effects excluded. Shaded bands indicate simultaneous 95% confidence intervals across 0–3 s (FWER-controlled). Simultaneous confidence intervals are shown to assess whether the High–Low contrast can be localized to specific time points across the 0–3 s interval; intervals overlapping zero throughout indicate that no reliable localized divergence was identified.
Scalar summary metrics (early and late mean differences based on a common boundary t*, peak constriction, time-to-peak, half-recovery time, and their differences) all had 95% CIs that included zero (Supplementary Table S5 online). This is consistent with the possibility that trait differences may be difficult to isolate using a single scalar metric, while still emerging in the temporal structure of the full waveform.
To complement this, we computed functional effect sizes summarizing the difference waveforms over 0–3 s. Sensation Avoiding (S4) showed relatively larger differences, with a maximum absolute difference of 0.041 (95% CI [0.016, 0.088]), an RMS difference of 0.026 (95% CI [0.009, 0.056]), and |Δ|AUC of 0.071 (95% CI [0.021, 0.160]) (Supplementary Table S5 online). In contrast, these indices were relatively small for Sensation Seeking (S2).
Overall, the influence of AASP quadrants appears more as modulations of the waveform shape over 0–3 s than as mean-level differences in PPC, suggesting that quadrant traits may be associated with differences in PLR temporal profiles.
The present study examined associations between sensory processing traits grounded in Dunn’s model and the temporal profile of the pupillary light reflex (PLR) in a nonclinical young-adult sample. Across models, trait-related associations were expressed more clearly as differences in fitted waveform shape over the 0–3 s post-stimulus interval than as shifts in average response level. Because the PLR is a multiphasic response—with parasympathetically mediated constriction followed by recovery (redilation) involving parasympathetic withdrawal and sympathetic contributions^7,8^—single scalar indices such as latency or peak constriction inevitably capture only part of the information contained in the full time course waveform. The present findings are therefore best interpreted as model-based, global differences in fitted waveform shape, rather than as evidence for reliable divergence at specific time points.
In this context, the significant time-by-trait smooths should be understood as indicating variation in the fitted response functions across the analyzed interval, whereas the High–Low contrast curves with simultaneous confidence intervals address whether such variation can be localized to a specific time window. Rather than attempting to localize trait effects to specific time points, we summarized overall waveform differences across the 0–3 s interval using functional effect-size measures. This emphasis is consistent with the uncertainty inherent in estimating the timing of maximal differences and with the possibility that peak timing varies across individuals and state factors. Accordingly, the results suggest that sensory processing traits may be associated with global post-stimulus temporal profiles in the PLR, without implying a single fixed time window in which differences must occur.
Among the four Adolescent/Adult Sensory Profile quadrants, the clearest pattern was observed for Sensation Avoiding (S4): higher S4 scores were associated with a systematic shift toward shallower constriction across a broad portion of the waveform. Functional effect-size measures likewise indicated that S4 showed the largest overall differences among the quadrants. In Dunn’s framework, Sensation Avoiding reflects a low neurological threshold combined with a more proactive tendency to reduce sensory input to prevent overactivation^4^. The present findings are therefore consistent with a relatively shallower constriction pattern in individuals with higher S4 scores under the present stimulation conditions, but they do not by themselves establish the physiological basis of that pattern. Any interpretation in terms of autonomic mechanisms should therefore be regarded as tentative and will require future studies specifically designed to separate parasympathetic and sympathetic contributions.
The derived LowTh and Active scores should be understood as theory-guided summaries of Dunn’s two conceptual dimensions, whereas the quadrant-based models preserve the original AASP scoring structure. In this sense, the supplementary quadrant analyses were used to assess whether the waveform-level findings observed in the two-axis analyses were broadly consistent with the original four-quadrant representation.
Within this two-axis representation, the active-strategy dimension also showed waveform modulation toward shallower constriction at higher trait levels. Comparisons with the supplementary quadrant models suggested that this Active-axis pattern aligned most closely with the waveform pattern observed for Sensation Avoiding (S4). Although Sensory Sensitivity (S3; a low-threshold but less proactive pattern) also showed some waveform differences, these appeared smaller than those for S4. More broadly, the present findings suggest that sensory traits may be associated with differences in the temporal profile of the PLR, rather than being fully captured by a single scalar index of response magnitude.
In contrast, the hypothesized pattern for Sensation Seeking (S2)—specifically, delayed recovery or response maintenance—was not clearly supported. One possibility is that the present paradigm did not strongly engage the motivational and sustained engagement features that are often emphasized for Sensation Seeking within Dunn’s framework^4,5^. The stimulus was a brief (100 ms) full-field flash followed by darkness, and pupil dynamics are known to be sensitive to task demands, motivation, and sustained attention^9,12^. Future studies may therefore benefit from using contexts more likely to recruit sustained engagement (e.g., reward-related or affectively salient stimulation) or designs that place greater demands on persistence (e.g., repeated stimulation sequences).
Importantly, time-varying trait associations remained statistically supported in robustness analyses that adjusted for co-occurring psychological traits (AQ, ASRS, BDI-II, STAI) and in a physical-parameter model that explicitly accounted for stimulus physical properties (luminance and melanopic DER). These adjustments reduce the likelihood that the observed time-varying effects are attributable solely to correlated psychological symptoms or to physical differences between conditions. Within the present conditions and analytic framework, the results therefore suggest that sensory processing traits explain unique variance in PLR temporal dynamics beyond these measured covariates.
The present findings may also help organize prior mixed results linking sensory traits to pupillary measures. Many studies have relied on scalar metrics (e.g., peak constriction or latency), which may be insensitive to differences expressed in time-dependent dynamics, particularly recovery-related features^9,10,15^. By modeling the full waveform and accounting for artifact handling and within-trial autocorrelation (AR(1)), the current approach provides a principled framework for detecting trait-related modulation of temporal profiles while avoiding over-interpretation of localized time windows. At the same time, this flexibility comes at the cost of increased model complexity, and the present sample size limits how confidently fine-grained waveform-level effects can be interpreted.
Several limitations should be noted. First, the study used a cross-sectional design; thus, causal direction and mediating mechanisms cannot be established. Second, the sample was limited to a relatively homogeneous group of young adults, and the sample size was modest relative to the flexibility of the GAMM framework used here. Accordingly, although the model was designed to capture waveform-level structure, effect estimates remain uncertain in terms of stability, precision, and interpretability, and the present results should be regarded as provisional pending replication in larger and more diverse samples. Although the sensitivity analyses and subject-level cluster bootstrap provided supportive evidence for the stability of the main waveform-level patterns, they do not eliminate the uncertainty associated with the present sample size and model complexity. Third, the primary analyses focused on the 0–3 s interval to target constriction and early recovery under the present design. Although later post-stimulus dynamics may contain additional information, including possible melanopsin-related influences, the present stimulus set was not optimized to isolate such late responses; this was one reason why later intervals were treated as sensitivity analyses rather than as the primary inferential target. Fourth, although preprocessing excluded blinks/tracking loss and large gaze deviations, residual measurement confounds cannot be ruled out. In particular, we did not correct for pupil foreshortening error^26^, and state factors such as attention and arousal may have influenced the temporal profile^27^. Additional factors such as sleep, caffeine/nicotine intake, medication, ocular surface conditions, dark-adaptation time, and menstrual cycle phase can also influence pupil measurements^28^ and warrant more systematic assessment in future work.
Building on these limitations, the present study yields testable predictions. If trait-related waveform modulation—particularly in recovery-related features—is partly mediated by state factors (attention/arousal), then experimentally manipulating attentional load should systematically alter recovery dynamics and the magnitude of overall waveform differences. If the expression of Sensation Seeking depends on engagement, then adding reward value or affective salience, or using repeated stimulus sequences, should increase sensitivity to trait-related temporal profiles. Such manipulations could help clarify the relative contributions of central adjustment processes and peripheral reflex mechanisms underlying shape-based individual differences in the PLR.
In a nonclinical young-adult sample, the present study found model-based associations between sensory processing traits and the temporal profile of the pupillary light reflex (PLR), with trait-related differences appearing more clearly in waveform shape than in mean-level shifts. These findings suggest that waveform-based analyses may help characterize global temporal response profiles associated with sensory traits. At the same time, the present results are best interpreted as model-based, global differences in fitted waveform shape, rather than as evidence for reliable divergence at specific time points, and they do not establish the underlying physiological mechanisms. Larger and more targeted studies will be needed to clarify the size, timing, and mechanisms of these effects.
We recruited 44 undergraduate and graduate students aged 18–24 years (M = 20.61, SD = 1.57). To contextualize the sample size, we summarized typically developing sample sizes reported in a systematic review of pupillometry studies (mean = 28.37, SD = 18.39)^15^. Participants were recruited from Utsunomiya University, Japan. Sex was recorded as sex assigned at birth (22 female, 22 male). After preprocessing, five participants were excluded based on predefined data-quality criteria, leaving 39 participants for analysis (20 female, 19 male; M = 20.69 years, SD = 1.59). Although this final sample size is within the range commonly reported in related pupillometry studies, it remains modest relative to the flexibility of the present waveform-based GAMM analyses; the resulting effect estimates should therefore be interpreted cautiously. All participants had normal or corrected-to-normal vision and provided written informed consent prior to the experiment. The study was conducted in accordance with the Declaration of Helsinki and was approved by the Ethics Committee for Human Research at Utsunomiya University (approval H22-022).
Before the experiment, participants completed the following Japanese versions of questionnaires. Sensory processing traits were assessed with the Adolescent/Adult Sensory Profile (AASP^6^; Japanese version^25^). ASD traits were assessed with the Autism-Spectrum Quotient (AQ^29^; Japanese version^21^). Adult ADHD symptoms were assessed with the Adult ADHD Self-Report Scale v1.1 (ASRS-v1.1^30^; Japanese version^22^). Anxiety was assessed with the State–Trait Anxiety Inventory (STAI^31^; Japanese version^24^). Depression was assessed with the Beck Depression Inventory-II (BDI-II^32^; Japanese version^23^). Reliability and validity of these Japanese versions have been established. For reference, prior studies reported internal consistency (Cronbach’s α) of 0.81 for AQ^21^, 0.83 for ASRS-v1.1^22^, 0.87 for BDI-II^23^, and 0.91–0.93 for STAI^24^. For the AASP, reliability coefficients have also been reported for each quadrant (0.75–0.83)^25^.
In the analyses, the four AASP quadrant scores (Low Registration, Sensation Seeking, Sensory Sensitivity, and Sensation Avoiding) were used as primary predictors. To adjust for psychological traits that may co-occur with sensory processing traits, AQ, ASRS, BDI-II, and STAI were included as covariates, along with age and sex assigned at birth. Continuous participant-level variables (questionnaire scores and age) were z-standardized across participants.
To operationalize Dunn’s two conceptual dimensions more directly, the four AASP quadrant scores were first z-standardized across participants and then combined into a low-threshold axis (LowTh) and an active-strategy axis (Active). These derived scores were intended as theory-guided summaries of the threshold and strategy dimensions, not as replacements for the original quadrant scores. Accordingly, we also analyzed the four quadrant scores themselves in supplementary models (S1–S4) to examine whether the results converged with the two-axis representation. We then computed the following axis \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ LowTh=(\left(SSen{s}{z}+SAvoi{d}{z}\right)-(LRe{g}{z}+SSee{k}{z}))/2
\usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ Active=(\left(SSee{k}_{z}+SAvoi{d}_{z}\right)-(LRe{g}_{z}+SSen{s}_{z}))/2 $$\end{document} Here, LReg, SSeek, SSens, and SAvoid denote the Low Registration, Sensation Seeking, Sensory Sensitivity, and Sensation Avoiding quadrant scores, respectively. Both axis scores were subsequently z-standardized across participants and used as predictors. ### Apparatus Visual stimuli were presented on a 24.1-inch LCD monitor (ColorEdge CS2420-Z; EIZO Corporation, Ishikawa, Japan; resolution 1920 × 1200 pixels; display area 51.98 × 32.54 cm) at a refresh rate of 60 Hz. Participants were seated in a dark room with head stabilization using a chin rest, at a viewing distance of 57 cm. Room illuminance was confirmed to be 0 lx using a lux meter (CL-200A; KONICA MINOLTA, Inc., Tokyo, Japan). Stimulus presentation and experimental control were implemented in MATLAB R2022b (MathWorks, Inc., Natick, MA, USA) using Psychophysics Toolbox^33–35^, running on a MacBook Pro (13-inch, 2017; Apple, Inc., Cupertino, CA, USA; macOS Mojave 10.14.6). Binocular pupil diameter and gaze position were recorded using a video-based eye tracker (iRecHS2 Ver.660^36^) at 500 Hz. Under optimal conditions, the reported gaze accuracy of the system is 0.028° horizontally and 0.096° vertically. Because the iRecHS2 outputs pupil size in pixel units (radius), we normalized the preprocessed pupil time series relative to baseline and analyzed proportional pupil change (PPC). ### Stimuli The stimulus colors were defined by the display’s RGB values (8-bit; 0–255), with four white [255 255 255], red [255 0 0], green [0 255 0], and blue [0 0 255]. Stimuli were presented as full-screen uniform rectangular patches, subtending 49.0° × 31.9° at a viewing distance of 57 cm. Stimulus duration was 100 ms. The PLR is influenced by the physical properties of the stimulus (photometric and spectral components)^28^. In particular, contributions from intrinsically photosensitive retinal ganglion cells (ipRGCs)—the melanopsin pathway—have been reported to affect sustained constriction and recovery processes^37,38^. Therefore, photometric properties for each condition (luminance Lv, CIE 1931 2° chromaticity coordinates (x, y), and melanopic daylight efficacy ratio (melanopic DER)) were calculated from measured spectra and are shown in Table 2 (see Supplementary Methods S1 online for measurement procedures and calculation details). In the statistical analyses, in addition to modeling the condition as a categorical factor, we also estimated an alternative model (M3) that included Lv and melanopic DER as continuous, condition-level covariates to explicitly adjust for physical influences. Although melanopic DER was quantified and included in the physical-control model, the present stimulus set was not designed to isolate melanopsin-driven late PLR responses, because the color conditions were not constructed as isoluminant or melanopsin-isolating stimuli. Note that black [0 0 0] and gray [128 128 128] were not test stimuli; they were measured as reference values for the background color and fixation-point display conditions before and after stimulation, respectively.Table 2Measured luminance (Lv), CIE 1931 (x, y) chromaticity, and melanopic DER for each condition (mean ± SD).ConditionLv (cd/m^2^)x (CIE1931 2°)y (CIE1931 2°)melanopic DERWhite228.1 ± 0.00.3142 ± 0.00020.3299 ± 0.00010.93765157 ± 0.00074392Red59.2 ± 0.10.6868 ± 0.00060.3095 ± 0.00010.00957944 ± 0.0013346Green152.2 ± 0.10.2063 ± 0.00010.7331 ± 0.00020.69160999 ± 8.92E-05Blue17.2 ± 0.00.1481 ± 0.00000.0584 ± 0.00026.32937268 ± 0.01344843Black0.20 ± 0.000.2572 ± 0.00000.2818 ± 0.00011.23667891 ± 0.00355007Gray49.4 ± 0.00.3141 ± 0.00010.3283 ± 0.00010.94182227 ± 0.00038424melanopic DER was computed in accordance with CIE S 026/E:2018 ^41^ as a dimensionless melanopic-to-photopic ratio obtained by applying the melanopsin action spectrum to each condition’s spectral radiance and normalizing the resulting melanopic component by the photopic component based on the CIE standard luminous efficiency function V(λ) (see Supplementary Methods S1 online for details). ### Procedure The trial sequence is shown in Fig. 7. Each trial began with an instruction screen. Participants were instructed to maintain fixation on the central point, view the upcoming stimulus passively, and minimize blinking and body movement during the trial. The instruction screen was presented for at least 5 s and remained visible until the participant pressed a button (self-paced). After the button press, the instruction screen disappeared and a gray fixation point (0.4° diameter) was presented on a black background. The fixation duration was jittered between 4.8 and 5.0 s to reduce the predictability of stimulus onset. After fixation, a color stimulus (White, Red, Green, or Blue) was presented for 100 ms, followed by a 5-s black screen. The main experiment consisted of 40 trials (10 trials per color condition) preceded by 5 practice trials. Trial order was randomized within participants.  ### Eye-tracking and pupil preprocessing Pupil diameter and gaze position were recorded from both eyes at 500 Hz. Trials were segmented based on the stimulus trigger, with stimulus onset set to 0 ms, creating epochs from − 200 to 5000 ms. Samples were marked as missing when they met any of the following (i) blinks or tracking loss (pupil value of 0 or NaN), (ii) out-of-range values (< 10 px or > 400 px), or (iii) gaze deviation > 3°. For blink/tracking-loss and out-of-range segments, padding of 50 ms before and 150 ms after was applied. Missing segments were interpolated using piecewise cubic Hermite interpolation (PCHIP), and a third-order Butterworth low-pass filter (4 Hz) was applied with zero-phase filtering. Pupil values were normalized as relative change from the baseline mean (− 100 to 0 ms) to compute proportional pupil change (PPC). Trials for which the baseline mean could not be computed were treated as invalid. Rules for trial rejection and binocular combination (window-wise) are provided in Supplementary Methods S2 online. ### Dataset construction and participant inclusion After preprocessing, PPC was combined across eyes on a trial-by-trial basis. When both eyes were valid, the left–right average was used; when only one eye was valid, that eye’s data were used. For gaze time series (visual angle), we did not average across eyes; instead, we selected the best eye for each participant and used data from that eye (definition in Supplementary Methods S2 online). After combination, the time-series data were converted to long format and aggregated with a 10-ms bin width, yielding downsampled data at 100 Hz. Each trial was labeled with the stimulus condition, and participant-level covariates (age, sex assigned at birth, AQ, ASRS, BDI-II, STAI) and the four AASP quadrant scores were merged. Participant inclusion criteria were defined based on the primary analysis window (0–3 s): participants were included if they had at least six valid trials (≥ 6/10) in the White condition (10 trials). This window was chosen to focus the primary analysis on the canonical constriction and early recovery phases of the PLR while maintaining a common window for trial validity and participant inclusion. This criterion excluded five of the 44 participants, leaving 39 for analysis. Sensitivity analyses for the inclusion thresholds (≥ 5/10 and ≥ 7/10), as well as for alternative analysis windows (0–2, 0–3, 0–4, and 0–5 s), and an examination of rejection bias are reported in Supplementary Methods S3 and Supplementary Tables S1–S2 online. ### Reliability of psychophysiological measures To assess measurement consistency of PPC, we computed split-half reliability over the 0–3 s post-stimulus interval. Within each stimulus condition, trials were split into odd and even trials, the correlation (Pearson’s *r*) between the two odd–even mean PPC values (averaged across the 0–3 s interval) was computed, and a Spearman–Brown correction was applied. The resulting reliability coefficient was *r* = 0.89. ### Statistical generalized additive mixed models (GAMM) To examine associations between PLR temporal dynamics and sensory processing traits, we used GAMM. Analyses were conducted in R (version 4.5.2) using the mgcv package (version 1.9–4)^20^. The analysis included 39 participants who met the inclusion criterion (≥ 6 valid trials out of 10 in the White condition). The dependent variable was the pupil change time series (proportional pupil change; PPC) over 0–3 s after stimulus onset (10-ms bins), with constriction defined as negative values. Rather than reducing the PLR to a single summary metric, our modeling framework directly modeled the entire 0–3 s time series to evaluate modulations of waveform shape. Continuous predictors (questionnaire scores and age) were z-standardized across participants (mean 0, SD 1). Smoothing parameters were estimated using fast restricted maximum likelihood (fREML) to reduce overfitting. To improve computational efficiency for the large time-series dataset, models were fitted using the bam() function in the mgcv package with discretization enabled^20,39^. Details of the estimation settings (discretization and basis-dimension *k* choices) and model diagnostics (basis-dimension checks and concurvity) are reported in Supplementary Methods S4 online. To address residual autocorrelation, we introduced a first-order autoregressive (AR(1)) error structure within trials^40^. The autocorrelation coefficient ρ was estimated from residuals of a provisional model without AR(1), and the estimated ρ was then used when fitting the AR(1) models. Details of the procedure and residual diagnostics are provided in Supplementary Methods S5 online. ### Model construction (M1–M3, S1–S4) Across all models, the average post-stimulus PLR waveform was represented as a nonlinear smooth function of time. To account for individual differences, we included both a random intercept capturing between-subject differences in mean level and a subject-specific random smooth (factor-smooth interaction) capturing individual differences in time-varying patterns. The primary inferential target was shape modulation represented by nonlinear time-by-trait interactions. In the primary analysis model (M1), stimulus condition (White, Red, Green, Blue) was included as a fixed effect, with age and sex included as covariates. To allow time profiles to differ across conditions, we included a time-by-condition term, *s*(Time, by = Condition). Sensory processing traits were operationalized as the two axes based on Dunn’s model (low LowTh; active Active), and we included their main effects as well as nonlinear time-by-trait interactions, *ti*(Time, trait). This allowed us to test whether PLR waveform shape varied as a function of traits. In this framework, the two-axis model (LowTh and Active) served as the primary theory-guided representation of Dunn’s conceptual dimensions, whereas the quadrant-based models were used as complementary analyses to assess convergence with the original AASP scoring structure. In the robustness model (M2), we retained the model structure of M1 (condition, time terms, time-by-condition, LowTh and Active main effects and time-by-trait interactions, age, sex, random effects, and AR(1)) and added psychological covariates (AQ, ASRS, BDI, STAI-S, STAI-T). This model was used to reduce the possibility that the time-by-trait effects observed in M1 could be explained by co-occurring psychological traits. To avoid underfitting, we re-estimated the model with increased basis dimensions k for selected smooth terms (*k* settings are reported in Supplementary Methods S4 online). In the physical-parameter model (M3), we did not include stimulus condition as a categorical factor. Instead, we entered condition-level physical indices—luminance (Lv) and melanopic DER—as continuous predictors. Because the influence of physical indices may vary over time, we included varying-coefficient smooth terms, *s*(Time, by = *Lv*) and *s*(Time, by = *DER*), to model time-dependent effects. M3 also included LowTh and Active main effects and time-by-trait interactions (ti terms), age and sex, random effects, and AR(1), allowing us to test whether time-by-trait interactions remained after explicitly adjusting for stimulus physical properties. As supplementary analyses, we estimated models (S1–S4) using each of the four AASP quadrant scores (Low Registration, Sensation Seeking, Sensory Sensitivity, Sensation Avoiding) as the main predictor, with psychological covariates adjusted. In S1–S4, we included stimulus condition and the time-by-condition term, age and sex, random effects, and AR(1). Each model evaluated the main effect of the relevant quadrant score and the nonlinear time-by-quadrant interaction, *ti*(Time, quadrant), without simultaneously including LowTh and Active. ### Inference and model evaluation Fixed effects, covariates, smooth terms, random effects, and basis-dimension (k) settings for the model set (M1–M3, S1–S4) are summarized in Table 3. Procedures for inference and visualization (population-level predicted waveforms and difference waveforms; pointwise and simultaneous confidence intervals), auxiliary scalar summaries for interpretation (peak constriction, time-to-peak, half-recovery time T50, early and late interval mean differences), functional effect-size metrics (maximum absolute difference, RMS difference, and time integral of absolute difference), and supplementary robustness analyses (subject-cluster bootstrap) are detailed in Supplementary Methods S6 online .Table 3Summary of GAMM model specifications (M1–M3 and S1–S4), including smoothing basis dimensions (k) settings.ModelCondition termTime smoothTime-by-conditionTime × Trait (ti)Random effectsk settingsM1CondNames(Time)s(Time, by = CondName)ti(Time, LowTh) + ti(Time, Active)s(Subject) + s(Time, Subject)K_main = 60; K_cond = 20; K_ti = (60,5); K_fs = 10M2CondNames(Time)s(Time, by = CondName)ti(Time, LowTh) + ti(Time, Active)s(Subject) + s(Time, Subject)K_main = 100; K_cond = 40; K_ti = (60,5); K_fs = 20M3–s(Time)–ti(Time, LowTh) + ti(Time, Active)s(Subject) + s(Time, Subject)K_main = 100; K_phys = 100; K_ti = (60,5); K_fs = 20 (CondName excluded)S1–S4CondNames(Time)s(Time, by = CondName)ti(Time, Quadrant)s(Subject) + s(Time, Subject)K_main = 60; K_cond = 20; K_ti = (60,5); K_fs = 10Note. In this table, k denotes the basis dimension for smooth terms in mgcv and is not identical to the effective degrees of freedom (edf). We used K_main as the basis dimension for s(Time_s), K_cond for s(Time_s, by = CondName), K_fs for s(Time_s, Subject_ID, bs = “fs”), and K_ti for ti(Time_s, trait) with k = c(K_ti_time, K_ti_trait). In M3, CondName was not included; instead, condition differences were modeled via Lv and melanopic DER, including their time-varying terms (K_phys). In interpreting the GAMM results, we distinguish between global shape modulation and localized temporal separation. Significant nonlinear time-by-trait interaction terms indicate that the fitted PLR waveform functions vary as a function of trait level across the analyzed interval, but they do not by themselves establish a reliable difference at any single time point. To examine whether such differences could be localized, we also visualized High–Low difference waveforms with simultaneous confidence intervals across 0–3 s. If these intervals include zero throughout, the results are interpreted as supporting global waveform-level differences under the model without identifying a specific time window of reliable divergence. Given the flexibility of these waveform-based GAMMs, the resulting effect estimates should be interpreted primarily at the level of overall waveform structure rather than as highly precise estimates of localized temporal effects. ## Supplementary Information Below is the link to the electronic supplementary material. Supplementary Material 1