A Cross-Recurrence Analysis of the Pupil Size Fluctuations in Steady Scotopic Conditions

Pupil size fluctuations during stationary scotopic conditions may convey information about the cortical state activity at rest. An important link between neuronal network state modulation and pupil fluctuations is the cholinergic and noradrenergic neuromodulatory tone, which is active at cortical level and in the peripheral terminals of the autonomic nervous system (ANS). This work aimed at studying the low- and high-frequency coupled oscillators in the autonomic spectrum (0–0.45 Hz) which, reportedly, drive the spontaneous pupillary fluctuations. To assess the interaction between the oscillators, we focused on the patterns of their trajectories in the phase-space. Firstly, the frequency spectrum of the pupil signal was determined by empirical mode decomposition. Secondly, cross-recurrence quantification analysis was used to unfold the non-linear dynamics. The global and local patterns of recurrence of the trajectories were estimated by two parameters: determinism and entropy. An elliptic region in the entropy-determinism plane (95% prediction area) yielded health-related values of entropy and determinism. We hypothesize that the data points inside the ellipse would likely represent balanced activity in the ANS. Interestingly, the Epworth Sleepiness Scale scores scaled up along with the entropy and determinism parameters. Although other non-linear methods like Short Time Fourier Transform and wavelets are usually applied for analyzing the pupillary oscillations, they rely on strong assumptions like the stationarity of the signal or the a priori knowledge of the shape of the single basis wave. Instead, the cross-recurrence analysis of the non-linear dynamics of the pupil size oscillations is an adaptable diagnostic tool for identifying the different weight of the autonomic nervous system components in the modulation of pupil size changes at rest in non-luminance conditions.

Introduction

The pupil controls the amount of light radiations reaching the retina, by modulating its diameter through the interaction of two muscles under sympathetic-parasympathetic control. The pupil constriction is regulated by the contraction of the iris sphincter muscle receiving parasympathetic innervation mainly through cholinergic fibers. The pupil dilatation is instead related to the contraction of the radial muscle of the iris, under sympathetic control (Loewenfeld and Lowenstein, 1993). Due to the well-known neuroanatomical substrate, the clinical examination of the pupillary light reflex is considered an indicator of the optic nerve conduction, brainstem integrity, vigilance and coma. In recent years, studies in rodents and non-human primates found a tight coupling between pupil size and cortical state even during quiet wakefulness, suggesting a non-luminance mediated system for pupil size variations, associated to neural network oscillations. Studies combining electrophysiology, optical imaging and neural networks modeling, indicated that the link between brain state activity and pupil size is related to the neuro-modulatory effect of the noradrenergic and cholinergic systems (Murphy et al., 2014; Costa and Rudebeck, 2016; Joshi et al., 2016; Eckstein et al., 2017). In this respect, a direct relationship between pupil size and moment-to-moment fluctuations in the activity of noradrenergic neurons of the brainstem locus coeruleus (LC) has been verified (Aston-Jones and Cohen, 2005; Nassar et al., 2012). Other forebrain nuclei and cortical areas connected to LC are activated during spontaneous and event driven pupil size changes (Wang and Munoz, 2015; Joshi et al., 2016) suggesting a circuit for pupil response, linked to arousal, attention and perception systems (Jones, 2004; Naber et al., 2013; Wang and Munoz, 2015; Fazlali et al., 2016; Reimer et al., 2016; Larsen and Waters, 2018). Overall, these studies outline a new role for the pupil size monitoring as a reliable and non-invasive peripheral marker of rapid brain state changes (Hartmann and Fischer, 2014; Schwalm and Jubal, 2017).

From a methodological point of view, a challenge in the analysis of the pupil size variations is the identification of specific patterns that may be representative of changes in the cortical state activity. Different methods have been proposed to assess the pupillary spontaneous oscillations in isoluminant–non-accommodation inducing conditions or in the dark (Lüdtke et al., 1998; Pong and Fuchs, 2000; Zénon et al., 2014; Zénon, 2017). According to the assumptions those methods meet, we distinguish: stationary and linear assumption meeting methods, non-linear assumption meeting methods and non-linear and non-stationary assumption meeting methods. Like other physiological non-stationary signals, under steady stimulation, the pupillary oscillatory signal is expected to show non-linear and chaotic patterns (Poon and Merrill, 1997; Morad et al., 2000; Wilhelm et al., 2001; Merritt et al., 2004; Muppidi et al., 2013; Regen et al., 2013). The non-linear methods assume that the dynamics of the pupil size follow the rules of deterministic chaos rather than a stochastic or linear process (Rosenberg and Kroll, 1999). Common non-linear methods for the analysis of pupillary oscillations imply the use of the Short Time Fourier Transform (Nowak et al., 2008) and wavelets transformations (Henson and Emuh, 2010; Nowak et al., 2013; Reiner and Gelfeld, 2014). These methods assume an underlying stationary signal or require an a priori knowledge of the shape of the single basis wave; assumptions that do not well reflect the pupillary dynamics (Onorati et al., 2016). Among the most recent proposed non-linear and non-stationary meeting methods for the analysis of the pupil oscillations, there are the Hilbert Huang Transform, the EMD (Ruiz-Pinales et al., 2016; Villalobos-Castaldi et al., 2016), and the recurrence plots (Mesin et al., 2013, 2014; Monaco et al., 2014). The Hilbert-Huang transform is a frequency domain transformation, with the advantage of maintaining a good temporal and frequency resolution. Through the EMD, the original signal is split into components with slowly varying amplitude and phase, also known as IMFs. By applying a Hilbert transform to the IMF, instantaneous frequencies are generated as functions of time that give sharp identifications of embedded structures (Barnhart, 2011; Ruiz-Pinales et al., 2016; Villalobos-Castaldi et al., 2016). The RQA consists in taking single physiological measurements, projecting them into multidimensional space by embedding procedures and in identifying correlations that are not apparent in one-dimensional time series. This method provides quantitative indexes related to the number and duration of recurrences of the trajectory of a dynamical system in the phase space (Marwan, 2008; Webber and Marwan, 2015). Then, by applying the cross-recurrence analysis (CRQA) which is a bivariate extension of the RQA, we can investigate the dynamic interactions among the systems modulating pupil size oscillations. The use of CRQA has the advantage to better capture the recurring properties of a dynamic system given by the interaction over time of streams of information (Marwan, 2008; Coco and Dale, 2014). For this purpose, the EMD and CRQA were applied in succession. The main goal of our analysis was the identification of specific frequency components of the oscillatory signal comprised in the range of ANS, that could be quantified by couples of DET and ENT lying within the 95% prediction ellipse. Our result suggests that, in awake healthy subjects at rest, pupils oscillate in darkness with high frequency (HF) and low frequency (LF) components that are in the range of ANS, suggesting a balance between noradrenergic/cholinergic tone. Moreover, the position of the points on the ENT-DET plane seems to be related to the ESS score, and therefore, could give insights into the sleepiness state.

Materials and Methods

Participants

Twenty-six healthy subjects participated to the study (average age 36 ± 13 years old). The participants did not have neurological deficits or serious refractive problems. Moreover, the participants did not assume caffeine in the 2 h preceding the data collection (Wilhelm et al., 2014), and they reported to have slept more than 6 h in the night before the recording (average sleep hours 7.2 ± 0.1). The data collection was performed always between 3 and 6 pm. All subjects gave their written informed consent and the study respected the Declaration of Helsinki and was approved by the local Ethics Committee (Comitato Etico Locale Azienda Ospedaliera Universitaria Senese, EVAlab protocol CEL no. 48/2010).

Experimental Setting

Pupil diameter recordings were performed monocularly with an ASL 504 eye-tracker device (Applied Science Laboratories, Bedford, MA, United States) sampling at a mean frequency of 240 Hz. The remote eye-tracker was placed at 650 mm far from the eye of the participant. The relative position of the subject’s head with respect to the eye-tracker was kept still by mean of a chinrest.

Acquisition Protocol

Prior the data collection, the subjects were administered with an test ESS to investigate their vigilance state. ESS scores less than 11 are normally associated to subjects having normal sleepy state, while ESS scores greater than 11 suggested excessive daytime sleepiness (Parkes et al., 1998). ESS is a common used self-assessment questionnaire for the tiredness evaluation, hence it can turn to be a bias-prone measurement.

All the recordings were performed in a quiet light-controlled environment. To avoid the stimulation of pupillary light reflex, the subject was instructed to look straight for 15 min in a complete dark room (0 lux), similarly to the procedure adopted by Lüdtke et al. (1998). To reduce mental activity and cognitive load, subjects were instructed to try not to think to anything and to relax.

Data Processing

The flow chart of procedure employed to analyze the pupillary frequency balancing between the sympathetic and parasympathetic systems is shown in Figure 1. The pupil diameter data was exported in comma separated values format files and analyzed offline through Matlab (The Mathworks). The signal was de-blinked. Signal instances with the pupil diameter equal to zero were marked as blinks and removed from the signal. The remaining signal was then linearly interpolated. Moreover, machine artifacts introduced by the eye-tracker device due to failures to detect the pupil, were removed using Hampel filtering and low-pass filtered with a cut-off frequency (f₀) of 2 Hz. The Hampel function computes the median of the data within moving windows. The width of the filter window (w) was determined accordingly to the ratio of the sample frequency (f_s) over the cut-off frequency f₀ (Equation 1):

FIGURE 1

Figure 1. The flow chart presents the major procedures adopted for the analysis of the pupil size oscillation, from the data pre-treatment (deblink and artifact removal) and normalization, to the final drawing of the prediction ellipse. Data points of the prediction region in the entropy-determination plane underwent a further classification analysis and a pairwise comparison of the identified clusters was also done.

$$\text{w} \text{=} \text{0}.\text{44}\cdot {\text{f}}_{\text{s}}{\text{/f}}_{\text{0}}(1)$$

The variation of pupil size was computed with respect to a baseline value of the pupil estimated for each participant. Specifically, the baseline value of the pupil diameter signal was determined as the maximum value of the pupil size attained in the first 60 s of the signal in darkness (baseline), when the signal was expected to be more stable. The mean or the median of the pupil size were possible alternative reference values. However, taking the maximum value as reference enabled us to normalize the signal on the basis of a really observed value and to preserve the dynamics of the phenomenon. A baseline-corrected pupil diameter time series was then calculated as the diameter percentage change with respect to the value gathered in the basal condition (Equation 2).

$$\%change=\frac{X_t-Baseline}{Baseline}\cdot 100(2)$$

where X_t is the pupil diameter recorded at time t. The baseline correction provided the removal of inter-subject variability in pupil size percentage change of the pupil diameter signal (Lowenstein et al., 1963; Figure 2).

FIGURE 2

Figure 2. The time course of the pupil % change of an healthy subject. Negative values indicate a restriction of the pupil size with respect to the baseline value.

Data Analysis

A cubic spline interpolation was used for compressing the percentage change time series with a resolution of five data points per second, which satisfied the Nyquist criterion (for the given 2 Hz cut-off frequency). The EMD was applied to the cubic spline interpolation of the percentage change time series in the autonomic frequency band ranging from 0 to 0.45 Hz (Huang et al., 1998). Since we were interested in a global spectral characterization of the IMFs derived from the EMD, the spectral content of the IMFs was estimated through MUSIC algorithms (Schmidt, 1986). The IMFs having most of the power in the autonomic frequency band were retained. The IMFs were then aggregated accordingly to the HF (0.15–0.45 Hz) and LF (0–0.15 Hz) ranges related to the parasympathetic and sympathetic systems activity (Cabrerizo et al., 2014; Figure 3). A CRQA was performed (Marwan and Kurths, 2002; Marwan, 2016) to assess the similarity between the dynamics of the parasympathetic and sympathetic processes by comparing the interaction between the LF and the HF components in the phase space. Three hyper-parameters must be set in the CRQA: EmbDim, TD, and the neighborhood radius (R). A symplectic geometry-based algorithm was used for estimating the EmbDim (Lei et al., 2002). The TD value was chosen as the one within the range (0: w/EmbDim) that maximized the sample entropy of the percentage change of pupil size. A FAN was taken as the neighborhood criterion, such that the cross-recurrence point density had a fixed predetermined value of 20%.

FIGURE 3

Figure 3. The panels represent the HF and LF components extracted from the pupil % change time series of an healthy subject. The IMFs obtained through the application of the EMD technique whose frequency content was inside the range [0.15 – 0.45 Hz] were aggregated and form the HF component of the ANS activity (lower panel), while the IMFs in the range [0 – 0.15 Hz] gave rise to the LF component (upper panel).

Two main parameters from CRQA were considered: the determinism (DET), which quantifies the fraction of periodic structures in the trajectories of the LF and HF dynamics in the phase space, and the entropy (ENT), which is the Shannon entropy of the diagonal line length distribution. Periodic signals are expected to yield high values of DET and small values of ENT (Marwan et al., 2007). To enlarge the sample size enough to apply clustering procedures on the DET-ENT plane and to investigate more carefully for possible highlights of this method on the analysis of the balancing of the sympathetic and parasympathetic systems thorough the analysis of the oscillations of the pupil diameter, we employed the Gaussian-copula simulation approach. Hence, firstly the association among age, ESS, ENT, and DET was measured by the Pearson’s correlation matrix. Then, a Gaussian-copula which maintained the dependence structure was used to generate one-hundred correlated multivariate data of those variables.

The percentage of pupil size change of each of the simulated data was then represented as a point on the DET-ENT plane.

Statistical Analysis

On the simulated dataset the Doornik-Hansen multivariate normality test (Doornik and Hansen, 2008) was performed to verify the null hypothesis that the points in the DET-ENT plane were generated from a bivariate Gaussian distribution. The 95% prediction ellipse was calculated around the mean of observed points in the DET-ENT plane. Equations 3–4 indicate the formula for determining the length of the two semi-axes:

$$a_x=2\cdot \sqrt[2]{{\text{λ}}_2\cdot \frac{\left(n_{obs}-1\right)\cdot n_{var}\cdot f\left(1-\alpha ,n_{var},n_{obs}-n_{var}\right)}{n_{obs}-n_{var}}}(3)$$

$$a_y=2\cdot \sqrt[2]{{\text{λ}}_2\cdot \frac{\left(n_{obs}-1\right)\cdot n_{var}\cdot f\left(1-\alpha ,n_{var},n_{obs}-n_{var}\right)}{n_{obs}-n_{var}}}(4)$$

where a_x and a_y are the major and minor semi-axes of the ellipse, n_var is the number of variables (=2), n_obs is the number of the observations (=100), λ₁ and λ₂ are the eigenvalues (in descending order) obtained from the spectral decomposition of the covariance matrix of ENT and DET, f is the pdf of the F distribution for the given significance α level and degrees of freedom (n_var, n_obs-n_var). The orientation of the ellipse is given (in radians) by the direction of the eigenvector associated to the largest eigenvalue:

$$\text{θ}=atan\left(\frac{v_y}{v_x}\right)(5)$$

where atan is the inverse tangent function, and v_x and v_y are the components of the eigenvector corresponding to the largest eigenvalue. The coordinates of the points [P_x, P_y] laying on the ellipse contour are calculated as follows:

$$P_x=x_c+\left[\frac{a_x}{2}\cdot cos\left(t\right)\cdot cos\left(\text{ϑ}\right)-\frac{a_y}{2}\cdot sin\left(t\right)\cdot sin\left(\text{ϑ}\right)\right](6)$$

$$P_y=y_c+\left[\frac{a_x}{2}\cdot cos\left(t\right)\cdot sin\left(\text{ϑ}\right)+\frac{a_y}{2}\cdot sin\left(t\right)\cdot cos\left(\text{ϑ}\right)\right](7)$$

where x_c and y_c are the coordinates of the center of the ellipse, and t ranges in the interval [0, 2π].

The prediction region can provide the regulatory reference points for assessing if the underlying slow oscillations in the autonomic band of a new observed pupil size time series have the characteristics of a normal pattern.

Afterward, unsupervised clustering through K-means method with two clusters and a L1-norm distance function was applied within the elliptic prediction area. The two clusters were compared in covariance matrices and means vectors. Accordingly, the Box’s M-test was considered for verifying the homogeneity of the covariance matrices, and the Hotelling’s T² test was used for testing the means. The variables age, ESS and % change associated to each cluster were separately compared through the Mann-Whitney unpaired test.

All statistical tests were two-sided and performed on Matlab with a 5% level of significance.

Results

Self-organized adaptive systems like the brain generate complex signals which are inherently non-linear and non-stationary. Furthermore, unstable, weak, and state-dependent phase-locking characterizes the coupling between the biological oscillators (Shockley et al., 2002). Since the couplings between biological signals could also be predominantly transient, the canonical techniques of signal analysis, which basically rely on the assumption of stationary signals, are not appropriate. More importantly, the autonomic control of the spontaneous pupil fluctuations is expected to have non-linear/chaotic dynamics which can be well explored by recurrence analysis methods, whose domain is in the phase-space trajectories (Mesin et al., 2013). For these reasons, we chose the cross-recurrence method to analyze the spectral components of the ANS activity controlling the pupil fluctuations.

The EMD method was applied to the time series of pupil size variations to extract the low and high frequency components of the signal, which were found in the range of the ANS band. In fact, the EMD procedure, which is known to deal with non-linear and non-stationary signals like the pupil size oscillations, is a data driven method that overcomes the limitation of basis function shape typical of the wavelet decomposition method (Gonalves et al., 2007). The CRQA was then performed over the high and low frequency components and two parameters, i.e., entropy (ENT) and determinism (DET), were retained as the major features which quantified the non-linear dynamics of the high- and LF coupled oscillators in the autonomic band.

In Table 1 the sampling distributions of age, ESS scores, ENT, DET and average pupil change are reported. The sample declared a normal level of diurnal drowsiness (ESS: mean = 7.3; SD = 3.7). Five subjects (four of age lower than 30, and one of age greater than 60) reported relatively high ESS scores (>10). The cross-recurrence analysis returned low values both for ENT (mean = 0.92; SD = 0.09) and for DET (mean = 45.28%; SD = 7.46%). The percentage of pupil change in the sample (% change) (mean = -0.16%; SD = 1.91%) indicated an overall slight loss of pupil size with respect to the baseline, but high variability of the fluctuations as well.

TABLE 1

Table 1. Sampling distributions of age, Epworth Sleepiness Scale, entropy, determinism, and average pupil change.

We firstly analyzed the possible association among the observed age, ESS, ENT, and DET. Based on the results, the DET and ENT variables were not significantly correlated to the age and the ESS score of the participants. Instead, a significant correlation between ENT and DET was found (r = 0.58, p = 0.002).

The bivariate distribution of ENT and DET obtained from Gaussian-copula simulated points (Table 2) is depicted in Figure 4, together with the 95% prediction ellipse. The simulated values of ENT (mean = 0.90; SD = 0.09) and DET (mean = 43.66%; SD = 7.43%) were consistent with the values observed in the sample.

TABLE 2

Table 2. Values of age, Epworth Sleepiness Scale, entropy, determinism, and average pupil change generated from a Gaussian-copula.

FIGURE 4

Figure 4. The contour of the ellipse envelopes the region in the entropy-determinism space where new measured combinations of entropy (ENT) and determinism (DET) will fall with a 95% probability. The null hypothesis that the two variables are generated from a bivariate normal distribution was verified through the Doornik-Hansen test. Values of DET in the y-axis have been rescaled in the range [0 – 1]. A K-means algorithm was used to characterize sub-areas within the prediction ellipse and two distinct clusters were found out. The thick black crosses indicate the centres of the clusters. Cluster 2 (squares) exhibited significant higher combination of ENT and DET than cluster 1 (circles). In addition the ESS resulted greater in cluster 2 compared to cluster 1. This finding suggested that the ENT-DET bivariate distribution inherently conveys information about the drowsiness level of the subjects.

The major parameters of the prediction ellipse are displayed in Table 3. The coordinates of the center of the ellipse are the means of the simulated ENT and DET vectors. The axes of the ellipse indicate the magnitude of the inertia along the directions of ENT and DET. The interval estimations of ENT and DET were obtained through Equations 6 and 7. Table 4 displays the expected intervals of determinism for equally spaced intervals (0.05 bits) of entropy.

TABLE 3

Table 3. Parameters of the 95% prediction ellipse in the ENT-DET plane.

TABLE 4

Table 4. Normative intervals of determinism by ranges of entropy.

Through the K-means procedure, two clusters of points were identified within the prediction ellipse and their descriptive statistics is shown in Table 5.

TABLE 5

Table 5. Descriptive statistics of the clusters identified within the prediction ellipse.

The generated ENT-DET values underwent the Doornik-Hansen multinormality test. The hypothesis of bivariate normal distribution was not rejected (DH statistic = 6.97, p = 0.14).

The covariance matrices of the clusters were not significantly different (Box’s M-test = 3.2; p = 0.37). The result of the Hotelling T² test indicated that the bivariate ENT-DET means vectors between the clusters were significantly different (T² = 200.8; p < 0.0001). The two clusters exhibited also significant different ESS scores (U-test = 701.5; p = 0.002), whilst they did not differ in age (U-test = 1073; p = 0.64), nor in % change (U-test = 933.5; p = 0.14).

Discussion

The analysis of the pupil size oscillations is a promising diagnostic tool, enabling improvements in the identification of cortical state changes. Variations of cortical state activity during wakefulness have a strong influence on neural, perceptual and behavioral responses. Pupil diameter varies not only in response to variation of luminance and accommodation, but also during changes in alertness, attention, mental effort and decision making, suggesting a direct link between pupil size variation and cortical state changes (Preuschoff et al., 2011; Nassar et al., 2012; Naber et al., 2013; Alnæs et al., 2014; de Gee et al., 2014). Changes in the cortical state are associated to well characterized variations of the cortical signal frequency. Specifically, in awake rodents the investigation of local field potentials demonstrated the prevalence of LF fluctuations during periods of quiet resting. However, the initiation of locomotion or whisking was related to the suppression of low frequency components and increased high frequency oscillations (Poulet et al., 2012; Eggermann et al., 2014, McGinley et al., 2015b). This transition between slow and fast cortical activity was also observed across cortical regions (Poulet and Crochet, 2019). Electrophysiological studies have revealed that pupil constriction is associated to slow and synchronous cortical responses and inattentive behavior. Conversely, the cortical activation during task engagement or locomotion shows a persistent desynchronized neuronal activity associated to the dilatation of the pupil (Reimer et al., 2014, 2016; McGinley et al., 2015b; Schwalm and Jubal, 2017). Pupil size fluctuations and cortical state variations are modulated by the central noradrenergic and cholinergic pathways. Thus, monitoring pupil dynamics could be a reliable proxy of the changes in cortical states (Reimer et al., 2014, 2016; McGinley et al., 2015a,b). More specifically, the release of acetylcholine (Ach) from the basal forebrain and noradrenaline (NA) from LC have been shown to drive both the state of cortical connectivity and the pattern of the pupil size oscillations also in resting conditions (Reimer et al., 2016; Schwalm and Jubal, 2017). At the peripheral level, both Ach and NA are neurotransmitters of the ANS (parasympathetic and sympathetic systems, respectively) also controlling the pupil diameter. Overall, these premises encourage exploring new and reliable techniques for pupil dynamics monitoring that allow the identification of parameters attributable to NA and Ach modulatory effect in various cortical state changes.

We propose here, a method that can be used as a quantitative measurement of the non-linear dynamics of the pupil fluctuations. We applied a cross-recurrence technique for estimating determinism (DET) and entropy (ENT) features and their distribution, in order to quantify the degree of coupling between the oscillators of the low (LF) and high frequency (HF) components of the pupillary signal. To the best of our knowledge this is the first study on the use of the ENT-DET plane for analyzing the dynamical systems associated to pupil size fluctuation during stationary scotopic visual conditions.

In our cohort of subjects, we observed low levels of determinism (<60%) and entropy (<1). This is consistent with spontaneous physiological signals recorded from healthy subjects, which are expected to be highly complex. Actually, low determinism can be associated to increase in the uncertainty of the signals, and hence to increase in the signal chaotic properties (i.e., complexity). In facts, complex systems are typically highly ordered. Therefore, they tend to preserve low entropy and counteract the second law of thermodynamics (free energy principle). A de-complexification process occurs when free-running physiological signals present sustained loss of complexity. The loss of complexity leads to less ordered states with higher entropy and with stronger coupling of the oscillators controlling the expression of the signal. This degradation in complexity is typically observed in pathological conditions or advanced aging. Therefore, the major result of this study is the identification of a normative elliptical region in the ENT-DET plane for the pupillary oscillators that could be compared with data from group of patients with neurodegenerative diseases. We hypothesize that the occurrence of points outside of the defined elliptical prediction region may signal potential pathological conditions related to alterations in the ANS. As secondary outcome, we observed that, within the elliptical region of confidence, clusters of points with different characteristics of ENT-DET highly differed also in their ESS scores. This finding suggests that the location of the points in the ENT-DET plane can also reveal alterations in the sleepiness state.

Our results indicate that in resting wakefulness conditions, without the influence of light and accommodation, pupil size oscillations are under the effect of a balanced cholinergic/noradrenergic tone. We believe that the employed CRQA-based method may help to lay the groundwork for studying the LF and HF components of the pupil, which may be related to neuronal network state of the brain at rest. Importantly, it consists in a non-invasive procedure that could be easily adopted in clinical context and for diagnostic assessment such as neurodegenerative conditions. Furthermore, this method is adaptable to different experimental conditions (e.g., variations of the visual stimulus, recording during cognitive tasks, etc) provided that the opportune frequency components are dug out from the signal. The joint recording of the pupil size fluctuations along with other physiological signals (e.g., heart rate variability, EEG, etc) would improve the method, since the study of possible time-dependent and/or frequency-related changes in autonomic functions would be facilitated by this integration.

Ethics Statement

The study was approved by the local Ethical Committee Comitato Etico Locale Azienda Ospedaliera Universitaria Senese, EVAlab protocol CEL no. 48/2010.

Author Contributions

All authors conceived and designed the study, critically revised the manuscript, and approved the final version of the manuscript. FR and VS acquired the data. AR, PP, and VS involved in the analysis and interpretation of data, and drafted the manuscript. AR revised the scientific content of the study.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank particularly Dr. Gemma Tumminelli for the help in the recruitment of the participants and for the data collection.

Abbreviations

ANS, autonomic nervous system; CRQA, cross-recurrence quantification analysis; DET, percentage of determinism calculated from the cross-recurrence analysis; EmbDim, embedding dimension, a hyper-parameter to be estimated for the CRQA; EMD, empirical mode decomposition; ENT, entropy calculated from the cross-recurrence analysis; ESS, Epworth Sleepiness Scale; FAN, fixed amount of nearest neighbor; HF, high-frequency component in the ANS spectrum; IMFs, intrinsic mode functions extracted through the EMD; LF, low-frequency component in the ANS spectrum; MUSIC, multiple signal classification algorithm; R, neighborhood radius, a hyper-parameter to be estimated for the CRQA; RQA, recurrence quantification analysis; TD, time delay, a hyper-parameter to be estimated for the CRQA.

References

Alnæs, D., Sneve, M. H., Espeseth, T., Endestad, T., van de Pavert, S. H. P., and Laeng, B. (2014). Pupil size signals mental effort deployed during multiple object tracking and predicts brain activity in the dorsal attention network and the locus coeruleus. J. Vis. 14:1. doi: 10.1167/14.4.1

PubMed Abstract | CrossRef Full Text | Google Scholar

Aston-Jones, G., and Cohen, J. D. (2005). An integrative theory of locus coeruleus-norepinephrine function: adaptive gain and optimal performance. Annu. Rev. Neurosci. 28, 403–450. doi: 10.1146/annurev.neuro.28.061604.135709

CrossRef Full Text | Google Scholar

Barnhart, B. L. (2011). The Hilbert-Huang Transform: Theory, Applications, Development. Iowa: University of Iowa. doi: 10.1146/annurev.neuro.28.061604.135709

CrossRef Full Text | Google Scholar

Cabrerizo, M., Cabrera, A., Perez, J. O., de la Rua, J., Rojas, N., Zhou, Q., et al. (2014). Induced effects of transcranial magnetic stimulation on the autonomic nervous system and the cardiac rhythm. Sci. World J. 2014:349718. doi: 10.1155/2014/349718

PubMed Abstract | CrossRef Full Text | Google Scholar

Coco, M. I., and Dale, R. (2014). Cross-recurrence quantification analysis of categorical and continuous time series: an R package. Front. Psychol. 5:510. doi: 10.3389/fpsyg.2014.00510

PubMed Abstract | CrossRef Full Text | Google Scholar

Costa, V. D., and Rudebeck, P. H. (2016). Previews more than meets the eye: the relationship between pupil size and locus coeruleus activity. Neuron 89, 8–10. doi: 10.1016/j.neuron.2015.12.031

PubMed Abstract | CrossRef Full Text | Google Scholar

de Gee, J. W., Knapen, T., and Donner, T. H. (2014). Decision-related pupil dilation reflects upcoming choice and individual bias. Proc. Natl. Acad. Sci. U.S.A. 111, E618–E625.

Google Scholar

Doornik, J., and Hansen, H. (2008). An omnibus test for univariate and multivariate normality. Oxf. Bull. Econ. Stat. 70, 927–939. doi: 10.1111/j.1468-0084.2008.00537.x

CrossRef Full Text | Google Scholar

Eckstein, M. K., Guerra-Carrillo, B., Singley, A. T. M., and Bunge, S. A. (2017). Beyond eye gaze: what else can eyetracking reveal about cognition and cognitive development? Dev. Cogn. Neurosci. 25, 69–91. doi: 10.1016/j.dcn.2016.11.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Eggermann, E., Kremer, Y., Crochet, S., and Petersen, C. C. (2014). Cholinergic signals in mouse barrel cortex during active whisker sensing. Cell Rep. 9, 1654–1660. doi: 10.1016/j.celrep.2014.11.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Fazlali, Z., Ranjbar-Slamloo, Y., Adibi, M., and Arabzadeh, E. (2016). Correlation between cortical state and locus coeruleus activity: implications for sensory coding in rat barrel cortex. Front. Neural Circ. 10:14. doi: 10.3389/fncir.2016.00014

PubMed Abstract | CrossRef Full Text | Google Scholar

Gonalves, P., Abry, P., Rilling, G., and Flandrin, P. (2007). “Fractal dimension estimation: empirical mode decomposition versus wavelets,” in Proceedings of 32nd IEEE International Conference on Acoustics, Speech and Signal Processing - ICASSP ’07, Honolulu.

Google Scholar

Hartmann, M., and Fischer, M. H. (2014). Pupillometry: the eyes shed fresh light on the mind. Curr. Biol. 24, R281–R282. doi: 10.1016/j.cub.2014.02.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Henson, D. B., and Emuh, T. (2010). Monitoring vigilance during perimetry by using pupillography. Invest. Ophthalmol. Vis. Sci. 51, 3540–3543. doi: 10.1167/iovs.09-4413

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, N. E., Shen, Z., Long, S. R., Wu, M. C., Shih, H. H., Yen, N., et al. (1998). The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. Proc. R. Soc. A 454, 903–995. doi: 10.1098/rspa.1998.0193

PubMed Abstract | CrossRef Full Text | Google Scholar

Jones, B. E. (2004). Activity, modulation and role of basal forebrain cholinergic neurons innervating the cerebral cortex. Prog. Brain Res. 145, 157–169. doi: 10.1016/s0079-6123(03)45011-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Joshi, S., Li, Y., Kalwani, R. M., and Gold, J. I. (2016). Relationships between pupil diameter and neuronal activity in the locus coeruleus, colliculi, and cingulate cortex. Neuron 89, 221–234. doi: 10.1016/j.neuron.2015.11.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Larsen, R. S., and Waters, J. (2018). Neuromodulatory correlates of pupil dilation. Front. Neural Circ. 12:21. doi: 10.3389/fncir.2018.00021

PubMed Abstract | CrossRef Full Text | Google Scholar

Lei, M., Wang, Z., and Feng, Z. (2002). A method of embedding dimension estimation based on symplectic geometry. Phys. Lett. A 303, 179–189. doi: 10.1016/s0375-9601(02)01164-7

CrossRef Full Text | Google Scholar

Loewenfeld, I. E., and Lowenstein, O. (1993). The Pupil: Anatomy, Physiology, and Clinical Applications, Vol. 2. Hoboken, NJ: Wiley-Blackwell. doi: 10.1016/s0375-9601(02)01164-7

CrossRef Full Text | Google Scholar

Lowenstein, O., Feinberg, R., and Loewenfeld, I. E. (1963). Pupillary movements during acute and chronic fatigue a new test for the objective evaluation of tiredness. Invest. Ophthalmol. Vis. Sci. 2, 138–158.

Google Scholar

Lüdtke, H., Wilhelm, B., Adler, M., Schaeffel, F., and Wilhelm, H. (1998). Mathematical procedures in data recording and processing of pupillary fatigue waves. Vis. Res. 38, 2889–2896. doi: 10.1016/s0042-6989(98)00081-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Marwan, N. (2008). A historical review of recurrence plots. Eur. Phys. J. Spec. Top. 164:3. doi: 10.1140/epjst/e2008-00829-1

CrossRef Full Text | Google Scholar

Marwan, N. (2016). CRP Toolbox For Matlab Toolbox. Available at: https://www.academia.edu/14066129/CRP_Toolbox_for_MATLAB (accessed July 16, 2018).

Google Scholar

Marwan, N., and Kurths, J. (2002). Nonlinear analysis of bivariate data with cross-recurrence plots. Phys. Lett. A 302, 299–307. doi: 10.1016/s0375-9601(02)01170-2

CrossRef Full Text | Google Scholar

Marwan, N., Romano, M. C., Thiel, M., and Kurths, J. (2007). Recurrence plots for analysis of complex systems. Phys. Rep. 438, 237–329. doi: 10.1016/j.physrep.2006.11.001

CrossRef Full Text | Google Scholar

McGinley, M. J., David, S. V., and McCormick, D. A. (2015a). Cortical membrane potential signature of optimal states for sensory signal detection. Neuron 87, 179–192. doi: 10.1016/j.neuron.2015.05.038

PubMed Abstract | CrossRef Full Text | Google Scholar

McGinley, M. J., Vinck, M., Reimer, J., Batista-Brito, R., Zagha, E., Cadwell, C. R., et al. (2015b). Waking state: rapid variations modulate neural and behavioral responses. Neuron 87, 1143–1161. doi: 10.1016/j.neuron.2015.09.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Merritt, S. L., Schnyders, H. C., Patel, M., and Basner, R. C. (2004). Pupil staging and EEG measurement of sleepiness. Int. J. Psychophysiol. 52, 97–112. doi: 10.1016/j.ijpsycho.2003.12.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Mesin, L., Cattaneo, R., Monaco, A., and Pasero, E. (2014). “Pupillometric Study of the Dysregulation of the Autonomous Nervous System by SVM Networks,” in Recent Advances of Neural Network Models and Applications. Smart Innovation, Systems and Technologies, Vol. 26, eds S. Bassis, A. Esposito, and F. Morabito (Cha Springer), 107–115. doi: 10.1007/978-3-319-04129-2_11

CrossRef Full Text | Google Scholar

Mesin, L., Monaco, A., and Cattaneo, R. (2013). Investigation of nonlinear pupil dynamics by recurrence quantification analysis. BioMed. Res. Int. 2013:420509. doi: 10.1155/2013/420509

PubMed Abstract | CrossRef Full Text | Google Scholar

Monaco, A., Cattaneo, R., Mesin, L., Fiorucci, E., and Pietropaoli, D. (2014). Evaluation of autonomic nervous system in sleep apnea patients using pupillometry under occlusal stress: a pilot study. Cranio 32, 139–147. doi: 10.1179/0886963413z.00000000022

PubMed Abstract | CrossRef Full Text | Google Scholar

Morad, Y., Lemberg, H., Yofe, N., and Dagan, Y. (2000). Pupillography as an objective indicator of fatigue. Curr. Eye Res. 21, 535–542. doi: 10.1076/0271-3683(200007)2111-zft535

PubMed Abstract | CrossRef Full Text | Google Scholar

Muppidi, S., Adams-Huet, B., Tajzoy, E., Scribner, M., Blazek, P., Spaeth, E. B., et al. (2013). Dynamic pupillometry as an autonomic testing tool. Clin. Auton. Res. 23, 297–303. doi: 10.1007/s10286-013-0209-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Murphy, P. R., O’Connell, R. G., O’Sullivan, M., Robertson, I. H., and Balsters, J. H. (2014). Pupil diameter covaries with BOLD activity in human locus coeruleus. Hum. Brain. Mapp. 35, 4140–4154. doi: 10.1002/hbm.22466

PubMed Abstract | CrossRef Full Text | Google Scholar

Naber, M., Alvarez, G. A., and Nakayama, K. (2013). Tracking the allocation of attention using human pupillary oscillations. Front. Psychol. 4:919. doi: 10.3389/fpsyg.2013.00919

PubMed Abstract | CrossRef Full Text | Google Scholar

Nassar, M. R., Rumsey, K. M., Wilson, R. C., Parikh, K., Heasly, B., and Gold, J. I. (2012). Rational regulation of learning dynamics by pupil-linked arousal systems. Nat. Neurosci. 15:1040. doi: 10.1038/nn.3130

PubMed Abstract | CrossRef Full Text | Google Scholar

Nowak, W., Hachol, A., and Kasprzak, H. (2008). Time-frequency analysis of spontaneous fluctuation of the pupil size of the human eye. Opt. Appl. 38, 469–480.

Google Scholar

Nowak, W., Szul-Pietrzak, E., and Hachol, A. (2013). “Wavelet Energy and Wavelet Entropy as a New Analysis Approach in Spontaneous Fluctuations of Pupil Size Study – Preliminary Research,” in Proceedings XIII Mediterranean Conference on Medical and Biological Engineering and Computing, Seville, 807–810. doi: 10.1007/978-3-319-00846-2_200

CrossRef Full Text | Google Scholar

Onorati, F., Mainardi, L. T., Sirca, F., Russo, V., and Barbieri, R. (2016). Nonlinear analysis of pupillary dynamics. Biomed. Tech. 61, 95–106. doi: 10.1515/bmt-2015-0027

PubMed Abstract | CrossRef Full Text | Google Scholar

Parkes, J. D., Chen, S. Y., Clift, S. J., Dahlitz, M. T., and Dunn, G. (1998). The clinical diagnosis of the narcoleptic syndrome. J. Sleep Res. 7, 41–52. doi: 10.1046/j.1365-2869.1998.00093.x

CrossRef Full Text | Google Scholar

Pong, M., and Fuchs, A. F. (2000). Characteristics of the pupillary light reflex in the macaque monkey: discharge patterns of pretectal neurons. J. Neurophysiol. 84, 964–974. doi: 10.1152/jn.2000.84.2.964

PubMed Abstract | CrossRef Full Text | Google Scholar

Poon, C. S., and Merrill, C. K. (1997). Decrease of cardiac chaos in congestive heart failure. Nature 389, 492–495. doi: 10.1038/39043

PubMed Abstract | CrossRef Full Text | Google Scholar

Poulet, J. F., and Crochet, S. (2019). The cortical states of wakefulness. Front Syst. Neurosci. 8:64. doi: 10.1037/11149-006

CrossRef Full Text | Google Scholar

Poulet, J. F., Fernandez, L. M., Crochet, S., and Petersen, C. C. (2012). Thalamic control of cortical states. Nat. Neurosci. 15:370. doi: 10.1038/nn.3035

PubMed Abstract | CrossRef Full Text | Google Scholar

Preuschoff, K., t Hart, B. M., and Einhauser, W. (2011). Pupil dilation signals surprise: Evidence for noradrenaline’s role in decision making. Front. Neurosci. 5:115. doi: 10.3389/fnins.2011.00115

CrossRef Full Text | Google Scholar

Regen, F., Dorn, H., and Danker-Hopfe, H. (2013). Association between pupillary unrest index and waking electroencephalogram activity in sleep-deprived healthy adults. Sleep Med. 14, 902–912. doi: 10.1016/j.sleep.2013.02.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Reimer, J., Froudarakis, E., Cadwell, C. R., Yatsenko, D., Denfield, G. H., and Tolias, A. S. (2014). Pupil fluctuations track fast switching of cortical states during quiet wakefulness. Neuron 84, 355–362. doi: 10.1016/j.neuron.2014.09.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Reimer, J., McGinley, M. J., Liu, Y., Rodenkirch, C., Wang, Q., McCormick, D. A., et al. (2016). Pupil fluctuations track rapid changes in adrenergic and cholinergic activity in cortex. Nat. Commun. 7:13289. doi: 10.1038/ncomms13289

PubMed Abstract | CrossRef Full Text | Google Scholar

Reiner, M., and Gelfeld, T. M. (2014). Estimating mental workload through event-related fluctuations of pupil area during a task in a virtual world. Int. J. Psychophysiol. 93, 38–44. doi: 10.1016/j.ijpsycho.2013.11.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Rosenberg, M. L., and Kroll, M. H. (1999). Pupillary hippus: an unrecognized example of biologic chaos. J. Biol. Syst. 7, 85–94. doi: 10.1142/s0218339099000085

CrossRef Full Text | Google Scholar

Ruiz-Pinales, J., Salamanca, C., and De Santiago, V. (2016). “Pupillometric-based analysis of central autonomic levels using HHT,” in Proccedings of 2016 International Conference on Mechatronics, Electronics and Automotive Engineering (ICMEAE), Cuernavaca, 14–19.

Google Scholar

Schmidt, R. (1986). Multiple emitter location and signal parameter estimation. IEEE Trans. Antennas Propag. 34, 276–280. doi: 10.1109/tap.1986.1143830

CrossRef Full Text | Google Scholar

Schwalm, M., and Jubal, E. R. (2017). Back to pupillometry: How cortical network state fluctuations tracked by pupil dynamics could explain neural signal variability in human cognitive neuroscience. eNeuro 4:ENEURO.0293-16.2017. doi: 10.1523/ENEURO.0293-16.2017

PubMed Abstract | CrossRef Full Text | Google Scholar

Shockley, K., Butwill, M., Zbilut, J. P., and Webber, C. L. Jr. (2002). Cross recurrence quantification of coupled oscillators. Phys. Lett. A 305, 59–69. doi: 10.1016/s0375-9601(02)01411-1

CrossRef Full Text | Google Scholar

Villalobos-Castaldi, F. M., Ruiz-Pinales, J., Kemper, N. C., and Flores, M. (2016). Biomedical signal processing and control time-frequency analysis of spontaneous pupillary oscillation signals using the Hilbert-Huang transform. Biomed. Signal Process Control 30, 106–116. doi: 10.1016/j.bspc.2016.06.002

CrossRef Full Text | Google Scholar

Wang, C. A., and Munoz, D. P. (2015). A circuit for pupil orienting responses: implications for cognitive modulation of pupil size. Curr. Opin. Neurobiol. 33, 134–140. doi: 10.1016/j.conb.2015.03.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Webber, C. L., and Marwan, N. (2015). Recurrence Quantification Analysis: Theory and best practices, eds L. W. Charles Jr. and M. Norbert (Berlin: Springer). doi: 10.1016/j.conb.2015.03.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilhelm, B., Körner, A., Heldmaier, K., Moll, K., Wilhelm, H., and Lüdtke, H. (2001). Normwerte des pupillographischen schafrigkeitstests fur frauen und manner zwischen 20 und 60 jahren. Somnologie 5, 115–120. doi: 10.1046/j.1439-054x.2001.01156.x

CrossRef Full Text | Google Scholar

Wilhelm, H., Wilhelm, B., Stuiber, G., and Holger, L. (2014). The effect of caffeine on spontaneous pupillary oscillations. Ophthal. Physiol. Opt. 34, 73–81. doi: 10.1111/opo.12094

PubMed Abstract | CrossRef Full Text | Google Scholar

Zénon, A. (2017). Time-domain analysis for extracting fast-paced pupil responses. Sci. Rep. 7:41484. doi: 10.1038/srep41484

PubMed Abstract | CrossRef Full Text | Google Scholar

Zénon, A., Sidibé, M., and Olivier, E. (2014). Pupil size variations correlate with physical effort perception. Front. Behav. Neurosci. 8:286. doi: 10.3389/fnbeh.2014.00286

PubMed Abstract | CrossRef Full Text | Google Scholar