HormoneBayes: A novel Bayesian framework for the analysis of pulsatile hormone dynamics

The hypothalamus is the central regulator of reproductive hormone secretion. Pulsatile secretion of gonadotropin releasing hormone (GnRH) is fundamental to physiological stimulation of the pituitary gland to release luteinizing hormone (LH) and follicle stimulating hormone (FSH). Furthermore, GnRH pulsatility is altered in common reproductive disorders such as polycystic ovary syndrome (PCOS) and hypothalamic amenorrhea (HA). LH is measured routinely in clinical practice using an automated chemiluminescent immunoassay method and is the gold standard surrogate marker of GnRH. LH can be measured at frequent intervals (e.g., 10 minutely) to assess GnRH/LH pulsatility. However, this is rarely done in clinical practice because it is resource intensive, and there is no open-access, graphical interface software for computational analysis of the LH data available to clinicians. Here we present hormoneBayes, a novel open-access Bayesian framework that can be easily applied to reliably analyze serial LH measurements to assess LH pulsatility. The framework utilizes parsimonious models to simulate hypothalamic signals that drive LH dynamics, together with state-of-the-art (sequential) Monte-Carlo methods to infer key parameters and latent hypothalamic dynamics. We show that this method provides estimates for key pulse parameters including inter-pulse interval, secretion and clearance rates and identifies LH pulses in line with the widely used deconvolution method. We show that these parameters can distinguish LH pulsatility in different clinical contexts including in reproductive health and disease in men and women (e.g., healthy men, healthy women before and after menopause, women with HA or PCOS). A further advantage of hormoneBayes is that our mathematical approach provides a quantified estimation of uncertainty. Our framework will complement methods enabling real-time in-vivo hormone monitoring and therefore has the potential to assist translation of personalized, data-driven, clinical care of patients presenting with conditions of reproductive hormone dysfunction.

present hormoneBayes, a novel open-access Bayesian framework that can be easily applied to reliably analyze serial LH measurements to assess LH pulsatility.The framework utilizes parsimonious models to simulate hypothalamic signals that drive LH dynamics, together with state-of-the-art (sequential) Monte-Carlo methods to infer key parameters and latent hypothalamic dynamics.We show that this method provides estimates for key pulse parameters including inter-pulse interval, secretion and clearance rates and identifies LH pulses in line with the current gold-standard deconvolution method.We show that these parameters can distinguish LH pulsatility in different clinical contexts including in reproductive health and disease in men and women (e.g., healthy men, healthy women before and after menopause, women with HA or PCOS).A further advantage of hormoneBayes is that our mathematical approach provides a quantified estimation of uncertainty.Our framework will complement methods enabling real-time in-vivo hormone monitoring and therefore has the potential to assist translation of personalized, data-driven, clinical care of patients presenting with conditions of reproductive hormone dysfunction.

Introduction
Pulsatile hormone dynamics are ubiquitous and play a crucial role in the regulation of many bodily functions related to metabolism, stress, and fertility 1,2 .Hormones are typically secreted in both a basal manner to maintain steady state levels, as well as with superimposed interspersed transient bursts (pulses) 3 .It is now established that the pulsatile nature of hormonal secretion affects their interaction with receptors and downstream effector action [4][5][6] .
With regards to fertility, the hypothalamus is the central regulator of the reproductive endocrine axis.Notably, gonadotropin releasing hormone (GnRH) is secreted in a pulsatile manner, and seminal studies have demonstrated that this pulsatility is fundamental for its action to stimulate GnRH receptors on pituitary gonadotropes 5 .Moreover, disturbances in GnRH pulsatility are observed in common reproductive disorders including polycystic ovary syndrome (PCOS) in which GnRH pulsatility is increased 7 , and hypothalamic amenorrhea (HA) in which GnRH pulsatility is reduced 8 .However, despite this disparate alteration in GnRH pulsatility, differentiation of these two common reproductive disorders, which may both present similarly with menstrual disturbance, can be challenging 9 .LH is measured routinely in clinical practice using an automated chemiluminescent immunoassay method and is the gold standard surrogate marker of GnRH.Furthermore, LH can be measured at frequent intervals (eg 10minutely) to assess GnRH/LH pulsatility, and accurate assessment of LH pulsatility could help facilitate diagnosis and treatment of patients presenting with reproductive endocrine disorders 9 .However, this is rarely done in clinical practice because it is resource-intensive, inconvenient for patients, and there is a lack of user-friendly methods for computational analysis of the LH data available to clinicians.
Analysis of hormone pulsatility is a challenging computational problem, primarily due the stochastic nature of hormone dynamics and the consequent pulse-to-pulse variability, but also due to extrinsic factors (such as measurement error) obscuring the observed hormone dynamics 3 .Several computational methods for the analysis of endocrine data have been proposed in the literature 3,[10][11][12][13] , and deconvolution analysis, is the current gold-standard method for analyzing LH pulsatility in humans 3 .However, all methods lack an open-access, user-friendly interface that can be readily used by clinicians.
To meet these challenges, we have developed HormoneBayes, a novel, openaccess Bayesian framework for the analysis of hormone pulsatility data.Unlike previous methods, HormoneBayes uses a stochastic model, describing hormone levels in the circulation incorporating measurement error, and leverages Bayesian statistics 14 to infer model parameters and latent variable dynamics.We show that HormoneBayes can be used to accurately identify LH pulses and estimates clinically relevant measures such as inter-pulse intervals and secretion rates.The framework also provides a handle on estimation uncertainty via Bayesian posterior distributions.
We showcase how this feature can be used to enable the understanding of alterations in LH pulsatility by analyzing the effect of direct hypothalamic stimulation using the neuropeptide kisspeptin on a subject-by-subject basis.We also demonstrate that HormoneBayes can be used to analyze LH pulsatility in different clinical contexts/reproductive states (including healthy men, women before and after menopause, and women with reproductive disorders such as HA or PCOS).
Importantly, the framework comes with a user-friendly, open-access graphical interface that make the core functionality of the framework easily accessible to clinicians in the clinic.

Analyzing pulsatile hormone dynamics using the hormoneBayes framework
The hormoneBayes framework allows inference of key physiological parameters describing pulsatile hormone dynamics.The framework utilizes stochastic mathematical models describing circulating hormone levels and state-of-the-art Bayesian machinery to calibrate these models to data of hormone profiles and infer model parameters.Figure 1 presents a parsimonious model (a simple model with great explanatory power) describing circulating LH levels.The model assumes that there are two modes of LH secretion, namely pulsatile and basal.The former corresponds to an on/off signal that randomly switches between two states (corresponding to a high and a low activity) while the latter corresponds to a continuous noisy signal.Furthermore, the model incorporates LH clearance as a linear first order process leading to an exponential decay of LH levels following a pulse 3 .We note that more detailed clearance models, such as the bi-exponential model 3 , could be easily integrated.By considering the processes of LH release and clearance, the model predicts LH circulating dynamics in terms of five key parameters that can be recovered from data: 1. LH clearance rate; 2. maximum LH release rate; 3. strength of the pulsatile signal relative to the basal signal; 4. mean time in the on state; 5. mean time in the off state.Moreover, the model incorporates measurement error as an additional parameter that is determined based on the assay coefficient of variation (CV).
HormoneBayes relies on the Bayesian paradigm to extract information from the observed data.Using the Bayes theorem, the method revises our prior beliefs regarding model parameters by transforming the parameters' prior probability density distributions into posterior distributions.Parameter prior distributions enable the user to input context-specific information into the analysis, hence enhancing the flexibility of the method to handle different datasets.For example, when dealing with data from post-menopausal women the user could choose to adjust the parameter priors to acknowledge a higher LH secretion rate and/or more sustained basal secretion.

Identifying LH pulses
Using data to infer the latent hypothalamic signal provides a transparent way to identify pulses based on their likelihood under the model.As explained above, the model assumes LH pulses are partly driven by an on/off hypothalamic signal.This latent variable (i.e., inferred variable) takes two values indicating the 'on' (1) and 'off' (0) state of the hypothalamic pulse generator, and therefore the expected posterior estimate (inferred from LH profiling data) can be interpreted as the probability that at any given time the hypothalamic pulse generator is 'on'.This quantitative measure for accessing the likelihood of a pulse can significantly ease the job a clinician trying to decide whether an upstroke in the LH profile represents a pulse or not. Figure 2 illustrates a representative example of an LH trace with two obvious pulses occurring at times 100min and 300min.These are indeed identified by inspecting the hypothalamic signal profile which peaks at around those times.At time 200min a less pronounced bump in LH could be indicative of an LH pulse, however the inferred pulsatile hypothalamic signal remains well below 0.5, indicating that under the current model there is higher probability that the bump is a measurement artefact rather than a pulse.
To validate our pulse identification method, we used a database of LH profiles obtained from healthy pre-menopausal women and compared the identified pulses with those obtained by the current gold-standard deconvolution method, which uses a maximum likelihood approach to fit a secretion model to the data 3 .Here, we identify a pulse when the posterior probability that the hypothalamic signal is 'on' exceeds a threshold value.We use 0.5 as the threshold value, which signifies there is higher probability that the hypothalamic pulse generator is on (rather than off).This value yields the highest agreement between hormoneBayes and the deconvolution method (see fig S2; SI).We find that hormoneBayes agrees with the deconvolution method in 77 out of 87 identified pulses (89%).Moreover, 13 pulses identified by hormoneBayes were not identified by the deconvolution method.
Overall, this suggests a good agreement between the two methods, with hormoneBayes having the added advantage of providing a measure for the likelihood of each pulse that clinicians and researchers can use to inform their clinical decision making.

Variation of model parameter within and across groups
To test the applicability of hormoneBayes in different contexts we compile a database of LH profiles from four groups with diverse LH dynamics (men, healthy pre-menopausal women with regular menstrual cycles, post-menopausal women, women with HA, and women with PCOS).As illustrated in Fig. 3, the model successfully captures the differences in LH dynamics in all four groups.Moreover, the model allows us to summarize LH data through model parameters and assess the variability across and within groups.We find that two model parameters explain most of the variability between groups, namely the maximum secretion rate and pulsatility strength (Fig 3C).The first parameter describes how much LH can be secreted over time, whilst the second parameter quantifies the strength of the pulsatile signal relative to the basal signal.Based on these two parameters there is a strong distinction between women with PCOS and women with HA, who have lower LH secretion rates and/or diminished LH pulsatility strength (i.e., pulsatile signal is weak relative to the basal signal).Furthermore, post-menopausal women display higher secretion rates compared to pre-menopausal women but also reduced pulsatility strength (i.e., pulses are less pronounced when higher level of LH are established post menopause).Interestingly, healthy men and women illustrate a much lower parameter variability as compared to HA and post-menopausal women, which could be indicative of various degrees of severity of HA/PCOS and tighter LH regulation in healthy individuals of reproductive age.

Discussion
We have presented HormoneBayes, a novel computational framework for analyzing hormone pulsatility.The framework combines (i) mathematical models describing hormone dynamics with (ii) computational Bayesian machinery for inferring model HormoneBayes leverages the Bayesian paradigm to infer model parameters from the data.The Bayesian approach means that hormoneBayes outputs a (posterior) density distribution for model parameters that has an intuitive interpretation, i.e, it describes the probability that parameters attain some value given the data.This enables the use of posterior distributions for statistical testing, for example to assess the effect of hormonal interventions on LH secretion parameters at the population level but also on a subject-by-subject basis.We expect such personalized analysis to become mainstream as measurement technologies mature enabling cheap sampling of hormone levels in real time 15 .
Ultimately, data analysis using HormoneBayes is as credible as the underlying model used to describe hormonal dynamics.Here, we have used a parsimonious model that assumes two modes of LH secretion, pulsatile and basal.This assumption is in par with current physiological understanding of the system and the hypothalamic pulse generator hypothesis 16,17 .Furthermore, to model LH circulation levels the model assumes a linear clearance rate.At least one other model used for the analysis of LH pulsatility has used more complex (multiple timescale) clearance dynamics, however we expect this assumption should have a minimal impact for the purpose of assessing LH pulsatility.Nevertheless, the modular design of HormoneBayes allows future extensions of the model with the scope of comparing how well different models capture LH dynamics as well as enabling the analysis of hormone dynamics beyond LH.

Stochastic model of LH
We used a discrete-time, stochastic model to describe pulsatile LH dynamics.The model comprises of three dynamical variables, ܲ ௧ , ‫ܤ‬ ௧ , and ‫ܪܮ‬ ௧ , that describe the pulsatile and basal hypothalamic signals and the LH concertation in circulation, respectively.
The pulsatile hypothalamic signal ܲ ௧ can take two values: ‫ܪ‬ ௧ ൌ 0 corresponding to the ON state; and ‫ܪ‬ ௧ ൌ 1 corresponding to the OFF state.The stochastic dynamics of ‫ܪ‬ ௧ are governed by the following probability matrix i.e., parameters ߬ ைே and ߬ ைிி govern the probabilities that the value of ‫ܪ‬ will either flip or remain the same over the time interval ሺ‫,ݏ‬ ‫ݏ‬ ߜ‫ݐ‬ሻ.
The evolution of the basal hypothalamic signal, ‫ܤ‬ ௧ , is described using a discrete time autoregressive model obeying the following equation where ߝ ௧ is a normally distributed random variable with zero mean and unit variance.
The two hypothalamic signals drive LH secretion, and along with LH clearance dictate the circulating LH levels, ‫ܪܮ‬ ௧ .The equation describing the time evolution of ‫ܪܮ‬ ௧ , is where ݀ denotes the clearance rate, ݇ denotes the maximum secretion rate, and parameter ݂ (termed pulsatility strength) describes the relative strength of the two hypothalamic signals.Finally, the model assumes measurement error in the form: where ߟ ௧ is a normally distributed random variable with zero mean and std.deviation equal to the CV of the assay.Throughout our analysis we have used ‫ݐߜ‬ ൌ 1 min, hence, assuming that the system dynamics do not change significantly over shorter times.

Figure 1
Figure 1 illustrates how the Bayesian machinery allows us to calibrate the model to parameters from data.Using a parsimonious mathematical model of LH secretion, we have demonstrated the clinical utility of HormoneBayes in accurately analyzing LH pulsatility data, identifying pulses, and summarizing data in terms of model parameters to predict different clinical conditions/reproductive states on an individual patient level basis.
Sampling from the full posterior distribution is performed using a Gibbs sampler, which is an iterative Monte Carlo Markov Chain (MCMC) scheme.The algorithm is initialised with parameter values drawn from the prior distribution, i.e., Θ ~ܲሺΘሻ and each subsequent iteration, ݅ ൌ 1, … , ‫ܯ‬ involves two steps: (1) sampling latent variables ሺ‫ܪ‬ ௧ , ‫ܤ‬ ௧ ሻ given the data, D, and the current parameter values Θ ିଵ and (2)sampling new parameter values Θ given the D and the latent variables ሺ‫ܪ‬ ௧ , ‫ܤ‬ ௧ ሻ .The first step is performed using Sequential Monte Carlo (SMC) with ancestral sampling.The second step is further broken down into two parts, first parameters ሺ߬ ைே , ߬ ைிி ሻ are sampled using an adaptive Metropolis-Hastings sampler and then parameters ሺ݇, ݀, ݂ሻ are sampled using Hamiltonian Monte Carlo (HMC).To access the effect of pharmacological interventions on LH pulsatility, hormoneBayes allows the input of two LH profiling datasets, corresponding to control and perturbed conditions.The inference task is then modified to allow inference of parameters Θ ൌ ሺ߬ ைே , ߬ ைிி , ݇, ݀, ݂ሻ for the control case as well of four additional parameters Θ ൌ ሺ‫߬‬ ைே , ‫߬‬ ைிி , ‫,݇‬ ‫݂‬ሻ corresponding to the perturbed case.In Sampling from the posterior is performed as described above.We quantify the effect of the perturbation as the log-ratio of perturbed to control values.