Development of a Unified Cardiovascular-Pupillary Model for Interpreting Pupil Size Variability as an Autonomic Marker

Abstract

Background: Pupil size variability (PSV) has emerged as a potential non-contact indicator of autonomic nervous system (ANS) function; however, its physiological origins remain unclear. This study aims to develop a computational model to investigate the physiological foundations of PSV and assess how frequency-domain indices reflect cardiovascular autonomic balance. Methods: We integrated a well-established cardiovascular regulation model with a biomechanical pupillary muscle plant. The model simulates PSV alongside heart rate variability (HRV) by transmitting respiratory and baroreflex inputs through an indirect neural pathway to the pupillary muscles. Frequency-domain analyses were conducted to compare simulated PSV and HRV across different autonomic states. Results: Simulations suggest that PSV arises from respiratory and baroreceptor inputs, with its classical range-nonlinearity (RNL) property emerging naturally from iris biomechanics. The model reproduces key physiological behaviors, including inspiration-linked dilation and parasympathetic modulation. Frequency-domain analyses reveal low- and high-frequency components in PSV that are similar to those found in HRV. However, the magnitudes of these PSV components depend heavily on the underlying autonomic state and the mean pupil size. Conclusion: These findings provide a mechanistic framework for interpreting PSV as a cardiovascular autonomic biomarker. This ultimately supports the development of non-invasive, wearable, and eye-based ANS monitoring systems.

Keywords

autonomic nervous system heart rate variability pupil size variability pupillary muscle plant biomechanical model

Introduction

The world’s population is ageing rapidly, marking an unprecedented and irreversible demographic shift. According to the World Health Organization, the number of adults aged 60 years or older will rise from 1 billion in 2020 to 1.4 billion by 2030, reaching 2.1 billion by 2050.¹ Countries worldwide are adapting their healthcare and social systems to support the ageing populations. The prevalence of cardiovascular diseases (CVDs) increases with age. Projections indicate substantial growth in CVD cases as the number of older adults expands.² CVDs remain the leading global cause of death, responsible for 19.8 million deaths in 2022, 32% of all deaths, with heart attacks and strokes accounting for 85% of the cases.³ This also results in huge financial burden for the society.

Continuous health monitoring has become essential for effective CVD management, and physiological variability measures play a key role. Heart rate variability (HRV), derived from beat-to-beat fluctuations in cardiac intervals, is a well-established indicator of autonomic function.⁴ Reduced HRV reflects autonomic impairment and is associated with cardiac mortality; it also predicts post-myocardial infarction outcomes and can reveal autonomic imbalance that precedes dangerous arrhythmias.⁴ Blood pressure variability (BPV) is an equally important predictor of cardiovascular risk.⁵ Long-term BPV is linked to higher risks of cardiovascular mortality, coronary heart disease, and stroke because it reflects instability in vascular and autonomic regulation. Power spectral analysis of HRV and BPV is a widely used method for assessing autonomic regulation.⁶ Their spectra contain two key bands: a high-frequency (HF) component (0.15-0.50 Hz) indicating parasympathetic activity, and a low-frequency (LF) component (0.04-0.15 Hz) reflecting combined sympathetic and parasympathetic influences. The LF/HF ratio serves as an index of sympathovagal balance. These frequency-domain measures provide practical markers of autonomic dynamics and baroreflex functions, and have become essential for detecting autonomic dysfunction and predicting cardiovascular risk, as reduced HF power or elevated LF/HF ratios are associated with higher cardiac morbidity and mortality. Continuous monitoring of HRV and BPV in conventional mobile health systems typically requires sensors to be attached at standardized body locations. Ambulatory electrocardiogram (ECG) uses chest electrodes, while photoplethysmogram (PPG) rely on fingertip, earlobe, or forehead probes. Each monitoring session usually involves manual setup, which can be difficult for older adults living alone. Moreover, many users find these sensors intrusive or uncomfortable, limiting their acceptance in daily life. Measuring the pupillary dynamic known as pupil size variability (PSV) offers a promising and unobtrusive method for monitoring autonomic function, as it can be captured using camera-equipped smart eyewear. PSV refers to the spontaneous fluctuations in pupil size that occur without changes in lighting or visual focus. Although once regarded as simple Gaussian noise arising from the pupillary muscles,^7,8 subsequent studies showed that PSV actually embeds respiratory rhythms⁹ and frequency components similar to those found in HRV and BPV.¹⁰ These shared spectral features arise from the common cardiovascular regulatory signals originating in the brainstem and coupled to the pupil via parasympathetic and sympathetic pathways. Therefore, PSV has been proposed as a non-invasive marker of autonomic nervous system (ANS) activity. However, early efforts to apply conventional HRV and BPV frequency-domain techniques to PSV produced inconsistent results, highlighting the need for deeper investigation into the physiological mechanisms underlying PSV.^11,12

This paper presents a computational model that integrates pupillary biomechanics with cardiovascular regulatory mechanisms to simulate PSV. The model is designed to investigate the physiological foundations of PSV and to support its potential use as an indirect, eye-based measure of ANS function. The main objective is to simulate PSV under different ambient brightness conditions and assess how effectively the resulting frequency-domain indices reflect cardiovascular variability and autonomic balance. The insights gained from this work will help advance PSV as an easily accessible biomarker for detecting cardiovascular abnormalities.

Related Works

Across decades of development, computational modeling of pupillary dynamics has evolved from simple linear descriptions to more sophisticated biomechanical and stochastic frameworks. Early works by Day and Dennison provided foundational models. Day et al proposed modeling the pupillary light reflex (PLR) loop as a second-order linear system,¹³ while Dennison introduced a parametric sphincter–dilator model grounded in empirical tension measurements from stimulated cat iris muscles.¹⁴ Semmlow and Stark further proposed a biomechanical model, which comprised spring and damper elements, to simulate PLRs.¹⁵ Hansmann expanded the physiological basis by detailing iris muscle mechanics, autonomic pathways, and key nonlinear behaviors.¹⁶ Building on these foundations, Usui and Stark developed a nonlinear stochastic model that reproduced random fluctuations with a range-nonlinearity (RNL) property.¹⁷ They presumed this behavior as Gaussian noise originating within the pupillary control system. RNL refers to fluctuations that are largest at mid-range pupil sizes and minimal at the extremes of constriction or dilation. Sun et al extended the systems perspective by linking pupil dynamics to broader visual and neurological processes, incorporating adaptive and nonlinear elements.¹⁸ More recently, Denison et al introduced a convolution-based model for explaining pupillary response to rapid sequences of visual events.¹⁹ Although these prior studies demonstrated a progression from simple empirical models to systems-level representations, they only focused on light- or visually-evoked responses, and did not address ANS-related PSV, which was generally regarded as noise rather than a physiological signal.

Physiological Mechanisms

This section introduces the physiology of pupillary dynamics and outlines the physiological mechanisms on which our computational model is based. Light enters the eye through the pupil, regulated by two antagonistic iris muscles: sphincter pupillae and dilator pupillae. As illustrated in Figure 1, the sphincter pupillae receives parasympathetic input via the oculomotor nerve (CN III). Preganglionic fibers from the Edinger–Westphal (EW) nucleus synapse in the ciliary ganglion (CG), where postganglionic fibers trigger contraction of the sphincter pupillae. Conversely, the dilator pupillae is controlled by a three-neuron sympathetic pathway originating in the posterior hypothalamus, as shown in Figure 2. These fibers descend to the ciliospinal center, synapse in the superior cervical ganglion, and innervate the dilator pupillae via long ciliary nerves.

Figure 1.

Parasympathetic innervation pathway from EW nucleus, through CN III, ciliary ganglion and short ciliary nerve, to the sphincter pupillae

Figure 2.

Sympathetic innervation pathway from the spine to the dilator pupillae

Previous studies have established that PSV reflects both respiratory and baroreceptor influences on the pupil size. In anesthetized cats, Borgdorff showed that inspiration produced pupil dilation and expiration produced constriction, with PSV persisting during apnea in synchrony with phrenic discharge after a 200ms to 1000ms delay and independent of sympathetic pathways.²⁰ Human studies by Ohtsuka later confirmed inspiratory dilation, expiratory constriction, and a proportional relationship between PSV amplitude and tidal volume.⁹ Censi et al further demonstrated that carotid baroreceptor stimulation mimicking elevated blood pressure (BP) induces pupil constriction.²¹

The origin of PSV involves indirect neural pathways projecting from the baroreceptors and brainstem’s respiratory center to the EW nucleus in the midbrain. Chen et al showed that midbrain neurons fire rhythmically with respiration, with a delay of more than one second from phrenic onset.²² More recently, Komisaruk and Bianca identified a vagal afferent pathway in rats in which stimulation produced pupil dilation through inhibition of EW output, and demonstrated projections from the nucleus tractus solitarius (NTS) to the parabrachial nucleus (PBN) and onward to the EW nucleus.²³ Together, these studies showed that the NTS is activated by increase in BP and during expiration, and inhibited by decrease in BP and during inspiration.

Consequently, when BP rises or during expiration, the NTS–PBN–EW–CG pathway becomes more active, leading to pupil constriction. When BP falls or during inspiration, this pathway is suppressed, resulting in pupil dilation. Figure 3 summarizes these neural pathways and explains how respiratory and BP fluctuations relate to PSV. Therefore, our computational model is developed under the following three assumptions: i.) PSV has two sources - centrally mediated respiratory activity and baroreceptor-driven activity. Both sources are transmitted from the NTS to the pupillary plant through the NTS–PBN–EW–CG pathway with a delay of more than 1 second; ii.) PSV primarily reflects modulation of the sphincter, with reciprocal innervation between the sphincter and dilator muscles; iii.) Because PSV is identical in both eyes, a single-eye pupillary plant model is sufficient.

Figure 3.

Neural projection pathways underlying PSV showing projections from the NTS to the pupillary plant via the PBN and EW nuclei. Increases in BP cause pupil constriction, whereas decreases cause pupil dilation

Model Development

To our best knowledge, this is the first computational model to describe PSV as a phenomenon mediated by cardiovascular regulatory signals. Implemented in MATLAB and developed based on established neural pathway physiology and experimental findings from previous studies, our model can simulate PSV, HRV, and baroreflex dynamics. For physiological relevance and expandability, the model consists of two parts: a cardiovascular regulation part and a biomechanical pupillary muscle plant.

Cardiovascular Regulation System

The cardiovascular regulation part of the model is adapted from the well-established Seidel–Herzel model, which is widely used in biomedical engineering and computational biology for studying cardiovascular dynamics and autonomic regulation.²⁴ The Seidel–Herzel model is chosen because it provides specific representations of baroreceptor and respiratory activities - the primary sources of PSV. As illustrated in Figure 4, blood pressure changes are detected by the baroreceptors; the resulting afferent signals are transmitted to the medulla, which then generates the sympathetic and parasympathetic outputs for regulating the heart and blood vessels. The model employs delay-differential equations to capture conduction delays in both neural and mechanical pathways. Additionally, medullary sympathetic and parasympathetic outputs are modulated by inputs from higher cerebral regions and respiratory center.

Figure 4.

Overview of the main components of the Seidel–Herzel model for simulating heart rate and blood pressure regulation

The interactions within the model are summarized as follows. Further details can be found in the works by Seidel and Herzel.²⁴ In response to fluctuations in BP, the baroreceptors activities signal is v_b(t):

v_{b} (t) = k_{1} (p (t) - p^{(0)}) + k_{2} \frac{d p (t)}{d t}

(1)

where p(t) is the BP, p⁽⁰⁾ is a constant pressure threshold, and k₁ and k₂ are fixed coefficients characterizing baroceptors’ response. The medulla receives v_b(t) and respiratory activities inputs, and generates medullary sympathetic and parasympathetic activities. Medullary sympathetic activity is expressed as:

v_{s} (t) = \max (0, {v_{s}}^{(0)} - k_{s}^{b} v_{b} (t) + k_{s}^{r} | \sin (π f_{r} t + Δ ϕ_{s}^{r}) |)

(2)

where v_s⁽⁰⁾ is resting sympathetic tone, k_s^b v_b(t) is inhibitory influence from baroreceptors, and k_s^b is coupling between baroreceptors and afferent sympathetic nerves. The sinusoidal term represents modulation by the respiratory neurons having f_r and Δϕ_s^r as respiratory frequency and phase, respectively. Similarly, medullary parasympathetic activity is expressed as:

v_{p} (t) = \max (0, {v_{p}}^{(0)} + k_{p}^{b} v_{b} (t) + k_{p}^{r} | \sin (π f_{r} t + Δ ϕ_{p}^{r}) |)

(3)

where v_p⁽⁰⁾ is resting parasympathetic tone, v_b(t) is baroreceptor activities, and k_p^b is coupling between baroreceptors and afferent parasympathetic nerves. Contractility of the heart is mainly influenced by the sympathetic nerves. The cardiac noradrenaline concentration is denoted as c_cNa(t). Its dynamics is described as:

\frac{d c_{c N a} (t)}{d t} = - \frac{c_{c N a} (t)}{τ_{c N a}} + k_{c_{c N a}}^{s} v_{s} (t - θ_{c N a})

(4)

where θ_cNa is the time delay of sympathetic activity between the medulla and the heart, and k_CcNa^s is the coupling between medulla and efferent sympathetic nerves to the heart. c_cNa(t) is then applied in the equations for describing heart contractility:

S^{'} (i) = S^{(0)} + k_{S}^{c} c_{c N a} (t) + k_{S}^{t} T (i - 1)

(5)

S (i) = S^{'} (i) + [\hat{S} - S^{'} (i)] \frac{{[S^{'} (i)]}^{n_{s}}}{{[S^{'} (i)]}^{n_{s}} + {\hat{S}}^{n_{s}}}

(6)

where S(i) is the contractility of the i^th heartbeat, T(i-1) is the duration of the heart cycle immediately preceding the systole; and and n_s are constants of a threshold function. Phase velocity of the heart’s sinus node is expressed as:

\frac{d ψ (t)}{d t} = \frac{1}{T^{(0)}} f_{s} (t) f_{p} (t)

(7)

where f_s and f_p are the sympathetic and parasympathetic influences on the phase velocity, respectively, and T⁽⁰⁾ is a time constant²⁴. Value of the sinus node phase, ψ, increases from 0 to 1 within a heart period. When it reaches 1, Ψ is reset to 0, contractility S_i is added to the diastolic pressure to give in a new systolic pressure, and a new heartbeat is generated. The BP, denoted by p(t) is described by two separate equations. During systole, of fixed duration of τ_sys, p(t) is described by Equation (8):

p (t) = d (i - 1) + S (i) \frac{t - t_{i}}{τ_{s y s}} \exp {1 - \frac{t - t_{i}}{τ_{s y s}}}

(8)

where d(i-1) is the diastolic pressure preceding the systole onset at beat time t_i, and τ_sys is duration of systole. During diastole, dynamics of p(t) is described by:

\frac{d p (t)}{d t} = - \frac{p (t)}{τ_{v} (t)}

(9)

where τ_v(t) is a fixed duration constant.

NTS-Pupillary Plant Pathway

Figure 5 shows the overall computational model for simulating pupillary dynamics. The model takes respiration frequency, autonomic tone and stimulus, and light condition and stimulus as inputs; and outputs HR, BP, and pupillary dynamics signals accordingly. The key connection between the cardiovascular regulation system and pupillary muscle plant is medullary parasympathetic activity, denoted as v_p(t). It represents neural projection via the NTS–PBN–EW–CG pathway. v_p(t) consists of a static part v_{p_stat} and dynamic part v_{p_dyn}(t) as shown in Equations (10) and (11), respectively.

v_{p_s t a t} = {v_{p}}^{(0)}

(10)

v_{p_d y n} (t) = k_{p}^{b} v_{b} (t) + k_{p}^{r} | \sin (π f_{r} t + Δ ϕ_{p}^{r}) |

(11)

Figure 5.

Functional blocks, inputs, and outputs of the overall computational model, which simulates pupillary dynamics in response to light stimuli and cardiovascular variabilities

While our current focus is on the parasympathetic-dominant fluctuations, the sympathetic pathway is explicitly modeled and contributes to the overall dynamics through reciprocal innervation. As illustrated in Figure 5, the pupil’s sphincter innervation inputs are associated with the parasympathetic activities v_{p_stat} and v_{p_dyn}(t) through the relationships:

e_{s_p_s t a t} = K_{s_s t a t} v_{p_s t a t}

(12)

e_{s_p_d y n} (t) = K_{s_d y n} v_{p_d y n} (t)

(13)

where and e_s__p__stat and e_s__p__dyn(t) are the static and dynamic sphincter innervations related to autonomic condition and PSV-associated fluctuations, respectively, while the constants K_{s_stat} and K_{s_dyn} represent coupling between NTS and sphincter. The overall sphincter innervation input comprises static part, e_s_ _stat, and dynamic part, e_s_ _dyn(t), as follows:

e_{s_s t a t} = {\begin{cases} e_{s_f_s t a t} + e_{s_p_s t a t} & , \\ 1 & , \end{cases} \begin{array}{l} e_{s_f_s t a t} + e_{s_p_s t a t} < 1 \\ e_{s_f_s t a t} + e_{s_p_s t a t} \geq 1 \end{array}

(14)

e_{s_d y n} (t) = e_{s_f_d y n} (t - t_{D_s f}) + e_{s_p_d y n} (t - t_{D_s p})

(15)

e_s__f__stat and e_s__f__dyn(t) are the static and dynamic terms related to response to ambient brightness condition and flash stimulation, respectively. t_{D_sf} and t_{D_sp} are the time delays for pupil flash response and PSV. Reciprocal relationships between the sphincter and dilator innervation are assumed. The expressions for the static and dynamic innervation inputs at the dilator is shown in Equations (16) and (17), respectively:

e_{d_s t a t} = λ_{s t a t} (1 - e_{s_s t a t})

(16)

e_{d_d y n} (t) = - λ_{d y n} e_{s_d y n} (t)

(17)

where λ_{s_stat} = λ_{d_stat} = 1 and 0 ≤ λ_dyn ≤ 1.

Pupillary Muscle Plant

The biomechanical pupillary muscle plant within our model is adopted from the framework proposed by Usui et al. The model parameters were determined from experiments in which smooth muscle tensions were measured under various muscle lengths and stimulation intensities.^25-27 As detailed in Figure 6, the overall innervation inputs, e_s(t) and e_d(t), drive the corresponding input-dependent signals, g_s(t) and g_d(t), which then modulate the pupil size according to tension relationships. Static parts of the input-dependent signals are derived by multiplying the corresponding innervation inputs by constant coupling coefficients, γ_{s_stat} and γ_{d_stat,} as follows:

g_{s_s t a t} = γ_{s_s t a t} e_{s_s t a t}

(18)

g_{d_s t a t} = γ_{d_s t a t} e_{d_s t a t}

(19)

where g_{s_stat} and g_{d_stat} are the static sphincter and dilator, input-dependent terms, respectively. g_{s_dyn}(t) and g_{d_dyn}(t), are related to the dynamic innervation inputs e_{s_dyn}(t) and e_{d_dyn}(t) according to the second-order differential Equations (20) and (21), as follows:

\frac{d^{2} g_{s_d y n}}{d t^{2}} + (2 α_{s} + β_{s}) \frac{d g_{s_d y n}}{d t} + α_{s} (α_{s} + β_{s}) g_{s_d y n} = β_{s} e_{s_d y n} (t - t_{D s})

(20)

\frac{d^{2} g_{d_d y n}}{d t^{2}} + (2 α_{d} + β_{d}) \frac{d g_{d_d y n}}{d t} + α_{d} (α_{d} + β_{d}) g_{d_d y n} = β_{d} e_{d_d y n} (t - t_{D d})

(21)

where the parameters α_s, β_s, α_d, and β_d, t_Ds, and t_Dd are constants obtained from previous empirical studies about impulses responses of iris muscles.¹⁴

Figure 6.

Biomechanical model of the pupillary plant consisting of the dilator and sphincter muscle pair. Pupil size, x, is determined by the combined P_pd, P_ad, P_ps, P_as, and P_v, representing the passive and active tensions of the dilator and sphincter muscles and the viscous element. The innervation inputs, e_d(t) and e_s(t), generate input-dependent signals g_d(t) and g_s(t), which produce muscle tensions according to Carlson’s length–tension relationship

As shown in Figure 6, the acceleration of pupil size x(t) is represented as the sum of the active and passive tensions in the sphincter and dilator muscles. These tensions are calculated by substituting g_s(t) and g_d(t), into the Carlson’s length–tension relationship, which characterizes how smooth muscle tension varies with innervation level and muscle length.²⁷ The parameters in the relationship were empirically derived from results of experiments conducted by Dennison and Usui.^14,27 The pupil size and its rate of change are dependent on the tensions according to Equation (22):

m \frac{d^{2} x (t)}{{d t}^{2}} = P_{p d} + P_{a d} - P_{p s} - P_{a s} - P_{v}

(22)

where x(t), m, P_pd, P_ad, P_ps, P_as, and P_v are pupil size, mass, dilator passive tension, dilator active tension, sphincter passive tension, sphincter active tension, and viscous element, respectively.

Simulation Experiments

The model framework described in the previous sections was implemented in MATLAB using forward Euler method with 2-ms time intervals. The model simulated cardiovascular dynamics and PSV under different scenarios.

Simulation of Cardiovascular Dynamics

To test the model’s ability to synthesize cardiovascular dynamics, initial simulations involved substituting constant and time-varying parameters into Equations (1) through (9) in iterations, each corresponding to increment of 2ms of time, from 0s to 240s. Values of constant parameters in the model were as follows.²⁴ BP parameters: p⁽⁰⁾ = 50.0, τ_sys = 0.125, p_initial(t) = 70mmHg, d_initial= 60mmHg; baroreceptor activity parameters: k₁ = 0.02 mmHg^-1, k₂ = 0.00125 mmHg^-1; respiration parameters: k^r_s = 0.1, k^r_p = 0.1, Δϕ_s^r = 0, Δϕ_p^r = 0; medullary autonomic activities parameters: v_s⁽⁰⁾ = 0.8, v_p⁽⁰⁾ = 0, k^b_s = 0.7, k^b_p = 0.3; transmitter concentration parameters: k^s_CcNa = 1.2, τ_cNa = 2.0s; cardiac contractility parameters: S⁽⁰⁾ = 25mmHg, k^c_s = 40mmHg, k^t_s = 10mmHg s^-1, = 70.0mmHg. Respiration rates of f_r = 0.2Hz and f_r = 0.3Hz were chosen.

The model was able to simulate cardiovascular variabilities, baroregulatory signals, and respiratory sinus arrythmia (RSA) successfully. The results were comparable with those reported in previous works.²⁴ As shown in Figures 7 and 8, the synthesized signals included BP, baroreceptors activities, medullary sympathetic activity, medullary parasympathetic activities, cardiac noradrenaline concentration, cardiac vagal activity, respiratory oscillator, and heart period, denoted as p(t), v_b(t), v_s(t), v_p(t), c_cNa(t), c_vNa(t), ψ(t), and T(i), respectively. The heart period oscillated at the respiratory frequency, corresponding to Mayer waves.

Figure 7.

Results from simulation session with respiration rate of 0.2Hz. Top to bottom trace: p(t), v_b(t), v_s(t), v_p(t), c_cNa(t), c_vNa(t), ψ(t), and T(i)

Figure 8.

Results from simulation session with respiration rate of 0.3Hz. Top to bottom trace: p(t), v_b(t), v_s(t), v_p(t), c_cNa(t), c_vNa(t), ψ(t), and T(i)

Simulation of Pupil Size Variability (PSV)

To simulate PSV, the medullary sympathetic and parasympathetic activity signals, denoted as v_s(t) and v_p(t), were first synthesized using the cardiovascular regulation part of the model. v_p(t) was then input into the pupillary plant, which subsequently produced the corresponding innervation inputs e_{s_stat}, e_s__dyn(t), e_{d_stat}, e_d__dyn(t), input-dependent signals g_{s_stat,} g_{s_dyn}(t), g_{d_stat}, g_{d_dyn}(t), and pupillary dynamics accordingly to Equations (12) to (22). The simulations were conducted using model parameters as listed in Table 1, and across scenarios of increasing respiratory influence (Sessions F1 to F5) and increasing parasympathetic activity (Sessions G1 to G5) as summarized in Table 2. Figures 9 and 10 show twenty seconds of the simulated signals, including PSV, v_p(t), BP, heartbeat period, and respiratory activity under a fixed ambient brightness condition of e_{s_f_stat} = 0 in Sessions F2 and F4, respectively. Similarly, Figures 11 and 12 show the simulated signals under the same lighting condition for Sessions G2 and G4, respectively. It can be seen that, each inspiration was followed by pupil dilation, and each expiration by pupil constriction. Delay between inspiration onset and PSV-associated pupil dilation was about 1 second. These behaviors agree with the observations made in previous animal and human studies.^{9,10,20,21,28}

Table 1.

Model Parameters for Simulating PSV

Parameter	Value (range)	Unit
t_{D_sp}	0.5	s
t_{D_dp}	0.9	s
t_{D_sf}	0.202	s
t_{D_df}	0.705803	s
λ_{dyn_f}	0.989052
λ_{dyn_p}	0.05
θ_vNa	1.65	s
θ_cNa	2.4	s
v_s⁽⁰⁾	0.8
f_r	0.25	Hz
v_p⁽⁰⁾	0 to 0.2 (default = 0.1)
k^r_p	0 to 0.3 (default = 0.1)

Table 2.

Settings for Simulating PSV. Sessions F1 to F5 Correspond to Increasing Respiratory Influence. Sessions G1 to G5 Correspond to Increasing Parasympathetic Activities

Session	v_p⁽⁰⁾	k^r_p	Description
F1	0.1	0.10	Investigate effect of k^r_p (increasing respiratory influence/magnitude)
F2	0.1	0.15
F3	0.1	0.20
F4	0.1	0.25
F5	0.1	0.30
G1	0.1	0.10	Investigate effect of v_p⁽⁰⁾ (Increasing parasympathetic activity)
G2	0.12	0.10
G3	0.14	0.10
G4	0.16	0.10
G5	0.18	0.10

Figure 9.

Simulated signals. Top to bottom: v_p(t), BP, T, respiratory activity, and PSV under condition F2 and e_{s_f_stat} = 0. Delay between vagal onset and pupillary constriction indicates parasympathetic dominance during expiration

Figure 10.

Simulated signals. Top to bottom: v_p(t), BP, T, respiratory activity, and PSV under condition F4 and e_{s_f_stat} = 0

Figure 11.

Simulated signals. Top to bottom: v_p(t), BP, T, respiratory activity, and PSV under condition G2 and e_{s_f_stat} = 0

Figure 12.

Simulated signals. Top to bottom: v_p(t), BP, T, respiratory activity, and PSV under condition G4 and e_{s_f_stat} = 0

Across Sessions F1 to F5, the higher range of k_rp produced a peak-to-peak PSV amplitude that increased linearly with the strength of the respiratory influence. This aligns with Ohtsuka’s finding that, once above a certain threshold, pupil fluctuation amplitude was proportional to electrical phrenic discharge in humans.⁹ In Sessions G1 to G5, stronger parasympathetic activity resulted in smaller mean pupil sizes, reduced HR, and lower BP, consistent with the well-established physiological actions of the parasympathetic nervous system (PSNS).

To further examine how ambient brightness conditions influence PSV, each simulation session was repeated with e_{s_f_stat} varied from −0.4 to 0.4 in increments of 0.2 to represent increasing ambient brightness. Figures 13 and 14 show the resulting PSV waveforms. Even with ANS-related parameters held constant, higher ambient brightness produced smaller mean pupil sizes in every session. The PSV peak-to-peak amplitude was smallest when the pupil was either highly constricted or dilated, and largest when the pupil size was in the mid-range. This finding provides a new interpretation of the RNL property first described by Stark et al.^7,8 They attributed the phenomenon to random Gaussian noise in the reflex arc, whereas our model suggests that PSV originates from ANS activities and that its statistical characteristics are shaped by the biomechanical properties of the pupillary muscle plant.

Figure 13.

Simulated PSV’s in Sessions F1 to F5: Increasing respiratory influence (k^r_p = 0.1, 0.15, 0.2, 0.25, 0.3). Top to bottom: PSV under increasing ambient brightness (e_{s_f_stat} = -0.4, -0.2, 0, 0.2, 0.4). The plots demonstrate the range-nonlinearity (RNL) effect, showing that PSV amplitude is biomechanically constrained at extreme pupil sizes

Figure 14.

Simulated PSV’s in Sessions G1 to G5: Increasing parasympathetic activity (v_p⁽⁰⁾ = 0.1, 0.12, 0.14, 0.16, 0.18). Top to bottom: PSV under increasing ambient brightness (e_{s_f_stat} = -0.4, -0.2, 0, 0.2, 0.4)

Consistent with Usui’s observations on pupillary fluctuations, the standard deviations of the simulated PSV signals were highest at mid-range pupil sizes (7-9 mm) and lowest at the lower and upper extremes (2 mm or 12 mm). However, unlike Usui’s findings, the resulting distributions were not Gaussian. This discrepancy likely arises because, rather than assuming Gaussian noise, our PSV signals were driven by ANS activity containing respiratory and baroregulatory components. Figure 15 shows the standard deviations of the simulated PSV plotted against mean pupil size, with the bottom-to-top traces corresponding to Sessions F1 through F5. These simulated behaviors establish the physiological validity of the model and set the stage for analyzing PSV and HRV in frequency domain.

Figure 15.

Plot of PSV standard deviation vs. PSV mean for F1 (bottom trace) to F5 (top trace)

Frequency Analysis of Simulated HRV and PSV

Analysis of power spectral density (PSD) of HRV is commonly performed to provide indices of autonomic balance, with high-frequency power (HF: 0.15-0.5 Hz) generally associated with parasympathetic activity and low-frequency power (LF: 0.04-0.15 Hz) reflecting both parasympathetic and sympathetic activities.⁶ To investigate whether the same interpretations can be applied to PSV, each simulated PSV signal in our study was low-pass filtered (FIR, f_c = 0.5 Hz, p = 200) and resampled at 1.25Hz before parametric spectral estimation using autoregressive (AR) modeling with the Burg method (p = 12). HRV signals were generated by constructing event series of delta functions at simulated heartbeat times, low-pass filtering them (FIR, f_c = 0.5 Hz, p = 200) for spectrum-of-counts analysis, and resampling at 1.25 Hz prior to Burg spectral estimation.

From the resulting PSV and HRV spectra, total power (TP: 0.04-0.5 Hz), LF power (LFP: 0.04-0.15 Hz), and HF power (HFP: 0.15–0.5 Hz) were computed as the areas under their respective frequency bands, with normalized indices defined as nLFP = LFP/TP and nHFP = HFP/TP. These indices were derived under the conditions corresponding to the different autonomic states in simulation Sessions F1–F5 and G1–G5, and under the various ambient brightness conditions. The aim is to explore, through simulation, whether the frequency characteristics of PSV and their associated indices can be interpreted in the same way as HRV, thereby supporting the use of PSV as a potential marker of autonomic function.

Figure 16 shows the frequency spectra of the simulated HRV signals under the five autonomic conditions of increasing respiratory influence (F1 to F5). A clear respiratory component at 0.25Hz was present in each HRV spectrum, with its peak amplitude increasing as the respiratory influence, k_rp, was raised from 0.1 to 0.3. A 0.1Hz component corresponding to the Mayer wave was also evident, although it was weaker and showed decreasing amplitude as k_rp increased. The corresponding HRV indices nLFP, nHFP, TP, and LF/HF for F1 to F5 are shown in Figure 17. As k_rp increased, nLFP decreased, nHFP increased, TP increased, and LF/HF decreased, indicating enhanced sympathetic and parasympathetic activities, and a shift in autonomic balance toward greater parasympathetic dominance.

Figure 16.

Frequency spectra of simulated HRV signal under five different autonomic conditions. Top to bottom traces: increasing respiratory influence (F1 to F5: k^r_p = 0.1, 0.15, 0.2, 0.25, 0.3)

Figure 17.

Top to bottom: nLFP, nHFP, TP, and LF/HF for simulated HRV under five different autonomic conditions of increasing respiratory influence (F1 to F5: k^r_p = 0.1, 0.15, 0.2, 0.25, 0.3). These trends reflect a shift in autonomic balance toward greater parasympathetic dominance as respiratory influence increases

Figure 18 shows the frequency spectra of the simulated HRV signals under the five conditions of increasing parasympathetic activity (G1 to G5). Respiratory and Mayer wave components were also evident at 0.25Hz and 0.1Hz, respectively. As parasympathetic activity, v_p⁽⁰⁾, increased from 0.1 to 0.18, the 0.1Hz component weakened while the 0.25 Hz component became more pronounced. The corresponding HRV indices nLFP, nHFP, TP, and LF/HF for G1 to G5 are shown in Figure 19. As v_p⁽⁰⁾ increased, nLFP decreased, nHFP increased, TP decreased, and LF/HF decreased. These changes indicate reductions in both sympathetic and parasympathetic activities, together with a shift in autonomic balance toward greater parasympathetic dominance. This transition from G1 to G5 resembles the trend typically observed during the recovery of parasympathetic activity.

Figure 18.

Frequency spectra of simulated HRV signal under five different autonomic conditions. Top to bottom traces: increasing parasympathetic activity (G1 to G5: v_p⁽⁰⁾ = 0.1, 0.12, 0.14, 0.16, 0.18)

Figure 19.

Top to bottom: nLFP, nHFP, TP, and LF/HF for simulated HRV under five different autonomic conditions of increasing parasympathetic activity (G1 to G5: v_p⁽⁰⁾ = 0.1, 0.12, 0.14, 0.16, 0.18)

Figure 20 shows the frequency spectra of the simulated PSV signals for autonomic condition F1, with v_p⁽⁰⁾= 0.1 and k_rp = 0.1, under five ambient brightness levels in which e_{s_f_stat} was varied from −0.4 to 0.4 in increments of 0.2. Similar to the HRV spectrum, the PSV spectrum exhibited a prominent respiratory component at 0.25Hz and a weaker Mayer wave component at 0.1Hz. Figure 21 presents the PSV spectra for the remaining autonomic conditions, F2 through F5, and different ambient brightness levels. As the respiratory influence increased across these sessions, the 0.25Hz respiratory component in the PSV became progressively stronger, whereas the 0.1Hz Mayer wave component weakened. This trend closely paralleled that observed in the corresponding HRV spectra.

Figure 20.

Top: frequency spectra of simulated PSV under autonomic condition F1 (v_p⁽⁰⁾ = 0.1, k_r^p = 0.1) and under five different ambient brightness conditions. Solid-thick (e_{s_f_stat} = -0.4), Dotted (e_{s_f_stat} = -0.2), Dashed (e_{s_f_stat} = 0), Dash-dot (e_{s_f_stat} = 0.2), Solid-thin (e_{s_f_stat} = 0.4). Bottom: Same frequency spectra magnified

Figure 21.

Frequency spectra of simulated PSV under different autonomic conditions and five ambient brightness conditions: solid-thick (e_{s_f_stat} = -0.4), dotted (e_{s_f_stat} = -0.2), dashed (e_{s_f_stat} = 0), dash-dot (e_{s_f_stat} = 0.2), solid-thin (e_{s_f_stat} = 0.4). Top-left: F2 (v_p⁽⁰⁾ = 0.1, k_r^p = 0.15), Top-right: F3 (v_p⁽⁰⁾= 0.1, k_r^p = 0.20), Bottom-left: F4 (v_p⁽⁰⁾ = 0.1, k_r^p = 0.25), Bottom-right: F5 (v_p⁽⁰⁾= 0.1, k_r^p = 0.30)

It is interesting to note that even when the autonomic condition remained unchanged, the shape of the PSV spectrum varied with the mean pupil size, which was altered by increasing ambient brightness from e_{s_f_stat} = −0.4 (darkest) to e_{s_f_stat} = 0.4 (brightest). The spectral power at 0.1Hz and 0.25Hz increased to a maximum as e_{s_f_stat} was raised from −0.4 to −0.2; however, as e_{s_f_stat} was further increased in steps of 0.2, both components diminished and reached their minimum at e_{s_f_stat} = 0.4. This pattern was likely a consequence of the pupil-size RNL effect: the biomechanics of the pupillary plant produce the largest fluctuations at mid-range pupil sizes and the smallest fluctuations when the pupil is at its extreme constriction or dilation. This RNL-driven modulation of the PSV spectra was observed consistently across all sessions from F1 to F5. The simulation results therefore suggest that the PSV spectral indices would be affected by both autonomic activity and ambient brightness.

Figure 22 further illustrates how nLFP, nHFP, TP, and LF/HF of the PSV changed under the different autonomic and ambient brightness conditions. All four PSV indices followed the same trends as their corresponding HRV indices (in Figure 17), but only when the ambient brightness condition (e_{s_f_stat}) was held constant. Under a fixed autonomic state, the indices varied systematically with e_{s_f_stat}. Specifically, nLFP, TP, and LF/HF increased to a maximum as e_{s_f_stat} was raised from −0.4 to −0.2, then decreased to a minimum as e_{s_f_stat} was further increased to 0.4, while nHFP exhibited the opposite pattern. These results suggest that the LF/HF ratio of PSV would be a reliable indicator of sympathovagal balance only when the mean pupil size remains unchanged. In the PSV model, mean pupil size is determined by both static autonomic tones and ambient brightness conditions.

Figure 22.

nLFP, nHFP, TP, and LF/HF for simulated PSV under five different autonomic conditions of increasing respiratory influence (F1 to F5: k^r_p = 0.1, 0.15, 0.2, 0.25, 0.3) and five different ambient brightness conditions. The five bars in each group correspond to e_{s_f_stat} = -0.4 (leftmost bar), -0.2, 0, 0.2, 0.4 (rightmost bar)

Figure 23 shows the frequency spectra of the simulated PSV signals for autonomic conditions G1 to G5 under the five ambient brightness levels. The PSV spectra contained a prominent respiratory component at 0.25 Hz and a weaker Mayer wave component at 0.1 Hz, similar to the HRV spectra. As parasympathetic activity, v_p⁽⁰⁾, increased from 0.1 to 0.18, the 0.25Hz respiratory component in the PSV strengthened, while the 0.1Hz Mayer wave component attenuated, consistent with the trend observed in the corresponding HRV spectra. The RNL-induced modulation of the PSV spectra was also evident across sessions G1 to G5: the spectral power at both 0.1Hz and 0.25Hz increased to a maximum as e_{s_f_stat} rose from −0.4 to −0.2, then decreased as e_{s_f_stat} was further increased, reaching a minimum at e_{s_f_stat} = 0.4.

Figure 23.

Frequency spectra of simulated PSV under different autonomic conditions and five ambient brightness conditions: Solid-thick (e_{s_f_stat} = -0.4), Dotted (e_{s_f_stat} = -0.2), Dashed (e_{s_f_stat} = 0), Dash-dot (e_{s_f_stat} = 0.2), Solid-thin (e_{s_f_stat} = 0.4). Top-left: G2 (v_p⁽⁰⁾ = 0.12, k_r^p = 0.10), Top-right: G3 (v_p⁽⁰⁾ = 0.14, k_r^p = 0.10), Bottom-left: G4 (v_p⁽⁰⁾ = 0.16, k_r^p = 0.10), Bottom-right: G5 (v_p⁽⁰⁾ = 0.18, k_r^p = 0.10)

Figure 24 shows that all four PSV spectral indices followed the same trends as their corresponding HRV indices (in Figure 19) when e_{s_f_stat} was held constant. Under a fixed autonomic condition, nLFP, TP, and LF/HF increased to a maximum as e_{s_f_stat} was raised from −0.4 to −0.2, then decreased to a minimum as e_{s_f_stat} was further increased to 0.4, while nHFP exhibited the opposite pattern. The results from sessions G1 to G5 therefore reinforce that the PSV LF/HF ratio can serve as a reliable indicator of sympathovagal balance only when the mean pupil size remains fixed.

Figure 24.

nLFP, nHFP, TP, and LF/HF for simulated PSV under five different autonomic conditions of increasing respiratory influence (G1 to G5: v_p⁽⁰⁾ = 0.1, 0.12, 0.14, 0.16, 0.18) and five different ambient brightness conditions. The five bars in each group correspond to e_{s_f_stat} = -0.4 (leftmost bar), -0.2, 0, 0.2, 0.4 (rightmost bar)

Frequency Response of the Pupillary Muscle Plant

With the autonomic condition fixed and ambient brightness varied, the patterns observed in nLFP, nHFP, TP, and LF/HF suggested that the RNL characteristics of the pupillary plant may introduce a frequency-dependent filtering effect on PSV. A plausible explanation is that when the pupil is near its minimum or maximum size, biomechanical constraints restrict PSV amplitude, resulting in relatively less attenuation of higher-frequency components than of lower-frequency components compared with mid-range pupil sizes. To investigate this filtering mechanism within our model, the input-dependent signals g_s(t) and g_d(t) were replaced with sinusoidal test inputs of amplitude 0.2 and frequencies spanning 0.04-0.25 Hz in 0.01Hz steps. To establish different mean pupil sizes, DC offsets were added to these test inputs: g_{s_stat} was varied from −0.4 to 0.4 in increments of 0.1, and g_{d_stat} was determined according to full reciprocal relationship. This procedure is analogous to determining the frequency response of the pupillary plant under different scenarios of ambient brightness conditions.

Figure 25 shows the frequency-response plots of the simulated PSV peak-to-peak amplitudes for g_{d_stat} values of −0.4, −0.3, −0.2, −0.1, 0, 0.1, 0.2, 0.3, and 0.4, corresponding to mean pupil sizes ranging from largest (darkest ambient condition) to smallest (brightest ambient condition). Although all traces displayed frequency-filtering characteristics, direct comparison was difficult because each began with a different magnitude at 0.04Hz. To enable clearer comparison, each frequency-response curve was normalized to its own magnitude at 0.04Hz and re-plotted, as shown in Figure 26.

Figure 25.

Peak-to-peak amplitudes of simulated PSV when g_s(t) and g_d(t) were sinusoidal testing signals of frequency ranging from 0.04Hz to 0.25Hz. Left (top to bottom): g_{d_stat} = -0.4 (solid), -0.3 (dashed), -0.2 (dotted), -0.1 (dash-dotted), 0 (thick). Right (top to bottom): g_{d_stat} = 0.1 (solid), 0.2 (dashed), 0.3 (dotted), 0.4 (dash-dotted)

Figure 26.

Peak-to-peak amplitudes of simulated PSV when g_s(t) and g_d(t) were sinusoidal testing signals of frequency ranging from 0.04Hz to 0.25Hz. All values are plotted relative to the peak-to-peak value at 0.04Hz. Left (top to bottom): g_{d_stat} = -0.4 (solid), -0.3 (dashed), - 0.2 (dotted), -0.1 (dash-dotted), 0 (thick). Left: g_{d_stat} = 0.1 (solid), 0.2 (dashed), 0.3 (dotted), 0.4 (dash-dotted)

In Figure 26, the frequency-response curves show that for pupil sizes ranging from largest to mid-range (corresponding to darkest to moderate ambient brightness), the PSV peak-to-peak amplitude decreases as frequency increases, demonstrating a clear low-pass filtering behavior. The slope of each curve depends strongly on the mean pupil size: as g_{d_stat} increases from −0.4 to −0.2, the rate of frequency roll-off becomes progressively steeper, with the steepest slope occurring at g_{d_stat} = −0.2. Beyond this point, the roll-off becomes more gradual as g_{d_stat} increases toward 0.1. From a filter design perspective, the pupillary plant behaves as a low-pass filter whose cutoff frequency is highest when the pupil is at its extremes and lowest when the pupil is in the mid-range. Notably, for the smallest pupil sizes (brightest ambient condition), this low-pass behavior largely disappears.

The abovementioned filtering characteristics of the pupillary muscle plant align with the ambient brightness-dependent patterns previously observed in nLFP, nHFP, TP, and LF/HF. Mid-range pupil sizes produced the highest nLFP and lowest nHFP because the lower cutoff frequency led to stronger attenuation of higher-frequency components. Conversely, pupil sizes at the extremes resulted in the lowest nLFP and highest nHFP, reflecting weaker attenuation of higher-frequency components due to the higher cutoff frequency.

Discussions

This paper presents the first computational model that explicitly couples cardiovascular autonomic regulation with a biomechanical pupillary plant to explain PSV, cardiovascular variability, and autonomic regulation. The model incorporates established neural pathways linking the cardiovascular system to a biomechanical pupillary muscle plant model. This overall framework successfully reproduced key empirical phenomena, and has clarified several unresolved issues regarding the origin and interpretation of PSV. The model accurately simulated the coupling between respiration, baroreflex dynamics, and pupillary dynamics. PSV amplitude increased during inspiration and decreased during expiration with a latency, consistent with findings from previous animal and human studies. The proportional increase in PSV amplitude with stronger respiratory influence reproduced the relationship between phrenic discharge and pupil oscillations reported previously. The model also reproduced the expected parasympathetic effects on pupil size and cardiovascular system. A key finding is the mechanistic explanation for the classical RNL phenomenon. While prior models attributed RNL to random Gaussian noise, our simulations across multiple ambient brightness conditions demonstrated that RNL emerges naturally from the nonlinear biomechanical properties of the iris muscles. The model reproduced the characteristic pattern where PSV amplitude is smallest at extreme pupil sizes and largest at mid-range. This behavior was observed consistently in both the time domain (PSV amplitude) and the frequency domain (LF/HF components), reinforcing the idea that RNL is a biomechanical property rather than a random artifact.

The frequency-domain results further revealed that PSV spectral indices, namely nLFP, nHFP, TP, and LF/HF, depend jointly on autonomic inputs and mean pupil size. When ambient brightness was held constant, PSV indices closely followed the trends of HRV indices under both respiratory and parasympathetic modulations. However, when pupil size changed due to ambient brightness levels, these indices shifted, even when the autonomic state was unchanged. This observation aligns with the RNL behavior and underscores the need to consider ambient brightness or mean pupil size when interpreting PSV-based autonomic metrics. The frequency-response experiments performed on the pupillary plant showed that the pupillary plant behaves as a low-pass filter with a cutoff frequency that depends on mean pupil size: attenuation of higher-frequency components is strongest at mid-range sizes and weakest at the extremes. This filtering profile matched the patterns observed in the PSV spectra and their indices, confirming that biomechanical filtering within the pupillary muscle plant is a primary determinant of the frequency characteristics of PSV.

Limitations

This study is entirely simulation-based. Our proposed model provides a theoretical framework for understanding the origins of PSV, rather than direct empirical validation. While the model successfully reproduces features comparable to established HRV and BPV dynamics, further experimental work involving concurrent in vivo measurements of PSV, HRV, and BPV is required to validate our findings against clinical data. Several limitations of our model should be acknowledged. First, there is lack of direct experimental validation using concurrent physiological recordings. Second, the model utilizes idealized neural pathways that simplify the highly complex subcortical interactions governing autonomic control. In real physiology, sympathetic drive and neural noise are significant contributors to PSV, particularly in the LF band. Our current model treats these sympathetic influences primarily as reciprocal to parasympathetic inputs; future iterations should incorporate independent sympathetic noise and drive to fully capture LF dynamics. Third, the simulations do not account for measurement noise or artifacts, such as blinks and saccades, which are prevalent in real-world situations. Finally, our current model focuses strictly on autonomic cardiovascular coupling and does not yet account for higher-level cortical influences such as cognitive or emotional load.

There is growing evidence that pupillary dynamics are highly influenced by cortical activity.²⁹ For instance, recent clinical hypotheses highlight the role of the locus coeruleus (LC) in higher-level neural control of pupil dynamics. Dysregulation of descending cortical inputs to the LC can alter sympathetic outflow and the inhibition of the EW nucleus, leading to abnormal pupillary responses in conditions such as ischemic brain injury.³⁰ Cortical lesions may also induce changes in pupillary reflex dynamics.³¹ These top-down are major factors in real-world PSV that extends beyond the cardiovascular-pupillary coupling modeled in this study.

PSV-Based Autonomic Monitoring

Designing a compensation technique for the low-pass filtering effect of the pupillary plant is technically challenging due to several physiological and computational complexities. The frequency response of the pupillary plant is highly nonlinear because of the iris muscles’ biomechanical properties, and its cutoff frequency varies continuously with mean pupil size, which depends on both ambient brightness and autonomic tone. As a result, the same physiological input may be attenuated differently under different lighting or autonomic conditions, making static linear compensation ineffective. The filtering effect is also asymmetric across frequencies, particularly between 0.1Hz and 0.25Hz, where Mayer waves and respiratory oscillations occur. Attenuation changes disproportionately across these bands, violating the assumptions needed for classical inverse filtering and causing inverse approaches to amplify noise or become unstable. In addition, any compensation method must handle biological variability, including inter-subject differences in iris stiffness, nonlinear sphincter–dilator interactions, and individual luminance adaptation. A robust solution therefore requires adaptive, model-based, or machine-learning-assisted methods that can estimate current pupil biomechanics and adjust the compensation dynamically.

Despite these challenges, our simulation results indicate that PSV-based autonomic monitoring is still feasible when measurements are taken under dark or moderately bright conditions, corresponding to larger or mid-range pupil sizes. In these regimes, the pupillary plant exhibits weaker or more predictable low-pass filtering, allowing respiratory and baroreflex-related components of PSV to remain largely intact. This suggests that reliable autonomic indices can be extracted from PSV in practice when monitoring durations are sufficiently long and lighting conditions are properly standardized.

Conclusion

This study introduced a novel, physiologically grounded computational model that integrates cardiovascular regulation, autonomic control, and pupillary biomechanics to simulate PSV and its interaction with HRV. The model reproduces major empirically observed phenomena, including respiration-linked oscillations, baroreflex effects, parasympathetic modulation, and the RNL characteristic, supporting the concept that PSV arises from autonomic inputs transmitted through the NTS–PBN–EW–CG pathway and shaped by intrinsic iris biomechanics. A key contribution of this work is providing mechanistic evidence that the pupillary plant behaves as a pupil-size-dependent low-pass filter whose cutoff frequency shifts with ambient brightness. This biomechanical filtering explains why PSV spectral indices depend not only on autonomic activity but also on mean pupil size. Consequently, PSV-derived metrics such as LF, HF, and LF/HF can be meaningful indicators of autonomic balance only when pupil-size effects are controlled. By providing a unified and physiologically consistent framework, this model lays important groundwork for using PSV as a biomarker in eye-based wearable systems.

Footnotes

ORCID iD

Kevin Hung

Author Contributions

Kevin Hung contributed to project management and funding resources, research design, model design, model implementation, data analysis, and manuscript writing. Gary M.-T. Man contributed to data analysis and manuscript writing. Daniel H.-K. Chow reviewed and edited the manuscript and provided valuable guidance for the interpretation and discussion of the findings.

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was fully supported by the Research Grants Council of the Hong Kong Special Administrative Region, China, under Grant UGC/FDS16/E15/21.

Declaration of Conflicting Interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

References

World Health Organization . Ageing and health. 2025. Available at: https://www.who.int/news-room/fact-sheets/detail/ageing-and-health. accessed 28 February 2026.

Townsend

Timmis

Erglis

, et al. Challenges in the prevention, treatment and management of cardiovascular disease among older adults. The Lancet Regional Health - Europe. 2025;56:101372.

World Health Organization . Cardiovascular diseases (CVDs). 2025. Available at: https://www.who.int/news-room/fact-sheets/detail/cardiovascular-diseases-(cvds. accessed 28 February 2026.

Villareal

Liu

Massumi

. Heart rate variability and cardiovascular mortality. Current Atherosclerosis Reports. 2002;4(2):120-127.

Stevens

Wood

Koshiaris

, et al. Blood pressure variability and cardiovascular disease: systematic review and meta-analysis. British Medical Journal. 2016;354:i4098.

Pagani

Lombardi

Guzzetti

, et al. Power spectral analysis of heart rate and arterial pressure variabilities as a marker of sympatho-vagal interaction in man and conscious dog. Circulation Research. 1986;59(2):178-193.

Stark

Campbell

Atwood

. Pupil unrest: an example of noise in a biological servomechanism. Nature. 1958;182(4639):857-858.

Stanten

Stark

. A statistical analysis of pupil noise. IEEE Transactions on Biomedical Engineering. 1966;13(3):140-152.

Ohtsuka

Asakura

Kawasaki

Sawa

. Respiratory fluctuations of the human pupil. Experimental Brain Research. 1988;71(1):215-217.

10.

Calcagnini

Giovannelli

Censi

Bartolini

Barbaro

(2001) Baroreceptor-sensitive fluctuations of heart rate and pupil diameter. In: 2001 Conference Proceedings of the 23rd Annual International Conference of the IEEE Engineering in Medicine and Biology Society; 25-28 October 2001; Istanbul, Turkey, Volume 1, pp. 600-603. IEEE.

11.

Koç

Baysal

Özbek

. Does a change in iris diameter indicate heart rate variability? Biomedical Research. 2018;29(7):1340-1344.

12.

Park

Kim

. Acupuncture has differential effects on HRV and pupil size in postures of healthy subjects. Autonomic Neuroscience. 2007;135(1):53.

13.

Day

. A Model of Accommodative-Pupillary Dynamics. Dissertation. Iowa, USA: Iowa State University; 1969.

14.

Dennison

. A Mathematical Model for the Motor Activity of the Cat Iris. Dissertation. Massachusetts, USA: Worcester Polytechnic Institute; 1967.

15.

Semmlow

Stark

. Simulation of a biomechanical model of the human pupil. Mathematical Biosciences. 1971;11:109-128.

16.

Hansmann

. Human Pupillary Mechanics: Physiology and Control. Dissertation. Berkeley, California, USA: University of California; 1972.

17.

Usui

Stark

. A model for nonlinear stochastic behavior of the pupil. Biological Cybernetics. 1982;45:13-21.

18.

Sun

Krenz

Stark

. A systems model for the pupil size effect. Biological Cybernetics. 1983;48:101-108.

19.

Denison

Parker

Carrasco

. Modeling pupil responses to rapid sequential events. Behavior Research Methods. 2007;52:1991-2007.

20.

Borgdorff

. Respiratory fluctuations in pupil size. American Physiological Society Journal. 1975;228(4):1094-1102.

21.

Censi

Calcagnini

De Pasquale

, et al. Baroreceptor-sensitive fluctuations of human pupil diameter. Computing in Cardiology. 1999;26:233-236.

22.

Chen

Eldridge

Wagner

. Respiratory-associated rhythmic firing of midbrain neurons in cats: relation to level of respiratory drive. The Journal of Physiology. 1991;437:305-325.

23.

Bianca

Komisaruk

. Pupil dilatation in response to vagal afferent electrical stimulation is mediated by inhibition of parasympathetic outflow in the rat. Brain Research. 2007;1177:29-36.

24.

Seidel

Herzel

. Modelling heart rate variability due to respiration and baroreflex. In: Mosekilde

Mouritsen

, eds. Modelling the Dynamics of Biological Systems: Nonlinear Phenomena and Pattern Formation. Berlin, Heidelberg: Springer; 1995:205-229.

25.

Bahler

. Modeling of mammalian skeletal muscle. IEEE Transactions on Biomedical Engineering BME-. 1968;15(4):249-257.

26.

Meiss

. Some mechanical properties of cat intestinal muscle. American Physiological Society Journal. 1971;220(6):2000-2007.

27.

Hirata

Usui

. A nonlinear dynamical model for human pupillary muscle plant. IEICE Transactions on Information and Systems J77-D-. 1994;2(1):170-180.

28.

Yana

Okuyama

Yoshida

, et al. An Evidence of Blood Pressure and Respiration Originated Variations in Fluctuations of Pupil Diameter. In: Proceedings of the Annual International Conference of the IEEE Engineering in Medicine and Biology Society; Orlando, Florida, USA; 31 October-3 November, 1991, pp. 472-474. New York: IEEE.

29.

Peinkhofer

Knudsen

Moretti

Kondziella

. Cortical modulation of pupillary function: systematic review. PeerJ. 2019;7:e6882.

30.

Scala

Miccoli

Pafundi

, et al. Automated pupillometry is able to discriminate patients with acute stroke from healthy subjects: an observational, cross-sectional study. Brain Sciences. 2024;14(6):616.

31.

Scala

Miccoli

Di Giovanni

, et al. Automated pupillometry is a predictor of outcome of stroke patients: an observational, prospective, cohort study. Brain Communications. 2025;7(1):fcaf079.