Reduced model for female endocrine dynamics: Validation and functional variations
RReduced model for female endocrine dynamics:Validation and functional variations
Erica J. Graham a, ∗ , No´emie Elhadad b , David Albers c a Mathematics Department, Bryn Mawr College, Bryn Mawr, PA 19010, USA b Department of Biomedical Informatics, Columbia University, New York, NY 10032, USA c Pediatrics Department, University of Colorado Denver–Anschutz Medical Campus, Aurora, CO 80045, USA
Abstract
A normally functioning menstrual cycle requires significant crosstalk between hormones originating in ovar-ian and brain tissues. Reproductive hormone dysregulation may cause abnormal function and sometimesinfertility. The inherent complexity in this endocrine system is a challenge to identifying mechanisms ofcycle disruption, particularly given the large number of unknown parameters in existing mathematicalmodels. We develop a new endocrine model to limit model complexity and use simulated distributions ofunknown parameters for model analysis. By employing a comprehensive model evaluation, we identify acollection of mechanisms that differentiate normal and abnormal phenotypes. We also discover an intermedi-ate phenotype–displaying relatively normal hormone levels and cycle dynamics–that is grouped statisticallywith the irregular phenotype. Results provide insight into how clinical symptoms associated with ovulatorydisruption may not be detected through hormone measurements alone.
1. Introduction
Female endocrine physiology is relatively poorly understood and is related to several high-impact diseases—including polycystic ovary syndrome (PCOS), endometriosis, and diabetes-related illnesses like non-alcoholicfatty liver disease—along with general reproductive health and mental health. Here we are primarily mo-tivated by two of these: PCOS, a condition where menstrual periods are infrequent or prolonged and/orwhere excess male hormone (androgen) levels are present, and endometriosis, a condition where uterine-lining tissue grows throughout the body. Both conditions lead to severe complications and reduced qualityof life. The systems physiologic understanding of such conditions are poorly understood, and few data andfew models exist to explain their origin or to provide a basis for high-fidelity phenotypes.We want to forge the path for better understanding of female endocrine physiology anchored with math-ematical physiology. The pathway for using quantitative descriptions of physiology to impact human healthis non-trivial: (1) physiologic systems are complex; (2) mathematical models of physiology do not rely onfirst principles, but are rather idealizations of how we imagine the physiologic subsystems to function; (3)model verification is also complex, requiring data and a means of evaluating the model’s ability to representthose data; (4) models and parameters must be personalized and synchronized with patient data; and, (5)translation of the model-based information is difficult and complicated.Within clinical practice, the most common usage of physiology is for reasoning by analogy or via qual-itative relationships. For example, in an ICU we inject insulin to reduce dangerously high glucose levels,and this is done incrementally with a dose chosen to be small enough to avoid hypoglycemia. We want tomove beyond this reason-by-analogy approach; but, in addition to the complexity of clinical data, there areroadblocks, such as a lack of data for validation, no direct translation between complex quantitative physio-logic information and the clinical setting, and the subsequent lack of clinical entry points. Nevertheless, the ∗ Corresponding author. Address: Mathematics Department, 101 N. Merion Ave., Bryn Mawr, PA 19010, USA.
Email address: [email protected] (Erica J. Graham)
Preprint submitted to XXXX a r X i v : . [ q - b i o . T O ] J un roblems we choose to tackle, along with the basic construction and usage of related models, can be directedby the real potential to impact our understanding of health, especially as far as informed decision-makingand plausible data collection are concerned.Here we focus on the subsystem of the female endocrine system that is thought to control the ovula-tory cycle. In the qualitative, physiologically mechanistic, mathematical context, descriptions of ovulatorydysfunction are complicated. First, the physiology is complex: there are multiple mechanisms, both in thebrain and in the ovaries, that can alter ovulatory function. Second, and more problematic, we are extremelydata-limited: even qualitatively, there is no systematic, comprehensive way to classify dysfunction becausedynamical time scales are long (months to years), intra- and inter-personal variability is high, and observablemanifestation of dysfunction can have many sources whose delineation can be difficult to resolve with data.The invasiveness of collection procedures adds to the limited data at our disposal.The data that can exist include primary reproductive hormone measurements, which are useful fordelineating broadly defined clinical abnormalities and quantifying generalized ovulatory states. For example,two prototypical data sets reported in the literature include pituitary and ovarian hormones collected dailyover the course of a typical cycle [1, 2]; we use the data in [1] in this paper. The challenge is thesedata may only provide a partial view to more subtle abnormalities. For example, PCOS can result in thecomplete absence of, or sporadic, ovulation. But, distinguishing between mechanisms governing these twoobservable clinical manifestations is difficult because clinically feasible diagnostic tools rely on measurementstaken either at a single time point or over the course of a few hours [3]. We would require data spanningmultiple months in order to build a comprehensive hormone profile with any hope of revealing importantreproductive features, especially in the absence of clearly identifiable ovulatory states. Still more confoundingare conditions such as endometriosis, which lack a clear etiology, yet are linked to circulating hormone levels[4]. In the present context, it is important to note that a high-fidelity, data-driven, robust and expansivedefinition of normal ovulatory function does not currently exist. This makes defining ‘normal’ and ‘dysfunc-tional’ a complex task, as dysfunction is usually defined as a deviation from normal. Because of this, we willadopt a narrow definition of normal and consequently limit our ability to discover a different-from-normalphenotype. Ideally, we seek an alternative to patterns in hormone dynamics to distinguish between ovu-latory phenotypes, with the hope that identifying underlying mechanisms of dysfunction lies in our abilityto connect clinical symptoms with mechanisms that may not be apparent in hormone measurements alone.Given the problems identified above, modeling and analysis at this level of detail is not possible.In this paper we do not attempt to overcome all problems between the development of a model and theuse of the model to help improve human health at once, but rather focus on two. First, we begin with amodel and reduce it, which has two consequences: (1) reduction of identifiable pathologies by decreasingthe number of states and unknown parameters; and, (2) induction of a physiologic hypothesis about whatvariables are important for representing the female endocrine system related to PCOS and endometriosis.The result is a new endocrine model of ovulation. Second, we evaluate the model’s ability to representdata, delineate the time-dependent differences between normal and abnormal cycles given a cycle length,and we examine emergent phenotypes through analysis of the parameter space. Although model evaluationis limited by data availability, we construct the evaluation methodology such that, given more data, a morepowerful, direct evaluation will be immediately possible. In Section 2, we review the Graham-Selgrade model of ovulatory dynamics [5], which is the modelingstarting point. We then reduce this model by removing testosterone to create a new model of ovulatorydynamics; this is our first result. In Section 3, we introduce the computational, evaluation, and data-relatedmachinery we use to validate the model and examine its ability to resolve and differentiate our narrowlydefined phenotypes. In Section 4, we follow with the computational validation of the new model as well asour investigation of the model parameters to study clinical phenotypes.2 . Ovulation Model and its Reduction
We develop a new endocrine model to describe essential processes in ovulation. This new model isa reduction of a model developed by Graham and Selgrade that uses ordinary differential equations todescribe the ovulatory cycle under the influence of elevated androgens, namely testosterone [5]. Testosteroneis a major element in the Graham-Selgrade model because androgen excess, or hyperandrogenism , andinsulin resistance are frequently associated with ovulatory dysfunction related to PCOS. Although PCOSis an important disorder with many open questions regarding its etiology [6], we presently aim to examinegeneralized ovulatory dysfunction, which may or may not stem from previously defined clinical phenotypes,e.g., PCOS. Moreover, while we reduce the model by eliminating testosterone as a state variable, testosteroneremains an implicit variable in the model; we verify testosterone-mediated dysfunction in the course ofvalidating the reduced model. To achieve this goal we work to reduce the model to limit the size ofparameter and state spaces, thereby reducing the complexity of the analysis and the data required to resolvephenotypes. We choose to begin with the Graham-Selgrade model because we deem it more amenable toreductions in the parameter space while retaining normal ovulatory dynamics, in contrast to other relatedmodels that are based on delay differential equations [7, 8, 9].
The Graham-Selgrade model [5] is divided into three major subsystems: pituitary regulation, follicledynamics, and ovarian steroidogenesis. Collectively, the model consists of 12 state variables, tracking serumconcentrations of five important reproductive hormones (follicle stimulating hormone (FSH), luteinizinghormone (LH), estradiol (E ), progesterone (P ), and testosterone (T)), along with precursors/intermediariesof LH, FSH, and T. It also describes the dynamics of three follicular stages and of the follicle response to LH,termed LH sensitivity . The final model contains 41 unknown parameters which are estimated by fitting themodel to data from the literature [1, 10]. The complete list of equations for the original Graham-Selgrademodel may be found in Appendix A. The Graham-Selgrade model uses a compartmental framework to examine changes in ovulation due toincreased androgens. The model follows the approaches of [7, 8, 9] and comprises three major subsystems,which describe changes in the pituitary-ovarian axis with mechanisms of steroidogenesis.I.
Pituitary regulation.
LH and FSH are the primary hormones produced by the pituitary gland. Syn-thesis and release of these hormones are regulated by ovarian steroid hormones, including E , P ,and T. The equations governing changes in FSH and LH are split between releasable (denoted F SH ρ and LH ρ ) and serum (denoted F SH and LH ) pools of the hormones and incorporate stimulatoryand inhibitory feedback by ovarian steroids. Using this compartmental approach, we can differentiatefeedback processes governing pituitary hormone synthesis versus release.Here we provide a generalized description of pituitary dynamics. Let H ( t ) denote the serum concen-tration of a pituitary hormone (either FSH or LH) and H ρ ( t ) its releasable amount at time t . For H = F SH, LH , the differential equations governing releasable and serum quantities have the formd H ρ d t = k synthesis ( · ) − k release ( E , P ) H ρ , (1)d H d t = k release ( E , P ) H ρ /V − δ H H. (2)Each k ( · ) term denotes a function of state variables and describes the change in hormone levels dueto the process indicated. Synthesis of FSH and LH is determined by different processes—with precise The model presented in [5] contains a typographical error in one of the equations, which omits one parameter ( c Φ ,T ) fromthe total parameter count cited. k synthesis omitted to reflect this—whereas their release is mediated solely by E and P .Release into the serum is scaled by the blood volume, V , and clearance of the hormones is assumedto be a first-order process, with rate constant δ H . Regardless of the highly nonlinear form of ovarianfeedback, the subsystem remains linear in H ρ and H . Collectively, the pituitary subsystem comprisesfour differential equations, with Equations (1) and (2) defined explicitly for both FSH and LH.II. Follicle dynamics.
Follicle growth, maturation, and differentiation are assumed to occur in a series ofthree sequential stages: (1) follicular , (2) ovulatory , and (3) luteal . We denote these using variablesΦ( t ) , Ω( t ), and Λ( t ), respectively. The follicular phase is characterized by recruitment and growth ofstimulated follicles. The ovulatory phase is characterized by ovum release from a designated follicle inresponse to a mid-cycle surge in LH. Finally, the luteal phase is characterized by the formation and,in the absence of fertilization, regression of the corpus luteum. The three follicular stages are modeledas follows: dΦd t = k recruitment ( T ) + k growth ( F SH, T )Φ − k ovulation ( F SH, LH )Φ , (3)dΩd t = k ovulation ( F SH, LH )Φ − k luteal ( S )Ω , (4)dΛd t = k luteal ( S )Ω − k regression ( S )Λ . (5)Transitions to subsequent stages are unidirectional and depend on pituitary hormone levels. Themodel also incorporates a role for T in follicle recruitment and growth. Graham and Selgrade furtherdefine a new LH support variable, S ( t ), to model the tonic LH-dependence of growth and prematureregression of the corpus luteum. Specifically, S decays exponentially (with rate δ S ) to 0 in the absenceof LH and approaches a maximal level of 1 for sufficiently large LH:d S d t = k activation ( LH )(1 − S ) − δ S S. (6)III. Ovarian steroidogenesis.
Throughout the ovulatory cycle, follicles may produce E , P , and T. In-tracellular steroid production is primarily FSH- and LH-dependent during a typical cycle and issubject to functional maturation of individual follicles. This subsystem exploits the two-cell two-gonadotropin theory of ovarian steroid production, which describes the differential functionality oftheca cells and granulosa cells within ovarian follicles [3]. The Graham-Selgrade model also introducesa semi-mechanistic description of testosterone production for examining a role for insulin in promotinghyperandrogenism. For T γ ( t ) denoting the ‘intermediate’ concentration of T destined to be convertedinto E , we write d T γ d t = k entry ( LH, α ) − k aromatization ( F SH ) T γ . (7)In a growing follicle, theca cells compose the outermost layers of cells surrounding the ovum andgranulosa cells the innermost layers. Importantly, theca cells possess androgen (i.e. T) productionmachinery and are stimulated by LH alone, whereas only neighboring granulosa cells can convert theseandrogens into estrogens, in an FSH-dependent process called aromatization . Therefore, we consider T γ to reflect the average concentration of T that enters granulosa cells from theca cells.Finally, we model the major ovarian outputs of the model: serum concentrations of E , T, and P :d E d t = k basal,E − δ E E + k aromatization ( F SH ) T γ · f E (Φ , Ω , Λ) , (8)d T d t = k basal,T − δ T T + (cid:20) k ovarianproduction ( LH, α ) + k peripheralproduction ( LH, α ) (cid:21) · f T (Φ , Ω , Λ) , (9)4 P d t = k basal,P − δ P P + k secretion ( LH ) · f P (Φ , Ω , Λ) . (10)The first two terms in Equations (8)–(10) represent basal secretion by the adrenal gland and first-order clearance of individual steroids, defined by rate constants k basal,I and δ I , respectively, where I = E, T, P . The last term in each equation defines secretion of steroid hormones into the circulation,which is assumed to occur immediately upon production. The average production rate per follicleis multiplied by a function f I (Φ , Ω , Λ) , I = E, T, P , that describes the relative contribution of eachfollicular stage to the production of a given steroid.Importantly, steroidogenesis is altered through feedback from FSH and LH, according to the twocell-two gonadotropin theory. Whereas LH is required almost exclusively for T (theca only) andP (theca and granulosa) production, FSH is entirely responsible for E (granulosa only). BecauseP is an androgen precursor in the theca, it is assumed that circulating P is produced primarilyby granulosa cells for modeling purposes. To address insulin’s influence in ovulatory dysfunction,the Graham-Selgrade model contains a detailed formulation of T production, wherein ovarian andperipheral conversion of T from its precursors are treated as two distinct processes. In Equations (7)and (9), the parameter α represents the relative degree to which insulin may increase T production. A natural course of action in determining salient model behavior is global sensitivity analysis (GSA)of parameters. This allows us to determine the relative sensitivity of model output to changes in theparameters. There are multiple challenges associated with the Graham-Selgrade model that make GSA asuboptimal next step in model analysis. First, the model contains considerably more parameters than thedata available for fitting. Second, coupling between state variables is highly nonlinear. Third, stable limitcycles are not guaranteed for all parameter combinations. Collectively, standard GSA approaches providelimited insight. In particular, a PRCC-based approached would be inappropriate, as simulations do not yieldmonotonic hormone responses that can be interpreted in any meaningful way (preliminary work, not shown).Alternatives such as the extended Fourier amplitude sensitivity test (eFAST) may also prove more useful,as discussed in [11]. However, selection of appropriate model output remains a challenge. In light of theseobservations, we hypothesize that structural reduction of the model may provide greater insight to relevantand essential processes governing the typical ovulatory cycle. We may also use the resulting framework toexamine more general questions of ovulatory function separate from the pathologies associated with PCOS.Here we introduce the major modifications to effectively reduce the number of unknown model parameters.
Complete data sets that track the pituitary hormones, LH and FSH, as well as ovarian steroids E andP during the course of an entire cycle are uncommon but not completely absent. What is missing is acomplete hormone profile that also includes androgen levels through the course of a normal ovulatory cycle.Based on this information, we eliminate T (and hence insulin) entirely from the Graham-Selgrade model.With this adjustment, we remove Equation (9) entirely and set α = 0 and T ( t ) = T for all t . We thenadjust the remaining differential equations as needed to eliminate T (and T γ ) coupling. Effect on steroidogenesis.
We assume that intermediate T transferred from theca to granulosa cells is im-mediately converted into E . Then using a similar reduction approach in [5] applied to a more mechanisticsteroidogenesis model, we set ˙ T γ = 0 and solve for T γ . Substituting the resulting expression and theoriginal parameters into Equation (8) gives ˙ E = e − δ E E + t g F ( LH ) · (Φ + η Λ S ) , where F ( LH ) = LH / [ κ LH + κ LH + κ ]. We further observe that for sufficiently large LH , κ (cid:28) κ LH + κ LH , sowe redefine F ( LH ) = LH/ ( LH + κ ) and rescale parameters accordingly (see Equation 19).Notably, the reduced version of ˙ E removes the FSH-dependence on E production from the originalframework. For a normally ovulating cycle, we may then consider FSH to be permissive for E synthesiswithin granulosa cells. Further, the time scale on which FSH alters follicular expression of aromatase—viathe steroidogenic acute regulatory protein— is significantly shorter ( ∼ minutes) than that of a typical cycle512]. As such, we assume any tonic level of FSH allows for proper steroid synthesis, and the degree towhich this occurs depends solely on the functional maturation of follicles in the reduced model. Therefore,we relegate FSH-dependent E production to a nonessential process and focus instead on other sources ofpathological behavior within the cycle. Effect on pituitary regulation.
In the original model, the influence of T is restricted to synthesis of releasableLH: T is assumed to increase basal LH production and to prevent P -mediated inhibition of all LH produc-tion. To eliminate T from the necessary terms, we redefine the affected parameters given our assumption ofrelatively constant T levels. In particular, we set v new0 L = v old0 L T / ( K L,T + T ) and K new iL,P = K old iL,P (1 + c L,T T )(see Equation 13). Effect on follicle dynamics.
In the original model, T serves two functions: (1) it influences the rate atwhich very immature follicles enter stage Φ to begin gonadotropin-dependent growth and differentiation,and (2) it increases follicle sensitivity to FSH signaling. To eliminate T from this subsystem, we notethat the originally estimated basal rate of T-mediated follicle recruitment is f ∼ O (10 − ). It is thereforeunlikely that this process contributes substantially to the function of a normal cycle, and so we simply set f = 0. Following our approach in the pituitary subsystem, we redefine the FSH sensitivity parameter to be h new1 = h old1 / (1 + c Φ ,T ), provided that follicular FSH sensitivity remains constant throughout the menstrualcycle (Equation 15). In addition to eliminating T as a state variable, we make two simplifying assumptions to further reducethe number of unknown parameters. These parameters are chosen in consideration of the important biologicalcomponents that must be maintained in order to consider the resulting model an accurate and usefulrepresentation of the menstrual cycle under physiological conditions. Rather than focusing on a topologicallyequivalent system, we focus on preserving plausible biological mechanisms.
FSH-dependent LH sensitivity.
An essential event in ovulation is the upregulation of LH receptors in thelate follicular stage (stage Φ). This is an FSH-dependent process occurring within sufficiently maturefollicles. The Graham-Selgrade model assumes FSH increases follicle sensitivity to LH during stage Φ. Inthe model reduction, we assume the maximal sensitivity parameter, h remains constant. This is a reasonablesimplification, as h ∼ O (10 ) and c Φ ,F · max t { F SH ( t ) } ∼ O (1) in the original model. LH-dependent P production. P conversion within theca and granulosa cells requires enzymes that areregulated by LH. The estimate for the half-maximal LH concentration, h p , that stimulates P production inthe original model is lower than the simulated LH concentration during the luteal phase (where P attainspeak concentration). Therefore, we assume that the steroid production per follicle is constant at the maximalrate p in LH. The new model is given by Equations 11–20 and contains 10 differential equations. With 27 unknownparameters, we have reduced the parameter space by more than a third. Terms in the model that havebeen altered due to removal of testosterone are boxed with a single line. Those that result from additionalsimplifying assumptions as described in Section 2.2.2 are boxed with a double line. For comparison, theoriginal model equations are listed in Appendix A.
Releasable FSH: d F SH ρ d t = v F c F,I S Λ K F,I + S Λ − k F c F,P P c F,E E F SH ρ (11) Serum FSH: d F SH d t = 1 V · k F c F,P P c F,E E F SH ρ − δ F F SH (12)
Releasable LH: d LH ρ d t = (cid:20) v L + v L E n K nmL + E n (cid:21) ·
11 + P /K iL,P − k L c L,P P c L,E E LH ρ (13)6 erum LH: d LH d t = 1 V · k L c L,P P c L,E E LH ρ − δ L LH (14) Follicular phase: dΦd t = f F SH h + F SH − f LH h + LH · Φ (15)
Ovulatory phase: dΩd t = f LH h + LH · Φ − wS Ω (16)
Luteal phase: dΛd t = wS Ω − l (1 − S )Λ (17) LH support: d S d t = ˆ s LH LH + h s (1 − S ) − δ S S (18) Serum E : d E d t = e − δ E E + t g LHLH + κ · (Φ + η Λ S ) (19) Serum P : d P d t = − δ P P + p Λ S (20)
3. Computational Methods and Model Evaluation
To discuss model evaluation and results, we explicitly distinguish between physiological and mathemat-ical notions of a ‘cycle’. When discussing properties of mathematical ovulation, we explicitly refer to theinter-ovulatory interval (IOI), which denotes the length of time between consecutive simulated LH surges.Physiologically, the IOI is equivalent to the time between two ovulatory cycles ; however, multiple IOIs maybe required before the solution completes a single mathematical (limit) cycle. For clarity and consistency,we restrict our generalized use of ‘cycle’ to refer to physiological ovulation and IOI to the calculated timesbetween these cycles.
We use two data sets, one synthetic and one real. The first data set, the synthetic data set, is generatedusing the original model [5]. This data set is used to show that the reduced model captures most of thedynamics of the original model. The second data set is hormone data available in [1]. These data containaverage daily measurements for 33 normally cycling women during the course of one complete ovulatorycycle for FSH, LH, E , and P . This second data set is used to demonstrate the both ability of the modelto estimate data well, and how to use the model to better understand physiology, given data. The available data have three primary limitations that influence our work. First, recall that normalis generally poorly defined, where ‘normal’ means no known pathophysiologic cycle features. Second, it isknown that there is substantial variation in IOIs even for an individual. For example, it is not uncommonfor the same person to have IOIs that vary from 20 to 40 days; these data obscure such intraindividualvariability by taking an average. And third, because the data are an average, they induce two potentialissues whose presence we may not be able to detect: (i) an average can fail to represent anyone if the meanis not representative of the population, and (ii) an average smooths variability.These data limitations impact the analysis in a fundamental way, most notably, the generalizability ofour results. We develop a phenotypic analysis subject to a standard IOI of 31 days and relative to averagehormone dynamics. It is surely possible that the average data do not represent individuals well. Moreover, itis surely possible that hormone dynamics of different IOIs will be different. Finally, within the model, thereare two ways of generating variable IOIs: (1) one can change parameters that alter a constant IOI, or (2)7ne may define parameter regimes that allow for variable IOIs between consecutive ovulatory events. Giventhat our data set is limited to a month and is an average, we cannot investigate the distinction betweenthese two model-based parameter differences. In short, the data set limits the generalizability of data-basedmodel validation, but our demonstration for how this pathway will work in the future remains relevant.
Equipped with computational machinery for estimating parameters given data, we want to demonstratethe capability of identifying and defining phenotypes that emerge from the model. As previously mentioned,we cannot move beyond what can be tested with data we have, but for our purposes this will not limit us.We focus on two model-defined phenotypes, physiologic and pathophysiologic, and codify these as regularor irregular cycle behavior, respectively. We assign these to cases of model-generated data by defining aset of attributes that may distinguish between physiological and pathological ovulatory function. Given thecomplex cross-talk in the reproductive endocrine system, analysis of hormone concentrations alone likelyprovides insufficient insight into the subtleties of ovulatory dysfunction. To overcome this challenge, weproceed by identifying a collection of parameters giving rise to predetermined regular or irregular cyclebehavior. To accomplish this, we implement an algorithm that allows us to carry out a comprehensiveevaluation of the reduced model.
We introduce an algorithm to compare the new endocrine model to the original Graham-Selgrade frame-work. In this section, we provide an overview (see Algorithm 1) and a detailed description of our implemen-tation of the five-step algorithm.
Step 1: Synthetic data.
We generate the synthetic data by numerically solving the system (A.1)–(A.12),using the parameters in [5], for a sufficiently long time to approach a stable limit cycle for normalovulation. We then align the trajectories so that the LH surge occurs at the end of day 15 of the firstcycle. Finally, we extract daily data between days 0 and 30 and then, to avoid propagated numericalinaccuracies, repeat the cycle twice more for each variable. We also expand the set of data to includeΦ, Ω, and Λ, under the assumption that follicular dynamics should follow a similar pattern to theoriginal model. With the exclusion of T from the model, we have a total of 7 state variables (includingFSH, LH, E , and P ) with n = 93 data points each. Step 2: Optimization.
To determine how well the new model compares to the original model, we estimatethe 27 parameters of the reduced model by fitting output to the synthetic data. We capture essentialcycle behavior with the optimized parameters using a weighted least squares approach. For a givenvariable X i ( t ), where i ∈ { F SH, LH, E , P , Φ , Ω , Λ } , we first assign default weights w i = 1 / Var( X i ( t ))to each data point at t j = 0 , , . . . ,
92. Because we cannot guarantee the expected behavior of follicular
Algorithm 1
Comprehensive model evaluation.
Step 1 . Generate synthetic data set using Equations (A.1)–(A.12).
Step 2 . Optimize reduced model parameters using weighted least squares and synthetic data.
Step 3 . Run N Monte Carlo simulations, initialized with perturbed best-fit parameters from Step 2 andrefit to clinical data.
Step 4 . Compute numerical solutions for each parameter profile generated in Step 3, and store resultinghormone data over multiple cycles.
Step 5 . Use results of Steps 3 and 4 to define salient phenotypes and distributions for each of the 27reduced model parameters. 8ynamics, we do not incorporate additional time-dependent weights for i = Φ , Ω , Λ. However, for thehormones we increase weights by variable factors at important peaks, troughs, and plateaus within thedata. These weights are adjusted to acquire the best qualitative fit to the data, with the understandingthat local minimization of the cost function may be sensitive to variation in weights and may notproduce a globally optimal solution.Let y i represent the vector of measurements corresponding to reduced model output variable x i ( ϕ ),defined by parameters ϕ . We define the optimization problem that minimizes the sum of the squarederror as min ϕ | V | · n (cid:88) i w i || y i − x i ( ϕ ) || , (21)where V = { F SH, LH, E , P , Φ , Ω , Λ } and denote the optimal parameter vector satisfying Equation(21) by ϕ ∗ . We use Matlab ’s fminsearch , which implements the Nelder-Mead simplex method, todetermine the optimal ϕ . In most cases, initial parameter guesses are taken from the original model.In others, they are derived from the adjustments made in the reduction process, as described in Section2.2. The best-fit parameters are listed in Table B.2 in the Appendix. Step 3: Monte Carlo simulations.
We determine the distributions of the 27 model parameters usingMonte Carlo simulations to generate a collection of best-fit parameters using various initial guesses inthe estimation scheme described in Step 2 and by comparing the output to data. We first assume thatthe values in ϕ ∗ represent mean quantities and that initial guesses, ϕ (0) , are uniformly distributedwithin ±
10% of the mean. That is, ϕ (0) k ∼ U (0 . ϕ ∗ k , . ϕ ∗ k ) for k = 1 , , . . . ,
27. To ensure a represen-tative sampling of N parameter combinations from each individual subinterval of length 1 /N rangingfrom 0 . ϕ ∗ k to 1 . ϕ ∗ k , we use Latin hypercube sampling (LHS) and randomly generate initial parame-ter guesses for the Monte Carlo simulations (see [13] for further discussion on LHS). For each initialparameterization we minimize a cost function similar to Equation (21), this time with measurements y i , i ∈ { F SH, LH, E , P } taken from clinical data available in [1]. Step 4: Range of simulated model output.
We numerically solve the reduced model over 93 days usingthe estimated parameters and generate an ensemble of these solutions via Monte Carlo sampling. Wealign each LH surge (assuming one exists) to day 15 and determine the length of each IOI. The
LHsurge is defined to be a peak LH concentration that is followed by an apparent luteal phase; any otherlocal maxima in LH failing to meet this criterion are ignored. Because we have restricted our samplingscheme in Step 3, we guarantee that the model does not approach a stable equilibrium. Although thislimitation does not capture complete ovarian failure (i.e., the absence of a cycle at all), it does allowfor reasonable comparisons in the presence of oscillatory dynamics.
Step 5: Phenotypes and parameter distributions. (a)
Defining phenotypes.
We implement a two-step process for determining distinct phenotypesusing model output. First, we use the values of extreme IOIs to ensure that the presence ofabnormally long or short IOIs at any time is considered pathological. In particular, we assign a regular phenotype to simulations resulting in both minimal and maximal IOIs between 25 and35 days, which is the textbook standard range for normal ovulatory cycles [3]. We assign an irregular phenotype to simulations failing to satisfy this criterion. Second, we compute the meansquared error (MSE) between the data (LH, FSH, E , and P ) and each simulation. Then weuse the minimal MSE attained by an irregular phenotype to define a threshold for secondaryregular phenotypes: regular + lo MSE refers to regular phenotypes with MSE strictly less thanthe computed threshold, and regular + hi MSE to regular phenotypes with MSE at or above thecomputed threshold.(b) Computing parameter distributions.
We construct empirical parameter distributions basedon the optimized parameter sets obtained from the Monte Carlo simulations. First, we apply theprimary (regular vs. irregular) phenotype classification criteria to the N Monte Carlo samples.9hen we normalize the population sizes of individual phenotypes by subsampling the associatedparameter distributions at their respective frequencies to an arbitrarily chosen size of 2000.
The addition of a phenotype classification generates several interesting questions that may be exploredwith the use of statistical and probabilistic tools.
Two-sample Kolmogorov-Smirnov (KS) test.
In the present work, we seek phenotypic differences determinedby model parameters. The KS test is used to determine whether two samples are drawn from the samedistribution [14, 15]. The test uses the Kolmogorov-Smirnov statistic, which is defined as the L ∞ norm ofthe distance between two cumulative probability distribution functions. For each parameter, we then apply ks.test , the R implementation of the two-sampled KS test, to analyze the phenotype-specific empiricaldistributions generated from our simulations. t-Distributed stochastic neighbor embedding (t-SNE). Beyond the structure manually imposed on the MonteCarlo dataset, we are interested in determining whether distinct phenotypes can be identified in anotherway. Patterns in the generated data may depend on any of 93 data points for each of four hormones, or anyof the 27 parameter estimates. Without a comprehensive understanding of the interplay between each ofthese elements, we seek a methodology that will answer the binary question of whether there are inherentdifferences (seen or unseen) between regular and irregular phenotypes. t-SNE is a machine learning tool forreduction of high-dimensional data to lower dimensions [16]. We use the
Rtnse package in R to apply thet-SNE and determine whether phenotypes can be clearly clustered by a profile of select model parameters.
4. Computational Results
Following Steps 1 and 2 of Algorithm 1, we simulate the reduced model and compare results to theoriginal model. We then use the parameterized model to simulate testosterone-mediated dysfunction, asfurther verification that the reduction is a plausible replacement of the original system.
Qualitative features.
In Figure 1 we numerically solve the reduced model using the best-fit parameters andcompare the result to output from the original Graham-Selgrade model. The qualitative dynamics are wellcaptured, with the primary quantitative discrepancy related to P . This arises due to an overshoot of thedata in the luteal stage Λ during the mid-luteal stage (roughly 3–7 days after the simulated preovulatory LHsurge, not shown). Since P levels are known to peak clinically around this time, we consider this behaviorto be within a physiologically relevant and normal range for the hormone. Further, because we assumethat the ovarian stages are crude approximations to actual follicular dynamics, there may be substantialvariability in the trajectories that may nevertheless yield normal ovulatory function, as illustrated in [10]. Verification of testosterone-mediated dysfunction.
A fundamental change in the reduced framework is theomission of T. Although absent from the model, we may still examine how T might influence pathologicalovulation. This approach also serves as proof of concept when using the reduced model in lieu of the originalone. To re-incorporate T into the present framework, we modify relevant parameters. Following [5], we let α denote the degree of insulin influence, where α = 0 reflects a normal state with basal insulin (and henceT) levels. Assuming T remains constant over time, we define its concentration using a linear function in α ,denoted T α : T α = T · [1 + ( δ T − · (1 + α )] /δ T , (22)where T is the initial T concentration in the absence of hyperinsulinemia and δ T is the first-order clearancerate of T from the blood, as defined originally. The parameters to be altered by T in the reduced modelare v L , K iL,P , and h . We only consider the case of normal luteinization (see [5] for details) because we10 lllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll FSH ( IU L ) llllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll LH ( IU L ) llllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll E ( pg mL ) Time (days) llllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll P ( ng mL ) Time (days)
Figure 1. Fit of reduced model to Graham and Selgrade model [5] over 61 days. FSH and LH are displayed in standardinternational units according to the 2 nd international reference preparation, where 1 IU FSH = 45 µ g and 1 IU LH = 15 µ g[17]. Conversion factors are based on the NIH preparation used in [1]. have omitted FSH-dependent upregulation of follicle LH receptors, which would impact parameter h . Toincorporate the necessary modifications to the current model, we redefine the parameters v L → v L ξ , K iL,P → K iL,P ξ , and h → h ξ for α >
0, where ξ = ( β + T ) · T α ( β + T α ) · T , (23a) ξ = 1 + β T α β T , and (23b) ξ = 1 + β β T α /T . (23c)The ξ i in Equations (23) determine the scaling of the model parameters as insulin influence increasesand are plotted in Figure 2. The constants β i are defined according to the original model, with the caveatthat bifurcation values of α may be shifted based on the values of these parameters. The derivation of the ξ i are given in Appendix C. In Figure 3, we plot the long-term local maximum and minimum values ofLH for 0 (cid:54) α (cid:54)
3. A stable limit cycle is roughly evident for α < .
2, with an apparent period doublingbifurcation giving rise to alternating LH surge amplitudes. Minimal LH levels remain relatively constantuntil another apparent bifurcation occurs at α ≈ .
7, where the system approaches a stable equilibrium.This suggests that the reduced framework responds to elevated T through modified and eventual loss ofovulatory function. Although the dynamic mechanisms governing ultimate dysfunction may differ from theoriginal model, we are able to capture disruptive behavior. Whereas LH and P are primary determinantsof disruption in the original model, subject to sustained oscillations, the reduced model exhibits disruptionthrough loss of periodicity. Nevertheless, there is some consistency among parameters v L , K iL,P , and h in regulating endocrine defects. From Steps 3 and 4 of Algorithm 1, we obtain an ensemble of model trajectories. Figure 4 showsprimary hormone trajectories over 186 days for two model solutions, one regular and one irregular, asdefined in Step 5(a) of the evaluation algorithm. For reference, the timing of the LH surge for the regularphenotype is indicated with a vertical line. Stable limit cycle behavior is exhibited for the regular cycle with acharacteristic length of 30.9 days. The irregular phenotype, however, consists of nonuniform behavior of themajor hormones. Specifically, the irregular cycle has a length of 80.7 days, with 19.5 and 61.2 days passingbetween consecutive LH surges. Although hormone levels are relatively normal through the course of the11 x x a Figure 2. Dimensionless functional forms used to incor-porate T into reduced model, as in Equations (22) and(23). Each ξ i contributes a T-dependent change (per-cent increase or decrease) in relevant parameters from theoriginal model [5]. ξ increases LH synthesis parameter v L , ξ increases P -mediated LH inhibition parameter K iL,P , and ξ descreases FSH sensitivity parameter h . α : degree of insulin influence. a L H ( I U / L ) PD HB
Figure 3. Simulated bifurcation diagram depicting adjusted rolefor T and insulin influence ( α ). Maximal and minimal LH con-centrations are shown for various values of α (cid:62)
0. For α (cid:47) . α > . α ≈ .
7, the system approaches an equilibriumwith an apparent Hopf bifurcation (HB). F S H ( I U L ) L H ( I U L ) LH surge(regular) E ( pg m L ) P ( ng m L ) Figure 4. Comparison of representative regular and irregular trajectories simulated by the reduced model. The regular cycledisplays a characteristic length of 30.9 days. The irregular cycle has a total length of 80.7 days, with IOIs of 19.5 and 61.2days.
Using the results from Step 5(b) of Algorithm 1, we can use the Kolmogorov-Smirnov test to assesswhether each parameter distribution differs from its counterpart in the opposing phenotype. Results forsubsampled parameter distributions are illustrated in Figure 5. Each box is shaded according to the minimallevel of significance that allows us to accept the alternative hypothesis, i.e. that regular and irregulardistributions are statistically different. Darker shaded squares correspond to higher levels of significance.Of the 27 parameters remaining in the reduced model, we identify eight that have significantly differentdistributions between regular and irregular phenotypes, with p < .
01 (indicated by ∗ ). These parametersare given in Table 1. Our remaining analysis focuses on these eight important parameters. To determine whether a correlation exists between important parameters and the accuracy of their ac-companying numerical solutions when fit to clinical data, we calculate the mean squared error between themodel output and the measured data. Although we fail to demonstrate a clear mechanistic relationshipbetween any of the eight relevant parameters (not shown) and their impact on hormone dynamics or phe-notypes, we do observe a threshold MSE value—estimated from the MC output—above which all irregularphenotype results lie and below which roughly 85% of regular results lie. We use this threshold to assign anadditional subcategory to simulations belonging to the regular phenotype; solutions are either ‘lo MSE’ or‘hi MSE’ when yielding MSE values below or above the computed threshold, respectively. It is importantto recall that there does exist a subset of parameters for which the IOI varies by 50%, where both regularand irregular IOIs are observed yet the limit cycle length is fixed. Because of this, lo MSE implies both lowintra-cycle hormone variability compared with data and also low IOI variability.In Figure 6, we compute 95% confidence intervals of simulated hormone concentrations over four monthsto examine how hormone profiles influence these refined phenotypes. As before, we align the simulated LHsurge of the first cycle at day 15. Low MSE simulations exhibit the least variation across all cycles (greenregions). Beyond the first LH surge, high MSE regular phenotypes (left panel, teal regions) have morevariation in the timing of characteristic ovulatory events (e.g. LH surge and luteal formation) than lowMSE, but considerably less variation than the irregular phenotypes (right panel, gray region). As a result,predictability of ovulation is reduced when we refine phenotypes.On the other hand, if we are less strict with our definition of ‘normal’, we find that it is more difficultto discern reproductive phenotypes. In Figure 7, we plot the 95% confidence intervals for all LH andP trajectories satisfying the criterion that at least one IOI is between 30 and 32 days long (teal). Forcomparison, we also include the 95% CI for applicable irregular trajectories (gray). Limited to informationon a single IOI, there is considerable overlap between opposing phenotypes, which may obscure our abilityto discern irregularities in hormone regulation. These results are important because they highlight how insufficient data can both mask ovulatory dysfunction and obscure phenotype definition, discovery, andanalysis. −16 −5 −4 −3 −2 −1 p−value * * * * * * * * v F K i F I k F c F P c F E c F I v v K m L K i L P k L c L P c L E f f h h w l s d S h k h s t g e p Figure 5. Two-sample Kolmogorov-Smirnov test. Shaded according p -value. ∗ p < . L H reg. + lo MSE reg. + hi MSE0510152025 0 15 30 45 60 75 90 105 120 P F S H E Time (days) 0510152025 0 15 30 45 60 75 90 105 120 L H irregular0510152025 0 15 30 45 60 75 90 105 120 P F S H E Time (days)
Figure 6. 95% confidence intervals of reduced model output over four regular cycles. Regular + lo MSE (green) compared to( left ) regular + hi MSE and ( right ) irregular phenotypes. Time-dependent regular + lo MSE means are indicated with blackcurves. able 1. Eight parameters identified as most important based on the Kolmogorov-Smirnov test. Parameters are ranked inorder from most (1) to least (8) significant, according to the p -value obtained. Rank Name Description1 η luteal E production;2 v F maximal FSH synthesis rate;3 h follicle sensitivity to FSH;4 ˆ s LH support maximal growth rate; Rank Name Description5 K mL half-maximal E stimulation level;6 δ s LH support decay rate;7 l maximal luteolysis rate;8 f maximal follicle growth rate. L H P Time (days)95% CI (all) irregular 95% CI
Figure 7. 95% confidence intervals (teal) of LH and P trajectories satisfying the criterion of at least one inter-ovulatoryinterval of 30–32 days. Light gray: 95% confidence interval for irregular phenotypes satisfying IOI criterion. Dark gray:maximal irregular % CI that fits within total CI (43% shown). In Figure 8, we examine the distribution of median IOIs for each phenotype. Median IOIs (mean ± s.d.)for the regular phenotype with lo and hi MSE are 30 . ± .
09 and 31 . ± .
72, respectively, whereas themedian irregular length is 32 . ± .
57 days. As with hormone levels, we observe little variation among IOIsin the lo MSE category and a significantly wider spread for the other categories. In terms of mathematicallyversus clinically cyclic behavior, we find that although most simulations result in oscillations, many do notexhibit limit cycle behavior with a characteristic IOI.
To distinguish between refined phenotypes based on parameter estimates, we implement a t-SNE of theparameter profiles according to the assigned phenotypes. We again limit our analysis to the eight significantparameters found in Section 4.1.3. As such, we aim to reduce an eight-dimensional system to two dimensionsand visualize relative distances between individual objects. In Figure 9, we find that the algorithm producestwo distinct clusters. Irregular cycles and regular cycles yielding a high MSE are considered to be closerto each other than to completely regular cycles with low MSE. Importantly, the MSE was not used as adimension in our analysis using t-SNE. This suggests that even with relatively normal ovulatory function, wecan cluster individuals exhibiting some level of abnormality based on a small number of model parameters.15 lllllllllllllllllllllllllllllll llllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll lllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll
20 25 30 35 40 45 5020003000400050006000 l reg. + lo MSEreg. + hi MSEirregular M SE Median inter−ovulatory interval (days) l
25 30 35 40Median IOI
Figure 8.
Main:
Computed mean squared error between data and MonteCarlo simulations vs. median inter-ovulatory intervals (IOI) for primaryand secondary phenotypes. Threshold between low and high errors is in-dicated by a horizontal dashed line.
Inset:
Box-and-whisker summary ofmedian IOI distributions by phenotype. Whiskers indicate the completerange of cycle medians, and box widths the middle 50% of the medians. l llll lll l lll l lll l lll ll l ll lll lll lll lllll ll ll l ll ll l llll lllll l ll lll lll lll ll l lll l ll llll lll lll ll llll ll l ll l lll ll l ll ll l ll l l ll ll l ll ll lll lll l ll lll ll l llll l ll ll l lll l llll lll ll ll ll ll lll lll llll ll ll ll ll lll ll ll l ll ll ll l ll l l llll lll llll lll llll lll l l l ll ll l llllll l ll l lll lll l ll lll ll lll lll l lll −20 −10 0 10 20−20−1001020 dimension 1 d i m en s i on l reg. + lo MSE reg. + hi MSE irregular Figure 9. t -Distributed Stochastic Neigh-bor Embedding of model results. Di-mensional reduction of identified phenotypesbased on the eight significant parameters η, v F , h , ˆ s, K mL , δ s , l, f gives atwo-dimensional embedding of model output.Specifically, distinct clusters emerge accordingto secondary phenotypes, which are separatedaccording to MSE levels.
5. Discussion
We introduce a new, reduced endocrine model that inherently demonstrates both regular and irregularphenotypes classified by the timing of ovulation. The model produces distinct phenotypes as a result ofaltered time-independent parameter regimes and in the absence of disease-specific factors, e.g. testosterone-mediated dysfunction in PCOS. Through a comprehensive model evaluation algorithm, we identify a subsetof model parameters that provide insight into physiological mechanisms of dysfunction. Further, the reducedframework provides a testable hypothesis of model prediction: consistently similar inter-ovulatory intervals(IOIs) between individuals likely reflect similar reproductive hormone dynamics. But, such consistency isa limiting factor in our ability to broadly predict ovulatory function and highlights the fact that a smallset of parameters can produce large variations in ovulatory profiles. These results also imply that there ispotentially a many-to-one relationship between endocrine states, e.g., physiologic parameters, and observ-able endocrine dynamics and dysfunction, e.g., hormone dynamics. This fuzzy causation is not uncommonin physiologic systems or in biomedicine broadly; but to develop better clinical treatment, it is critical tominimize the number of potential causes of an observable problem while maximizing the understanding ofthe physiologic mechanics driving endocrine dynamics. While we further clarify these issues below—we iden-tify potentially testable mechanisms that drive different endocrine dynamics and phenotypes—substantialproblems remain.Based on the most significant parameters identified by the present work, the model highlights mechanismsassociated with pituitary hormone synthesis ( v F , K mL ), follicle growth ( h , f ), luteal dynamics (ˆ s , δ s , l ),and ovarian E production ( η ). However, the redundancy in the biological processes associated with theseparameters allows us to more succinctly characterize sources of dysfunction based on two major processes:altered follicular growth and feedback associated with E concentrations. In vitro experiments suggest that granulosa cells may be more sensitive to FSH in PCOS, affecting bothfollicle growth stimulation [3]. Follicular growth is stimulated by FSH, and the maximal FSH synthesis rateparameter modulates pituitary stores of FSH. In the irregular phenotype, there is tendency toward increasedFSH synthesis. Since FSH is coupled to the ovarian subsystem only through the follicular stage, changesin FSH synthesis will primarily affect follicle growth in the reduced model. The present model reflects suchdysfunction, as the irregular phenotype is associated with both increased follicle growth rate and follicle16ensitivity to FSH.Variations in E are implicated in multiple manifestations of ovulatory dysfunction. For example, de-creased E is characteristic of menopausal women. Prolonged exposure to elevated E has been associatedwith ovulatory disruption in previous mathematical models [8, 18], and elevated E formation has beenfound in in vitro PCOS models [3]. In addition, Higher E responsiveness In the current work, parametersassociated with luteal stage dynamics are altered in the irregular phenotype. In particular, appearance anddisappearance rates of LH support are increased and decreased, respectively. This supports greater ovarianmass during the luteal phase, which contributes to significantly elevated E during this period. Simulatedirregular cycles are also associated with higher E production rates from functional luteal cells and increasedpituitary sensitivity to E , which can prematurely trigger the LH surge. Elevated subthreshold E prolongssuppression of FSH and LH release into the serum, thereby inhibiting follicle growth. In extreme cases, thisresults in two ovulation events close together, followed by an increased period of ovulatory suppression. Thisis exhibited in Figure 4, with a two-month lapse between ovulation events in the representative irregularphenotype.In the absence of testosterone-mediated modifications, all phenotypes in the new endocrine model exhibitsuccessful ovulatory events, albeit with varying frequencies. Further, the hormone concentrations arisingfrom irregular cycles lie within their respective physiological ranges. Interestingly, the range of IOI forirregular phenotypes are consistent with the ranges reported for individuals near menarche or approachingmenopause [3]. Although the ovulatory defects in the present framework do not reflect those specific todisorders like polycystic ovary syndrome, our results suggest that subtle mechanistic differences contribute toaltered ovulatory function and that hormone measurements alone are insufficient in identifying abnormalities.It also appears that our ability to identify defects via reproductive hormones depends on the samplingfrequency of data.The reduced framework is amenable to modifications allowing us to explore T-mediated ovulatory dys-function. Clinically, it remains unclear how disruptions propagate in the face of hyperandrogenism. We findthat when we alter pituitary-specific processes—particularly with respect to LH production—and folliclegrowth processes with linearly increasing levels of T, cyclic behavior ceases. Further, the steady state ap-proached for sufficiently large insulin influence includes a clinically low level of LH. In contrast, LH is oftenfound to be elevated in PCOS populations, but with high interindividual variability. These results suggestthat we may not associate the T-mediated disruptions within the reduced framework with specific PCOSsymptoms, but rather as part of a more generalized manifestation of ovulatory dysfunction due to abnormalresponses in the pituitary-ovarian axis.The over-arching goal is to use models for predictive decision support and to deepen our understandingof physiology. We wish to not only understand mechanisms of function but also the factors that differen-tiate those mechanisms. Endometriosis and PCOS are two high-impact disorders governed by physiology,both with incompletely understood etiologies. We wish to shed insight on these disorders to better informintervention and treatment decisions. The current model and evaluation process allows us to delineatedysfunction based on physiology. As constructed, the model is flexible enough to allow us to highlightimportant—generalizable or disorder-specific—mechanisms of dysfunction. ReferencesReferences [1] R. I. McLachlan, N. L. Cohen, K. D. Dahl, W. J. Bremner, M. R. Soules, Serum inhibin levels during the periovulatoryinterval in normal women: relationships with sex steroid and gonadotrophin levels, Clinical endocrinology 32 (1) (1990)39–48.[2] C. K. Welt, D. J. McNicholl, A. E. Taylor, J. E. Hall, Female reproductive aging is marked by decreased secretion ofdimeric inhibin, The Journal of Clinical Endocrinology & Metabolism 84 (1) (1999) 105–111.[3] J. F. Strauss, R. L. Barbieri, Yen & Jaffe’s Reproductive Endocrinology E-Book: Physiology, Pathophysiology, and ClinicalManagement, 7th Edition, Elsevier Health Sciences, 2013.[4] L. Lode, M. Often Sveen, M. Rudnicki, Abnormal pathways in endometriosis in relation to progesterone resistance: areview, Journal of Endometriosis and Pelvic Pain Disorders 9 (4) (2017) 245–251.
5] E. J. Graham, J. F. Selgrade, A model of ovulatory regulation examining the effects of insulin-mediated testosteroneproduction on ovulatory function, Journal of theoretical biology 416 (2017) 149–160.[6] A. S. Caldwell, M. C. Edwards, R. Desai, M. Jimenez, R. B. Gilchrist, D. J. Handelsman, K. A. Walters, Neuroendocrineandrogen action is a key extraovarian mediator in the development of polycystic ovary syndrome, Proceedings of theNational Academy of Sciences 114 (16) (2017) E3334–E3343.[7] P. M. Schlosser, J. F. Selgrade, A model of gonadotropin regulation during the menstrual cycle in women: Qualitativefeatures, Environmental health perspectives (2000) 873–881.[8] L. H. Clark, P. M. Schlosser, J. F. Selgrade, Multiple stable periodic solutions in a model for hormonal control of themenstrual cycle, Bulletin of mathematical biology 65 (1) (2003) 157–173.[9] A. O. Hendrix, C. L. Hughes, J. F. Selgrade, Modeling endocrine control of the pituitary–ovarian axis: Androgenic influenceand chaotic dynamics, Bulletin of mathematical biology 76 (1) (2014) 136–156.[10] C. C. Keefe, M. M. Goldman, K. Zhang, N. Clarke, R. E. Reitz, C. K. Welt, Simultaneous measurement of thirteensteroid hormones in women with polycystic ovary syndrome and control women using liquid chromatography–tandemmass spectrometry, PloS one 9 (4) (2014) e93805.[11] S. Marino, I. B. Hogue, C. J. Ray, D. E. Kirschner, A methodology for performing global uncertainty and sensitivityanalysis in systems biology, Journal of theoretical biology 254 (1) (2008) 178–196.[12] W. L. Miller, R. J. Auchus, The molecular biology, biochemistry, and physiology of human steroidogenesis and its disorders,Endocrine reviews 32 (1) (2011) 81–151.[13] S. M. Blower, H. Dowlatabadi, Sensitivity and uncertainty analysis of complex models of disease transmission: an HIVmodel, as an example, International Statistical Review/Revue Internationale de Statistique (1994) 229–243.[14] V. Rohatgi, A. Saleh, Wiley series in probability and statistics, Hoboken: John Wiley & Sons, Inc.[15] L. Mora-L´opez, J. Mora, An adaptive algorithm for clustering cumulative probability distribution functions using theKolmogorov–Smirnov two-sample test, Expert Systems with Applications 42 (8) (2015) 4016–4021.[16] L. van der Maaten, G. Hinton, Visualizing data using t-sne, Journal of machine learning research 9 (Nov) (2008) 2579–2605.[17] A. Labhart, Clinical endocrinology: theory and practice, Springer Science & Business Media, 2012.[18] L. A. Harris, J. F. Selgrade, Modeling endocrine regulation of the menstrual cycle using delay differential equations,Mathematical biosciences 257 (2014) 11–22.
Appendix A. Original Graham-Selgrade Model Equations [5]
Releasable FSH: d
F SH ρ d t = v F c F,I S Λ K F,I + S Λ − k F c F,P P c F,E E F SH ρ (A.1)Serum FSH: d F SH d t = 1 V · k F c F,P P c F,E E F SH ρ − δ F F SH (A.2)Releasable LH: d LH ρ d t = (cid:20) v L TK L,T + T + v L E n K nmL + E n (cid:21) ·
11 + P K iL,P ( c L,T T ) − k L c L,P P c L,E E LH ρ (A.3)Serum LH: d LH d t = 1 V · k L c L,P P c L,E E LH ρ − δ L LH (A.4)Follicular phase: dΦd t = f · TT + f F SH (cid:16) h c Φ ,T T/T (cid:17) + F SH − f LH (cid:16) h c Φ ,F FSH (cid:17) + LH · Φ (A.5)Ovulatory phase: dΩd t = f LH (cid:16) h c Φ ,F FSH (cid:17) + LH · Φ − wS Ω (A.6)Luteal phase: dΛd t = wS Ω − l (1 − S )Λ (A.7)LH Support: d S d t = ˆ sLH m h ms + LH m · (1 − S ) − δ s S (A.8)Serum T: d T d t = t − δ T T + [ t G ( F + c T,F F ) + t G G F ] · (A.9) · (cid:20) Φ + τ Ω + τ S Λ + τ (cid:18) − Φ + Ω + ΛΨ (cid:19)(cid:21)
Intermediate T: d T γ d t = t g G G F − t g F SHh + F SH T γ (A.10) erum E : d E d t = e − δ E E + t g F SHh + F SH T γ · [Φ + η Λ S ] (A.11)Serum P : d P d t = − δ P P + pLHLH + h p · Λ S (A.12) Functional Forms. • Insulin-stimulated conditions ( α > ) G = G ( α ) G = G ( α ) D ( α ) = LH [ G + A ] + LH [ G B + A · ( B + C )] + A · B · CF ( LH, α ) = LH / D ( α ) F ( LH, α ) =
LH/ D ( α ) • Basal conditions ( α = 0 ) G = G = 1 κ = 1 + Aκ = B + A ( B + C ) κ = ABC D = κ LH + κ LH + κ F ( LH ) = LH / D F ( LH ) = LH/ D Appendix B. Reduced Model Parameters
Pituitary Parameters
Parameter Value v F Ki F I k F c F P c F E c F I v L v L Km L Ki LP k L c LP c LE Parameter Value f f h h w l s δ S η κ h s t γ e p Table B.2. Optimal parameters generated from fitting the reduced model to original model output. Listed here are theestimated parameters. Other fixed parameters appearing in the model remain unchanged from [5].
Appendix C. Derivation of Testosterone-Dependent Terms