Ultradian rhythms in heart rate variability and distal body temperature anticipate onset of the luteinizing hormone surge

The menstrual cycle is characterized by predictable patterns of physiological change across timescales. Although patterns of reproductive hormones across the menstrual cycle, particularly ultradian rhythms, are well described, monitoring these measures repeatedly to predict the preovulatory luteinizing hormone (LH) surge is not practical. In the present study, we explored whether non-invasive measures coupled to the reproductive system: high frequency distal body temperature (DBT), sleeping heart rate (HR), sleeping heart rate variability (HRV), and sleep timing, could be used to anticipate the preovulatory LH surge in women. To test this possibility, we used signal processing to examine these measures in 45 premenopausal and 10 perimenopausal cycles alongside dates of supra-surge threshold LH and menstruation. Additionally, urinary estradiol and progesterone metabolites were measured daily surrounding the LH surge in 20 cycles. Wavelet analysis revealed a consistent pattern of DBT and HRV ultradian rhythm (2–5 h) power that uniquely enabled anticipation of the LH surge at least 2 days prior to its onset in 100% of individuals. Together, the present findings reveal fluctuations in distal body temperature and heart rate variability that consistently anticipate the LH surge, suggesting that automated ultradian rhythm monitoring may provide a novel and convenient method for non-invasive fertility assessment.

Premenopausal and perimenopausal estradiol, luteinizing hormone and progesterone metabolites. Participants monitored LH for all 55 cycles, whereas daily urine samples were collected by 20 women (n = 16 premenopausal, n = 4 perimenopausal) for the measurement of E2, α-Pregnanediol (αPg) and β-Pregnanediol (βPg). E2, αPg and βPg were collected to confirm that hormone concentrations were within healthy ranges for pre-menopausal women and that LH surges were followed by a rise in progesterone metabolites. For all 16 cycles, estradiol, αPg and βPg fell within normal ranges, with a pre-LH rise in E2 (2 days prior to LH onset through LH onset day were significantly greater than all other days, p < 0.01 in all cases). Likewise, LH surge onset was concomitant with a significant rise in αPg (p < 0.05 on LH onset, and < 0.01 6 days after LH onset and thereafter) ( Fig. 1A−C) and βPg (p < 0.005 on LH onset, and < 0.001 3 days after LH onset; data not shown for βPg). These hormonal changes were associated with a rise in temperature deviation above zero and a non-significant elevation of breathing rate around LH surge onset (Supplemental Fig. 1E-F). Consistent with previous findings 71 , LH surge length was variable, with 42% of individuals exhibiting supra-threshold LH concentrations 2 days following surge onset, falling to 26% of individuals 3 days after surge onset (Fig. 1A). LH was tonically supra-threshold in perimenopausal women (n = 10 cycles, Fig. 1D www.nature.com/scientificreports/ exhibited a significant increase in αPg and βPg only 6 days after midcycle (p < 0.05; data not shown for βPg), and a trend toward elevation of E2 prior to mid cycle (p = 0.176) (Fig. 1D−F).
Ultradian power of DBT, HRV, and LH surge onset. Ultradian (2−5 h) power of daytime DBT exhibited a stereotyped pattern preceding LH surge onset in premenopausal ( Fig. 2A,B), but not perimenopausal (Fig. 2C,D), cycles. Ultradian DBT power exhibited an inflection point a mean of 5.82 (± 1.82) days prior to LH surge onset and a subsequent peak a mean of 2.58 (± 1.89) days prior to the surge onset. A second trough in UR power occurred a mean of 2.06 (± 1.02) days after surge onset (χ 2 = 5.66, p = 0.0174,). These stereotyped changes were not present in perimenopausal cycles (χ 2 = 0.37, p = 0.5354, for the same comparisons) (Fig. 2C). Inflection point and subsequent peak of DBT and HRV ultradian power anticipate LH surge onset. In premenopausal women, the first inflection point of DBT and HRV (RMSSD) UR power occurred between -8 and -2 days prior to surge onset, whereas the subsequent peak in UR power for both metrics occurred between -6 days before to 2 days after LH surge onset (Fig. 4). 85% of cycles exhibited the first inflection point by 4 days prior to the surge, with 100% showing this inflection by 2 days prior to the surge. The peak of UR power occurred at least 1 day prior to the surge in 82% of cycles. Together, these inflection points and subsequent peaks in UR power of HRV (RMSSD) and DBT uniquely anticipated the LH surge days before its onset (see "Discussion" for potential relevance to fertile window). www.nature.com/scientificreports/

Discussion
The present findings reveal stereotyped fluctuations in DBT and HRV (RMSSD) UR power that anticipate 100% of LH surge onsets, a key component of female health and fertility. By contrast, changes in DBT circadian rhythm power were not predictive of the LH surge, suggesting that URs are uniquely coupled to the pre-ovulatory time of the menstrual cycle. Likewise, discrete, nightly behavioral and physiological measures did not anticipate the surge, suggesting that continuous measures of physiological output provide signals more amenable to LH surge anticipation. Finally, these features did not occur stereotypically in perimenopausal cycles with respect to mid cycle. These findings point to peripheral URs as oscillations that are coupled to menstrual cycle physiology and that have the potential to contribute to the development of tools for estimating the female fertile window. Although the underlying physiological mechanisms that lead to systematic changes in UR power in DBT require further investigation, much is known about general changes in body temperature across the menstrual cycle. Estrogens lower, and progestins raise, body temperature 53,72 . Accordingly, body temperature reaches a minimum, with minimum core circadian amplitude, during the late follicular phase and rises in the core, mouth, and skin following ovulation 73 . Body temperature also broadly reflects metabolic rate, which is elevated in the late follicular and luteal phases 74 . In mice, the structure of core temperature URs allows for the detection of female reproductive state, with a high plateau of temperature and trough of UR power during the active phase indicative of the LH surge and ovulation 13,75 . Most human studies to date have focused on core temperature, measured via an ingestible device that travels through the GI tract 73 , intravaginal or rectal sensor 76 , or oral thermometer 60 . However, ultradian, circadian, and ovulatory rhythms in temperature are readily observed at the periphery, providing several advantages: (1) DBT has higher amplitude fluctuations than core body temperature, making URs and CRs easier to detect 77 , (2) changes in DBT may correlate with sleep stage 78 , and (3) DBT is in circadian antiphase to core temperature, but shows the same general trend across the menstrual cycle, suggesting comparable reliability 77 . It is possible that rising UR power of DBT before the LH surge reflects higher UR power of reproductive hormones during this time. www.nature.com/scientificreports/ Whereas body temperature is the most commonly used non-hormonal output in menstrual cycle tracking, previous studies have found that HRV also changes by cycle phase and may therefore be a candidate for surge anticipation 79 . Parasympathetic input to the heart dominates during the follicular phase, lowering resting heart rate and elevating HRV (RMSSD) 57 . Sympathetic input to the heart dominates during the luteal phase, elevating heart rate and depressing HRV (RMSSD) 11,57 . Consequently, HRV (RMSSD) varies ~ 10 ms from the follicular to the luteal phase 11 , with a marked decrease in the latter portion of the cycle 52 . These fluctuations may be more difficult to detect during a short daytime recording window, and are impacted by daytime activities, making sleep an ideal window over which to look for unmasked features 45 . Natural negative controls illustrating reproductive and metabolic influences on HRV patterns are that (1) LH pulsatility is disrupted in obese and diabetic women 80,81 , and (2) mid cycle and luteal fluctuations in HRV are absent in polycystic ovarian syndrome (PCOS), a leading cause of female infertility 58,59 . In the present study, sleeping HRV (RMSSD) UR power rose in the late follicular phase, peaked near the LH surge, and dropped sharply before rising into the early luteal phase. Although the present study lacks sufficient power to evaluate other potential patterns that may be relevant to the menopausal transition, the preliminary absence of comparable features in perimenopausal individuals suggests that this group deserves further study. Together, signal processing of DBT and HRV could yield actionable information for individuals and clinicians wishing to estimate the "fertile window". However, there are several challenges inherent to accurately defining the female fertile window.
The fertile window (the time during which a woman may become pregnant) depends upon many factors, including (1) the timing of the LH surge, (2) the subsequent time of the release of the ovum or ovulation, (3) the presence of a viable corpus luteum releasing adequate progesterone 38 , (4) the duration of time sperm can survive in the female body, which is dependent both on sufficient number and quality of sperm and on the appropriate vaginal environment (e.g., pH) 22,82 , and 5) quality of the uterine environment. Most investigations report the highest probability of fertility as the 5 days preceding ultrasound-determined day of ovulation (USDO) 83 , but actual days on which an individual may become pregnant are much more variable, with pregnancy occurring up to 11 days prior to ovulation to 5 days after ovulation 71 .
Some of the reported variability in the fertile window likely results from discrepancies in language used to describe both human ovulation and the fertile window itself 84 . Despite their namesake, home "ovulation tests" that identify supra-threshold LH concentrations do not measure ovulation, which may occur many days after and occasionally a few days before LH surge onset 71 . Despite this variability, the fertile window is often treated as predictable, with definitions including the 5-6 days prior to the LH surge as a proxy for ovulation 85 , the first www.nature.com/scientificreports/ day of slippery clear cervical fluid through LH surge onset 86 , the total days of slippery clear cervical fluid 87 , day 10-17 of the cycle 88 , and retrospective measures of salivary ferning 85 , basal body temperature 64,89 , and progesterone metabolites (e.g., 90 ). Today, many online and app-based ovulation prediction algorithms are validated using day of cycle or LH data alone, in the absence of hormone measures or USDO 2,83,84,91 . Additionally, extant data sets regularly report excluding 20-50% of collected data due to cycle irregularities, without determining if given cycles were hormonally aberrant 71,[92][93][94] . Together, the confounding of the LH surge with ovulation and the variable criteria used to define the fertile window make it difficult to accurately determine the variance, and contributors to variability, of fertility relative to the LH surge or ovulation. Despite these discrepancies, the possibility that UR features anticipate the onset of the LH surge by a few to several days suggests applicability for family planning. When one considers the additional time between LH surge onset and ovulation, these features may anticipate much to all of the fertile window. If confirmed in larger cohorts, this method would constitute the earliest method of predicting a definitive event at any point within the fertile window. Open source, non-invasive methods for predicting the LH surge as a marker of likely future ovulation are not currently available 71 , but the present findings indicate that the onset of the LH surge may be anticipated days in advance by automated detection of changes in ultradian power of DBT and HRV (RMSSD). These changes consistently anticipate LH surge onset in women of a variety of ages, cycle lengths, surge timing and duration. The frequency band of 2-5 h examined in the present investigation was not specifically selected for the present group of participants but chosen based on the peak frequency band observed across physiological systems 14,25,56 , suggesting potentially broad applicability. Due to the high demand for accurate methods of fertility assessment, such novel methods carry the responsibility to clearly report the aspects of reproductive physiology that are detected and the methods by which detection is achieved once algorithms are tested on large populations 68,84,91 .
Future work will determine the extent to which these ultradian rhythm-based methods of menstrual cycle monitoring generate accurate predictions in larger, more diverse cohorts. In particular, the study of a greater number of cycles within individuals may enable personalization of relevant features. With these data, methods such as empirical mode decomposition for selection of tailored ultradian bands, or machine learning based methods for assessment across many different features at once, may result in greater specificity or longer predictive windows. These features could potentially be used on their own, with minimal user input (e.g., tracking of dates of menstruation), or in combination with other FAM methods. Ideally, such methods could be widely employed on wearable devices such as the Oura Ring, or on future generations of convenient and precise body temperature and HRV sensors. Together, these findings may guide further research aimed at understanding how hormones, metabolism, and the autonomic nervous system temporally interact; and may aid the development of open-source, non-invasive methods of fertility awareness.

Ethical approval. This study and all procedures were approved by the Office for the Protection of Human
Subjects at the University of California, Berkeley. All participants gave informed consent. All research was performed in accordance with relevant guidelines and regulations.
Participants and recruitment. Participants were recruited from the Quantified Self community, a global group of individuals interested in learning through self-measurement 95,96 . Individuals attended a prospective discussion about the project at the 2018 Quantified Self meeting in Portland, Oregon and contacted the experimenter via email if interested in participating. Prospective participants were contacted to discuss study structure, risks and benefits, and to review the informed consent form. Once informed consent was obtained, participants were instructed to complete an introductory questionnaire with their age, cycling status (regular, irregular, recovering from hormone/IUD use, perimenopausal, menopausal), and historical LH surge day(s), if known. Contact information was collected for the purposes of communication and delivery of study materials. Data from pregnancies (n = 3) that overlapped with the study were excluded from these analyses. Participants had not taken hormonal contraception within the prior year and did not have any known reproductive medical concerns. There were no age or parity restrictions, consistent with the principles of participatory research 96,97 . See Table 1 for participant demographics.

Study design.
Each of the 28 (n = 20 premenopausal, n = 5 perimenopausal, n = 3 premenopausal and became pregnant) participants collected 2 to 3 cycles of data for analysis. For all cycles, the Oura Ring, a DBT, HRV (RMSSD), HR, and sleep sensor, was worn continuously on the finger, as previously described 65,98 . For all cycles, LH was monitored via urinary test strips (Wondfo Biotech Co., Guangzhou, China) from day 10 (with first day of menstruation considered day 1) until a positive reading was detected, and subsequently until 2 days after LH fell below the limit of detection (see below for details on the Urinary Hormone Assay, Luteinizing Hormone). Of the 55 total cycles collected (45 = premenopausal, 10 = perimenopausal), 20 were paired with daily, morning urine tests for E2, αPg and βPg, the major urinary progesterone metabolites (Precision Analytical, McMinnville, OR). This study was designed using the principles of participatory research 96,99 in which individual participants maintain control of their own data prior to anonymization and came to the project with personal questions that could be answered with the data to be collected. All participants received a copy of their Oura Ring data.
Data collection and management. HR, HRV (RMSSD), DBT, sleep onset, sleep offset, sleep duration, breathing rate, and nightly temperature deviation (described briefly below) were collected using the Oura Ring (Oura Inc., San Francisco, CA; Oura Health Oy, Ltd., Oulu, Finland). The Oura Ring is a small, wireless sensor worn on the finger. By using an LED light source and LED sensor to measure reflection off the skin above www.nature.com/scientificreports/ the radial artery of the finger, the Oura Ring calculates HR, HRV (RMSSD), and breathing rate. The ring also contains 3 thermistors for detection of DBT. DBT is measured 24 h a day (binned in 1-min intervals). To avoid artifacts associated with activity, HR, and HRV (RMSSD) are only measured during sleep (binned in 5-min intervals), limiting our analyses of HR and HRV (RMSSD) to the sleeping period. All other metrics are calculated once per night. Briefly, the body temperature deviation for each night is the moving mean of nightly temperature between 10:00 pm and 8:00 am, minus the mean temperature of the previous 20 days. Oura Rings were loaned to the group by Oura Inc. The Oura Ring can be connected to a mobile phone application, Oura, via Bluetooth. At the start of the study, each participant downloaded the Oura application from either the Google Play Store (Google Inc., Mountain View, CA) or the Apple App Store (Apple Inc, Cupertino, CA) to their mobile phones and created an Oura account. Participants were able to view their own data provided by the application throughout the study. Participants were asked to synchronize data from the ring to the application each morning. Uploaded data was automatically transferred via the internet to the study database in the Oura cloud service. In order to access data from the cloud, data were imported into the Open Humans 100 framework, which provides encrypted, password protected data access to researchers, with the participants' revocable consent. In addition to data collected by the Oura Ring, participants uploaded personal spreadsheets that tracked days of menstruation, days of LH tests and results, days of urine collection, and notes (e.g., forgot to wear the Oura Ring) to Open Humans. Participants could opt out of the study and remove their data at any time. Data were anonymized by the researchers for analysis. Data, once anonymized at the end of the study, remained in the data set.
Hormone assays. For the assessment of E2, αPg and βPg, participants collected daily, first-morning urine samples across a cycle according to manufacturer's instructions (Precision Analytical, Willamette, OR). Briefly, a standardized piece of filter paper with an attached label was submerged in the urine sample and dried for 24 h. Filter paper was then frozen at ~ -18 C in participants' home freezers until analysis. E2, αPg, and βPg were analyzed using proprietary in-house assays referred to as Dried Urine Testing for Comprehensive Hormones (DUTCH) on the Agilent 7890/7000B GC-MS/MS (Agilent Technologies, Santa Clara, CA, USA). The equivalent of approximately 600 μl of urine was extracted from the filter paper using acetate buffer. In the first week of the cycle, and from 3 days after LH surge completion until the end of the luteal phase, samples were pooled every 2 days (a third day was pooled at the end of cycles in instances where the total number of remaining days after the surge was odd). Urine samples were extracted and analyzed as previously described, with previously established ranges of hormone concentrations expected in urine by phase of cycle and during menopause 90,101 . Briefly, creatinine was measured in duplicate using a conventional colorimetric (Jaffe) assay. Conjugated hormones were extracted (C18 solid phase extraction), hydrolyzed by Helix pomatia and derivatized prior to injection (GC-MS/ MS) and analysis. The mean inter-assay coefficients of variation were 8% for E2, 12% for αPg, and 13% for βPg. The mean intra-assay coefficients of variation were 7% for E2, 12% for αPg and 12% for βPg. Sensitivities of the assays were as follows: E2 and αPg, 0.2 ng/mL; βPg, 10 ng/ mL.
Luteinizing hormone was measured using the commercially available WondFo (Wondfo Biotech Co., Guangzhou, China) Luteinizing Hormone Urinary Test 102 , a validated at-home urine assay. Briefly, the strip was submerged by participants for 5 s in a fresh urine sample and laid horizontally for 5 min before reading. When samples were collected for E2, αPg and βPg, those same samples were used for LH testing. Each strip contains a positive control and a "test" line, indicating if LH is present in the urine at, or over a concentration of 25 MIU/ mL 102 . Test results were depicted as either + or -(no quantitative information provided) and were recorded in a personal spreadsheet by the participant. A photograph of each test was taken by participants to ensure accurate reading of the results.
Inclusion and exclusion criteria for collected cycles. Cycles were included in the premenopausal data set as likely ovulatory by four criteria a) one or more localized days of supra-threshold LH concentration, b) the presence of a rise in E2 (if collected) within typical range prior to or coincident with supra threshold LH, c) a subsequent rise in αPg and βPg (if collected), and d) positive values of DBT deviation, as previously described 98 , within 2 days of surge onset until the end of the cycle (See Supplemental Fig. 1). Cycles without E2, αPg, and βPg data were included by meeting criteria a and d only. Cycles with missing data within sixteen days of the of the LH surge (defined as no HR/HRV/DBT data for a given night) were omitted in order to avoid erroneous estimation of rhythmic power (see Data Analysis below). Cycles were defined as "perimenopausal" by the presence of positive LH measured at least every other day across the cycle and age > 45 years. Four such cycles were paired with daily urinary hormone analysis for E2, αPg, and βPg, as described above.
Data analysis. All code and data used in this paper are available at Open Science Framework (https ://osf. io/wzf47 /). Code was written in Matlab 2019b, Matlab 2020a and Python 3. Wavelet Transform (WT) code was modified from the Jlab toolbox and from Dr. Tanya Leise 103 . Briefly, data were imported from the Open Humans framework to Python 3, where HR, HRV (RMSSD), and DBT data were extracted. Data were cleaned in Matlab, with any data points outside + /-4 standard deviations set to the median value of the prior hour, and any points showing near instantaneous change, as defined by a local abs(derivative) > 10 5 as an arbitrary cutoff, also set to the median value of the previous hour.
Wavelet transformation (WT) was used to assess the structure of ultradian rhythms of DBT, HR, HRV (RMSSD), and circadian rhythms in DBT. As DBT shows high plateaus during the sleeping period, and URs during the day, DBT analyses here were used on data collected during the waking hours (see Supplemental  Fig. 2). Conversely, as indicated previously, because the Oura Ring only collects HR and HRV (RMSSD) during sleep, wavelet analyses were restricted to the sleeping window for these metrics. In either case, the excerpted www.nature.com/scientificreports/ data were compiled from all days of the cycle resulting in one continuous signal representing all days (DBT) or all nights (HR, HRV (RMSSD)). In contrast to Fourier transforms that transform a signal into frequency space without temporal position (i.e., using sine wave components with infinite length), wavelets are constructed with amplitude diminishing to 0 in both directions from center. This property permits frequency strength calculation at a given position. Wavelets can assume many functions (e.g., Mexican hat, square wave, Morse); the present analyses use a Morse wavelet with a low number of oscillations (defined by β and γ), analogous to wavelets used in previous circadian applications 103 . Morse Wavelet parameters of β = 5 and γ = 3 describe the frequencies of the two waves superimposed to create the wavelet; Additional values of β (3)(4)(5)(6)(7)(8) and γ (2-5) did not alter the findings (data not shown) 104 . This low number of oscillations enhances detection of contrast and transitions. The band of the wavelet matrix corresponding to 2-5 h rhythms were averaged in order to create a linear representation of UR WT power over time. This band corresponded to the timescale of ultradian rhythmicity observed across physiological systems 14,25,56 . Potential changes to circadian power of DBT (mean power per minute within the 23-25 h band) were additionally assessed prior to extracting days for ultradian-only analyses, but no significant changes across the cycle were detected (see Supplemental Fig. 3). Because WTs exhibit artifacts at the edges of the data being transformed, only the WT of the second through the second to last days of data were analyzed further.
To enable comparisons across cycles of different durations, premenopausal cycles were displayed from LH surge onset minus 7 days to LH onset plus 7 days. As perimenopausal individuals had tonically high LH, and therefore no surge onset to which all individuals could be aligned, each cycle's midpoint was chosen for alignment.
LH surge anticipation features. Wavelet power in the 2-5 h band was calculated as described above.
Extracted bands were smoothed using a daily moving average using the Matlab function "movmean". The Matlab function "findpeaks" was used to identify peaks as points at which either adjacent point had a lower UR power. This function was run on the negative of the signal to identify troughs. Points at which the derivative of the signal crossed zero, indicating a change in direction of UR power (i.e., either increasing to decreasing or vice versa), were found using the matlab function "diff ". The first time the derivative crossed zero in the cycle (i.e., the first inflection point), excluding the first five days of the cycle, during which LH is very unlikely to rise, was marked as the presence of the first feature for either HRV (RMSSD), DBT, or HR. Following this inflection point, the next peak identified by "findpeaks" was marked as the second feature. These methods of identifying peaks, troughs, and direction changes were used to ensure the diff function was identifying all visually identified peaks.

Statistical analyses.
Descriptive values are reported as means ± daily standard deviations (SD) unless otherwise stated. For statistical comparisons of average ultradian power in premenopausal and perimenopausal cycles, Kruskal Wallis (KW) tests were used instead of ANOVAS to avoid assumptions of normality for any distribution to assess the trend in average UR power leading up to the surge as compared to after the surge. For KW tests, χ 2 and p values are listed in the text. One-way repeated measures analysis of variance (rmANOVA) tests were used to compare peak average E2 to other days surrounding the surge, and baseline αPg and βPg (7 days prior to the surge) to other days surrounding the surge. For rmANOVAs, p values are listed in the text. Because the dominant trend was an inflection point in UR power followed by a peak, slopes of individual signals were compared rather than raw values at each timepoint. The same tests were applied to individuals, in addition to tests for significance of raw power differences on peak and trough days, using 25 min centered on peaks and troughs, respectively. To avoid multiple comparisons and chance of a type I error, differences between individually-determined peak and trough values of UR WT power found using "findpeaks" were assessed using a KW test, such that each cycle contributed only 1 peak value and 1 trough value (N = 45 data points per group).