Heart rate and age modulate retinal pulsatile patterns

Theoretical models of retinal hemodynamics showed the modulation of retinal pulsatile patterns (RPPs) by heart rate (HR), yet in-vivo validation and scientific merit of this biological process is lacking. Such evidence is critical for result interpretation, study design, and (patho-)physiological modeling of human biology spanning applications in various medical specialties. In retinal hemodynamic video-recordings, we characterize the morphology of RPPs and assess the impact of modulation by HR or other variables. Principal component analysis isolated two RPPs, i.e., spontaneous venous pulsation (SVP) and optic cup pulsation (OCP). Heart rate modulated SVP and OCP morphology (pFDR < 0.05); age modulated SVP morphology (pFDR < 0.05). In addition, age and HR demonstrated the effect on between-group differences. This knowledge greatly affects future study designs, analyses of between-group differences in RPPs, and biophysical models investigating relationships between RPPs, intracranial, intraocular pressures, and cardiovascular physiology.

I n vivo dynamic video-ophthalmoscopy (VO) provides a potential opportunity for a non-invasive and easily accessible evaluation of retinal hemodynamics. VO is inexpensive and well suited to become a diagnostic and disease monitoring tool for the real-time imaging of local microvascular blood supply and detecting various pathologies. VO is applicable in many fields of medicine such as ophthalmology (e.g., diabetic retinopathy 1 , glaucoma [2][3][4], neurology (e.g., Alzheimer disease 5 , multiple sclerosis 6 , stroke 7 ), or cardiology (e.g., coronary heart disease 8 , arterial stiffness [9][10][11] , hypertension 12 , diabetes 12,13 ).
Blood flow 14,15 , blood volume 4,16 and structural venous diameter changes [17][18][19] are the most commonly evaluated hemodynamic parameters from dynamic retinal imaging. The hemodynamic parameters are usually extracted in a manual or semi-automated fashion from specific morphological segments of a retinal vessel tree [16][17][18][20][21][22][23] . However, the reproducibility of the parameters remains a challenge due to the uniqueness of each individual's retinal vessel tree 24 and the bias introduced by a subjective inter-rater and inter-participant selection of analyzed morphological segments 25 . To increase the reproducibility of local retinal hemodynamics, we have recently implemented a blind source separation that automatically divided the optic nerve head (ONH) into 2-3 functionally distinct areas that emerged as specific retinal pulsatile patterns (RPPs) 26 . The reproducible RPPs were spontaneous venous pulsation (SVP) and optic cup pulsation (OCP) with mutually phase-shifted hemodynamics. OCP was postulated to represent arterial blood filling preceding the SVP outflow 26 .
Spontaneous venous pulsation is the most investigated hemodynamic phenomenon over the whole retinal background with the best detection capability and highest reproducibility in the area where the central retinal vein crosses the ONH. The etiology and visibility of the ONH SVP are likely related to the gradient between intracranial pressure (ICP) and intraocular pressure (IOP) waveforms [27][28][29][30] . Limited in vivo evidence of the SVP-ICP relationship exists [30][31][32][33][34] . The underlying biophysics is a subject of active investigation 35 .
Mathematical models, which were proposed to describe IOP-ICP-SVP relationships [35][36][37][38] , have not been verified by in vivo experiments yet. Levine and Bebie's SVP-ICP theory assumes an influence of heart rate (HR) on SVP amplitude, yet without providing the in vivo evidence 37,38 . Therefore, the influence of HR on retinal intra-vessel hemodynamics deserves further investigation. In addition to HR, demographic characteristics (e.g., age or sex) were linked to vessel stiffness 39,40 , vessel crosssection area 41,42 , pulse pressure 39,40 or cardiac cycle parameters (e.g., filling time, preload, stroke volume) and, consequently, modify intra-vessel hemodynamics including the blood flow/ volume in retinal mini-and micro-vessels 42 . Exactly how the in vivo video-ophthalmoscopic data (sensitive to blood flow or volume fluctuations) are affected by HR, age, or sex remains to be determined. If the impact of any variable is proved, a revisiting of clinical video-ophthalmoscopic outcomes demonstrating correlations between VO measurements and retinal neural fiber layer (RNFL) thickness 4,19,43 or other outcomes in populations of healthy subjects and those with disease conditions such as glaucoma is warranted.
In this study, we investigated human in vivo videoophthalmoscopic data and the effects of HR and age on the morphology of blood-volume-specific RPPs. We isolated two RPPs (i.e., SVP and OCP) from the video-ophthalmoscopic dataset, tested SVP-OCP phase shift, evaluated RPP reproducibility and morphology, and cross-correlated the morphology with HR, age, IOP, and RNFL thickness. Finally, we estimated the effects of HR and age with between-group comparisons in resulting morphological observations.

Results
Participant characteristics. Thirty-four retinal video-recordings (RVRs) were acquired, and exclusively left-eye RVRs were used in the analysis. HR estimated from SVP and OCP of all participants was 66 ± 13 min −1 (ranging 44-92 min −1 ). (Image analysis workflow estimating and segregating SVP and OCP is summarized in Fig. 1.) HR was significantly higher in patients with treated ocular hypertension (OHT) than in healthy participants and was significantly correlated with SVP and OCP morphologies ( Table 1, Figs. 2, 3). Due to this finding, HR was treated as a confounding variable in further ANCOVA between-group tests. Other physiological data such as age, IOP, and RNFL demonstrated no significant between-group differences ( Table 1). Refractive error, visual acuity, and perimetry were within physiological ranges without any significant pathology.
Phase shifted retinal pulsatile patterns in principal component space. In a total of 34 RVRs, we detected SVP in 33 RVRs (97%) and OCP in 31 RVRs (91%) when the first 12 principal components were visually inspected and evaluated in each RVR. Averaged SVP principal component index was 3 ± 2 (min 1; max 9), and averaged OCP principal component index was 3 ± 2 (min 1; max 11). Representative examples of SVP and OCP spatiotemporal patterns with detected and corrected control points are shown in Fig. 1a. An input-output workflow of the utilized principal component analysis (PCA) is schematically summarized in Fig. 1b.
Between-group differences in RPP morphology measurements. Lower SVP Ampl, SVP V T , and OCP t p , along with higher absolute OCP SlpD and OCP SlpU, were detected in treated OHT patients (two-sample t-test) ( Table 1). ANCOVA with age as a confounding variable preserved all significant between-group differences in OCP morphology, but none in SVP morphology (Table 1). In contrast, ANCOVA with HR as a confounding variable only showed trends of the higher absolute OCP SlpD and the lower OCP t p (Table 1). Group averaged OCP pulses with 25-75% confidence intervals are shown in Fig. 6 and demonstrate visible changes in group-specific pulse morphologies.

Discussion
Study novelty and practical impact. Our in vivo results demonstrate the modulation of the SVP and OCP pulse morphologies by HR and of the SVP morphology by age. Therefore, the unverified biological impact of HR and age at RPP morphology was in vivo confirmed and validated. The impact of HR and age on the SVP/OCP morphology has been neglected in previous clinical research 4,19,33,34,43 and even in theoretical SVPintracranial/intraocular pressure models 35,36 . The in vivo findings validate Levine and Bebie's theory, assuming that SVP Ampl is influenced by HR 37,38 . This discovery may also indicate descriptive vessel compliance characteristics as demonstrated by the uniformity of the SVP SlpD parameter. The studied morphology parameters should be investigated in future research as objective quantitative in vivo markers of retinal hemodynamics. Importantly, HR and age need to be considered and added as confounding variables in biological models or clinical studies evaluating between-group differences in RPP. The optimal strategy for the study design is a dataset free of between-group differences in HR and age, but such proper matching would be challenging and hardly achievable.
The analysis via spatial PCA of RVRs detected reproducible phase-shifted SVP and OCP patterns in healthy controls and treated OHT patients. SVP and OCP patterns demonstrated low inter-participant variation in eigenvector scaled time-courses, i.e., temporal pulses, and were highly reproducible. These results emphasize the feasibility of the spatial PCA for inferences achieved by RPPs and its applicability in future studies, especially those involving large population cohorts.
Although the etiology of the SVP phenomenon is still not clearly understood, the temporal uniformity in SVP SlpD can represent a novel insight into the SVP origin. This temporal uniformity may be due to similar inter-participant venous compliance characteristics that modulates venous capacitance. As the implication of vein capability to change its geometry, the vein compliance property enables a volume increase in local vein segments in response to local blood filling 45 . Since the vein resistance to local blood flow is minimal 45 , the blood filling corresponding to SlpD in the SVP morphology displayed temporal uniformity over investigated population.
We revealed the significant impact of HR on the morphology of reproducible RPPs. The higher HR resulted in lower SVP Ampl, V T , and t p and higher SlpU. Similarly, the higher HR led to the lower OCP V T , t p , and SlpD and higher SlpU. This observation of changes in V T linked to HR may reflect the impact of Starling's law. Higher HR causes a shorter period for cardiac blood filling and, consequently, lower cardiac preload leads to the lower blood volume ejected from the heart. All these characteristics of heart function imprint into the RPP morphology and may be evaluated non-invasively and in vivo with the VO.
Additionally, three SVP parameters, i.e., Ampl, SlpD, and V T , were significantly correlated with age. Higher SVP Ampl and V T with age can be influenced by peripheral pulse pressure pulsatility or vessel stiffening, which increase with age 11,46 . The increase of age-related vessel stiffness is directly linked to the elasticity loss 39,40 and is a potential factor underlying the steeper SVP SlpD trend observed in our older participants. Mechanistically, lesser venous capacitance leads to faster venous volume outflow. If the hypothesis is true, the SVP SlpD may become a non-invasive and fully automated measure proportional to vessel stiffness.
Previous studies showed that SVP were influenced by the IOP 28,35,37,47 and RNFL significantly correlated with RPP morphological measurements in the dataset, including population with retinal neurodegenerative disesase 4,19,43 . In our participants, the IOP and RNFL did not correlate with RPP morphology. A potential explanation may lie in well-controlled OHT condition in our patient's cohort who presented with normalized IOP, so no neurodegeneration assessed by RNFL was present.
SVP morphology is believed to originate from a gradient between IOP and ICP, and many theoretical models attempted to express the SVP-IOP-ICP relationships [35][36][37][38] . However, the majority of previously reported models have not accounted for the effects of HR or age on the SVP Ampl or timings [35][36][37][38] . Therefore, the current in vivo observation represents critical information that adds another piece to the puzzle of the SVP etiology conundrum. Therefore, future models and studies evaluating SVP should account for the effects of HR and age to avoid potential erroneous conclusions.
Similar conclusions apply to the OCP assessment. Although OCP morphology demonstrated significant differences between healthy controls and OHT patients, a considerable deal of uncertainty remains as the OCP morphology changes were HRrelated. The ANCOVA demonstrated the HR effect on betweengroup comparisons. It is important to note that the impact of HR has usually been omitted in previous clinical studies. In particular, Although the PCA is an automated method, the presented data processing pipeline utilized three manual interventions with the potential impact on outcome measures. The interventions were: (i) manual ONH segmentation; (ii) identification of SVP and OCP components; and (iii) control point corrections. As PCA integrates eigen time-course for each principal component estimated from the whole ONH area, minimal imperfections in ONH segmentation should not crucially impact outcome measures. But, a further study investigating the PCA method robustness is warranted to clarify this matter. The automated Fig. 2 Heart rate modulated morphology of spontaneous venous pulsations (SVP). a Visualization of mean single-pulses over all monocular (mono) and binocular (bino) retinal video-recordings from healthy controls and patients with medicated ocular hypertension (OHT). Graph lines are heart rate color-coded. b Evaluation of linear dependence between heart rate and SVP morphology measurements (Ampl amplitude of the SVP eigenvector, V T total relative pulse stroke volume in the eigenvector measures, SlpD slope-down from the eigenvector value at the period beginning to the negative eigenvector peak ≈ peak of the maximal absolute blood volume time-point, SlpU slope-up from the negative eigenvector peak to the period end, t p time to the negative eigenvector peak). Value r represents a corresponding Pearson correlation coefficient and value p the p-value of the correlation level. Fig. 3 Heart rate modulated morphology of optic cup pulsations (OCP). a Visualization of mean single-pulses over all monocular (mono) and binocular (bino) retinal video-recordings from healthy controls and patients with medicated ocular hypertension (OHT). Graph lines are heart rate color-coded. b Evaluation of linear dependence between heart rate and OCP morphology measurements (Ampl amplitude of the OCP eigenvector, V T total relative pulse stroke volume in the eigenvector measures, SlpD slopedown from the eigenvector value at the period beginning to the negative eigenvector peak ≈ peak of the maximal absolute blood volume time-point, SlpU slope-up from the negative eigenvector peak to the period end, t p time to the negative eigenvector peak). Value r represents a corresponding Pearson correlation coefficient and value p the p-value of the correlation level. Our study has several limitations to be addressed in future research. The sample size of our dataset is limited and should be extended to reproduce and validate our pilot observations. Still, the utilized correction for multiple comparison errors and the power analysis support that the presented correlation measures are related to the human body physiology rather than false positive observations. An optimal re-test should involve hundreds of participants with RVRs acquired at various videoophthalmoscopes of various assembling technologies. Therefore, we would like to initialize a multi-center RVR challenge collecting 10 s RVRs to re-test the impact of HR and age at RPP morphology. Four participants who had both monocular and binocular RVR, approximately two years apart, may decrease inter-participant variability in our dataset. We consider this effect of being minimal as HR differed between the first and second RVR in three of four participants and as the intra-participant variability mostly followed dataset trends of correlation measurements. Pharmacological treatment of OHT patients may have influenced the morphology measurements, as specific drugs can alter ocular hemodynamics. The high diversity of used drugs in our OHT group and our small sample size prevented the testing post-hoc differences between the untreated group and treated subgroups (e.g., beta-blocker versus prostaglandin treatment). The differences in vessel stiffness and pulse pressure related to sex 49,50 may also play a significant role in RPP morphology. Sex effects need to be validated in the larger cohort of participants.
In conclusion, we have demonstrated the in vivo evidence that heart rate and age modulate the morphology of retinal pulsatile patterns in humans. The observation corroborates Levine and Bebie's theory. The presented study will impact the design of future biological and clinical studies, future analyses of betweengroup differences in morphology of RPPs, and SVP-ICP-IOP biophysical modeling by emphasizing the necessity to include heart rate and age as important confounding factors.

Methods
Experimental design. In concordance with the Declaration of Helsinki, all participants signed an informed consent approved by the ethical committee at the Friedrich-Alexander University of Erlangen-Nürnberg. The participants are part of Fig. 4 Representative examples of heart rate (HR) dependent intra-participant variability of retinal pulsatile patterns in one healthy control and two medicated ocular hypertension (OHT) patients. For each participant, one retinal video-recording (RVR) was acquired with the monocular videoophthalmoscope (VO) utilizing the CCD camera chip and one with the binocular VO utilizing the CMOS camera chip. The between-acquisition time interval was about two years for each representative participant. In healthy control and one OHT participant (OHT1), intra-participant RPP morphologies were dissimilar at different HR while RPP morphology remained unchanged for comparable HR as captured for other OHT participant (OHT2). Another OHT participant had different HR over acquisitions with similar outcomes as presented for the OHT1 participant. SVP spontaneous venous pulsations, OCP optic cup pulsations, HR heart rate.
Participants enrolled in the NCT00494923 trial met the criteria of age range 18-65, open chamber angle and corrected visual acuity 0.7 or better when entering the registry. The trial excluded people with systemic disease and potential ocular involvement (e.g., diabetes mellitus), people with myopic or hyperopic refractive error >8.0D, and people with an eye disease (except for glaucoma). From the registry, healthy controls and patients with a history of ocular hypertension (OHT; i.e., intraocular pressure IOP > 21 mmHg, normal visual field and ONH appearance) were included in this study. Pharmacological treatment along with eye surgical intervention history are described in Table 2.
Monocular or binocular retinal video-recordings (RVRs) were acquired for 16 OHT patients (age 58.7 ± 12.9 years old, seven females) and 14 healthy controls (age 66.0 ± 13.2 years old, eight females) between January 2015 and December 2017. In three OHT patients and one healthy control, two RVRs (one monocular and one binocular) were obtained approximately two years apart. In total, 34 RVRs were acquired and analyzed.
Along with RVRs, refractive error, visual acuity, perimetry, IOP, and retinal nerve fiber layer (RNFL) thickness were obtained for each eye with standard clinical devices (white-on-white perimetry with computerized static projection perimeter, Octopus 500, Haag-Streit; Goldmann tonometry; Spectralis OCT, Heidelberg Engineering) followed by RVR acquisition. All measurements were acquired within a one-day session.
Video-ophthalmoscopic data acquisition. Each participant was examined while comfortably sitting with head rested and positioned on a video-ophthalmoscope chin holder to minimize participant's motion. Each participant was asked to fixate the eyesight at the target presented as a red LED light or cross in the videoophthalmoscope optical path. The principles of image acquisition were previously described 16 . In short, both available video-ophthalmoscope types (monocular and binocular) acquire images of the reflected light intensity modulated by heart rate induced attenuation changes 16,26 . Because such changes are caused by spatialtemporal retinal blood volume changes, the lowest pixel image intensity corresponds to the highest blood volume (and the highest attenuation) and vice versa.
Image analysis. The whole image analysis workflow is summarized in Fig. 1a. Motion artifacts in RVRs were suppressed with rigid image registration optimized for the RVRs 53 . A representative example of aligned monocular RVR is available at https://youtu.be/-CABIpjWX8Y and binocular RVR at https://youtu.be/ 4anapI0TZTQ. Next, the optic nerve head (ONH) area was manually segmented from the averaged aligned RVR image and defined as a region of interest (ROI) for further data analysis. Each ONH time frame was spatially smoothed with a 3 × 3 uniform convolution kernel to increase the signal-to-noise ratio (SNR) between local retinal hemodynamics and additive Gaussian noise in RVRs. Acquired relative blood volume changes were high-pass filtered in spectral domain with the cut-off frequency 0.12 Hz in each aligned pixel belonging to the ONH. The filter suppressed DC component and low-frequency drift but preserved pixel-specific pulsatile variance.
Spatial principal component analysis 26 (PCA) was estimated via singular value decomposition (Eq. 1) on each preprocessed ONH RVR.
Formation of PCA input (i.e., ONH RVR matrix X of m rows and n t columns where m is a number of analyzed pixels and n t is a number of time points) and extraction of PCA outputs from left and right eigenvector matrices U and V are SVP spontaneous venous pulsations, OCP optic cup pulsations, Ampl amplitude of the single-pulse in eigenvector space, V T total relative pulse stroke volume in the eigenvector measures, SlpD slope-down from the eigenvector value at the period beginning to the negative eigenvector peak ≈ peak of the maximal absolute blood volume time-point, SlpU slope-up from the negative eigenvector peak to the period end, t p time to the negative eigenvector peak. Z-scored spatial principal components (Fig. 1b) were thresholded to zero in each pixel where |Z | < 1 26 . The SVP 26 and OCP 26 spatiotemporal patterns were visually identified as a single component for each pattern from a set of the first 12 principal components. PCA eigenvectors characterizing pulsation time-courses (Fig. 1b) were de-trended, and outlier values were restored utilizing the k-means clustering algorithm, as both implemented and fully described in Labounkova et al. 26 . Control points defining continuous part of an RPP time-course with high SNR (see OCP control points in Fig. 1a) were automatically identified for each SVP or OCP eigenvector 44 (Fig. 1b) characterizing the relative blood volume changes 26 , and minor manual edits were done if needed. The automated identification of the control points is detailly described in the "Appendix C" of Labounkova et al. 26 .
Each RVR segregated several SVP or OCP single pulse repetitions whose beginning and end were well defined by the control points (Fig. 1a). Averaged single SVP or OCP pulse waveforms were derived for each RVR, and quantitative parameters describing its morphology were evaluated. The evaluated morphology parameters were pulse amplitude (Ampl; Fig. 1a), total relative pulse stroke volume is a volume in non-positive eigenvector values, T is pulse period, and t is time; Fig. 1a), slope-down to peak (SlpD; Fig. 1a), slope-up to baseline (SlpU; Fig. 1a), and time to peak (t p ; Fig. 1a). The averaged RVR HR = 60/T[min −1 ] was estimated from averaged SVP and OCP periods T[s].
The Ampl measurement is proportional to the maximum blood volume in the examined ROI during the cardiac cycle. V T is proportional to the quantity of total blood volume change in the ROI during the cardiac cycle. SlpD is proportional to the steepness of the blood volume filling in the ROI during the cycle and SlpU to the steepness of blood volume drainage outside the ROI. t p is proportional to the time of blood volume filling in the ROI during the cycle.
Statistics and reproducibility. Delay between overlapping high-quality SVP and OCP portions (Fig. 1a) was evaluated with maximal cross-correlation response function 55 . Pearson correlation coefficients between unaligned or aligned SVP and OCP patterns were quantified (Fig. 1a). The null hypothesis that SVP-OCP delay equals to 0 was tested with one sample t-test.
Pearson correlation coefficients and corresponding correlation p-values were evaluated between SVP or OCP morphological features (i.e., Ampl, V T , SlpD, SlpU, and t p ; Fig. 1a) and age, HR, IOP, and average RNFL thickness, respectively. Ten correlation effects were investigated for each variable (e.g., age, HR). False discovery rate correction (p FDR <0.05) was applied to the correlation p-values, and Benjamini-Hochberg adjusted p-values (p BH ) were computed to minimize the risk of the error type I. Due to the limited sample size, power was estimated for each significant correlation to assess the risk of error type II 56 . Two-sample t-test evaluated between-group differences in RPP morphology measurements, age, HR, IOP, and RNFL. ANCOVA (analysis of covariance) was the second statistical test evaluating between-group differences when HR or age were used as confounding variables. Results of two-sample t-tests and ANCOVA tests were compared, and the effects of HR or age on final between-group results were evaluated (Fig. 1a). Due to the limited sample size, we considered uncorrected p < 0.05 significant for t-tests and ANCOVA.
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
De-identified data sets can be made available upon a reasonable email request to Ivana Labounkova (ilabounk@umn.edu), Dr. Folkert Horn (folkert.horn@augen.imed.unierlangen.de), or another responsible personnel from the Department of Ophthalmology and University Eye Hospital, Friedrich-Alexander University Erlangen-Nürnberg at Erlangen, Erlangen, Germany.

Code availability
The MATLAB R2017b programming environment (MathWorks, Natick, USA) with academic license and the open-source Retina Imaging Toolbox (https://github.com/umnmilab/retinaimagingtoolbox; GNU GPL version 3 license) were used for all image and statistical analyses and visualizations. The table describes the active substances used to treat high intraocular pressure with the classification of the substance.