- Department of Molecular Cell Biology, Weizmann Institute of Science, Rehovot, Israel
Every menstrual cycle, many follicles begin to develop but only a specific number ovulate. This ovulation number determines how many offspring are produced per litter, and differs between species. The physiological mechanism that controls ovulation number is unknown; a class of mathematical models can explain it, but these models have no physiological basis. Here, we suggest a physiological mechanism for ovulation number control, which enables selection of a specific number of follicles out of many, and analyze it in a mathematical model of follicular growth. The mechanism is based on a signal, intra-follicular androgen concentration, that measures follicle size relative to the other follicles. This signal has a biphasic effect, suppressing follicles that are too large or too small compared to others. The ovulation number is determined by the androgen inhibitory thresholds. The model has a scaling symmetry that explains why the dominant follicles grow linearly with time, as observed in human ultrasound data. This approach also explains how chronic hyperandrogenism disrupts ovulation in polycystic ovary syndrome (PCOS), a leading cause of infertility. We propose specific experiments for testing the proposed mechanism.
Introduction
In mammals a large number of follicles start growing every menstrual cycle, but only M of the follicles ovulate, and the rest die in a process called atresia (Figure 1A) (2–6). The ovulation number M is species-specific. In humans and elephants M = 1 except for rare twin events. In young mice M ≈ 8 (7). The question of how M follicles are chosen is called the “choose M” problem (8).
Figure 1 Ovarian follicles compete for ovulation during the follicular phase of the menstrual cycle, under control of the Hypothalamus-Pituitary-Ovary (HPO) axis. (A) Follicles join the menstrual cycle and grow, M of them “win the race” and ovulate, and the rest die (atretic follicles). In this schematic figure, M = 1 is shown. (B) The HPO axis controls circulating hormone levels. H denotes hypothalamus, P pituitary, O ovaries (specifically the ovarian follicles), and E estradiol. Estradiol control of gonadotropin production changes sign from negative to positive at high sustained estradiol levels. (C) Follicles in the ovary produce steroid hormones. LH induces theca cells to convert cholesterol to androgen. Most of the androgen goes to the circulation, and a small amount is converted to estrogen by the granulosa cells under control of FSH. (D) The dominant human follicle size measured by ultrasound grows with an approximately constant velocity (linear growth) in the follicular phase, adapted from (1). Regression line for the first 13 days of the follicular phase is shown, after which the control of Estradiol changes sign.
The fitness of an organism depends on its ovulation number M. The higher M the more offspring, but M that is too high places a load on parental care, leaving fewer surviving offspring. The optimal M is at intermediate values, an observation known as Lack’s principle (9, 10). It is not clear how the ovulation number problem is solved physiologically: what determines M?
Ovulation is regulated by the endocrine Hypothalamus-Pituitary-Ovary (HPO) axis (Figure 1B) (11–14). In the HPO axis, the hypothalamus secretes gonadotropin-releasing hormone (GnRH), which promotes secretion of two gonadotropins from the pituitary, follicle stimulating hormone (FSH), and luteinizing hormone (LH). LH promotes production of androgen (A) in theca cells of the follicle. Androgen is converted to estradiol (E) and other estrogens by the granulosa cells of the follicle under control of FSH (Figure 1C). FSH also promotes follicular growth and survival.
The ovarian follicles implement a negative feedback loop, in which estradiol inhibits the production of upstream hormones in the HPO axis. At late stages of the follicular phase, the phase of the menstrual cycle in which follicles compete and in which M are chosen to ovulate, the feedback switches sign and becomes positive. High estradiol activates LH production, triggering the LH surge which causes ovulation. Ovulation is dysregulated in a prevalent disorder called polycystic ovaries syndrome (PCOS), linked with excessive androgen levels and impaired fertility (15–18). In cases of anovulatory PCOS, competing follicles stop growing prematurely and regress (17).
The mechanism that determines the ovulation number remains a mystery, and is a topic of current research (5, 6, 19, 20). The ovulation number is known to be regulated by circulating factors, since removing one ovary does not reduce the total number of eggs released during ovulation by half (7, 21, 22). An elegant theory of the ovulation number was developed by Lacker et al. (4, 23–25) in the 1980’s. Lacker’s model can provide a choice of specific M. In the model, circulating estradiol secreted by the follicles provides systemic control over ovulation number, with a biphasic effect in which follicles that are too small or too large are eliminated. The biphasic effect is central to the model, but, as acknowledged by Lacker, it lacks a physiological mechanism for this effect. Defining a physiological mechanism can advance our understanding, and offer experimental tests and therapeutic points of intervention.
Lacker’s model is also inconsistent with the dynamics of follicle growth. It shows super-exponential growth of the dominant follicles, reaching infinite size at a finite time. This is in contrast to more recent measurements of follicular growth profiles in women, not available at the time the model was formulated, that show approximately linear growth with time of the dominant follicle (Figures 1D and S1) (1, 26). Lacker’s model and later variants (27, 28) are also not consistent with follicle dynamics in PCOS (17), because they show follicles that are growth arrested and persist at an intermediate size, rather than follicles that grow and then shrink.
It would be important to develop a model for ovulation number control based on physiological mechanisms, which can explain the ‘choose M′ problem, linear follicle growth and the origin of conditions such as PCOS.
In this paper, we combine multiple lines of evidence to propose a physiological mechanism for ovulation number control, and to develop a minimal mathematical model for follicular growth. The main regulators of growth in the model are systemic FSH that enhances follicle growth, and local androgen in each follicle which has a biphasic effect on follicle growth. The model explains linear follicle growth based on an invariance property. It provides a mechanism for how high androgen levels cause PCOS and its anovulatory dynamics, and a framework to understand the frequency of dizygotic twin ovulations.
Results
Biphasic Model for Ovulation Control Based on Local Androgen
We sought a mechanism by which follicles that are too large or too small compared to the other follicles can be removed. Such a principle was assumed in Lacker’s mathematical model, but with no physiological underpinning. We suggest a candidate for such a biphasic control based on the local androgen concentration in each follicle, as we describe next.
Androgens are produced by follicles in the ovary and, to a lesser extent, by the adrenal glands. The androgens androstenedione and testosterone are produced in theca cells of the follicles. Small amounts of these androgens are used as the precursors for estradiol production by the adjacent granulosa cells of the follicle (Figure 1C). In sheep, humans and primates the main precursor for estradiol is androstenedione, whereas in rodents the main precursor is testosterone (29). Most of the androgens are released to the circulation, where they are diluted to low levels, while the androgen inside each follicle is at a much higher concentration.
Androgen is involved in follicular growth, as evident both by controlled experiments and disease states. The androgen receptor is expressed in granulosa cells (30). At low to physiological levels, androgen stimulates follicle growth and survival in vitro and in vivo (30). However, at high levels, androgen inhibits or interferes with follicular growth (31). This is seen in states with high androgen, such as polycystic ovary syndrome (PCOS), congenital adrenal hyperplasia (16, 32), hyperandrogenism in female-to-male transsexuals (33, 34), and aromatase deficiency (35). These states are related to a polycystic ovary morphology, in which ovarian follicles stop growing prematurely and ovulation is prevented (16). Excess external androgen can also induce PCOS symptoms in animal models (36).
These seemingly paradoxical effects of androgen led us to posit that androgen has a biphasic (inverse-U-shaped) effect on follicle growth rate (Figures 2A, B).
Figure 2 Model for ovulation number control based on a biphasic effect of local androgen. (A) Growth of individual follicles xi in the model is controlled by systemic feedback (circulating FSH and Estradiol, E) and local feedback (androgen concentration at each follicle, Ai). Follicles secrete estradiol and contribute to the total circulating estradiol, E, which inhibits FSH. Both local (intra-follicular) androgen and FSH control follicular growth, where FSH promotes it and androgen has a biphasic effect. (B) Biphasic control of follicular growth by androgen, which is proportional to the relative follicle size , is described by the function ϕ. Its zero-crossing points are and . (C) Simulation of the model with five initial follicles and M1 = 0.9, M2 = 3.5, in which a single follicle ovulates. (D) Simulation with five initial follicles, M1 = 2.9, M2 = 7.5, in which three follicles ovulate. (E) Simulation with M1 = 0.9, M2 = 10 and 15 initial follicles in which one follicle ovulates. (F) Simulation with M1 = 5.9, M2 = 13 and 15 initial follicles in which six follicles ovulate. Note that units of time can be rescaled by multiplying ϕ by a velocity constant (α), and were set so the simulation lasts 14 days by a suitable constant in each panel.
To model ovulation and follicle growth, we developed an equation for the size of follicle i, denoted xi. We next describe the reasoning step by step. Since follicle cells replicate to make more of themselves, we begin with the usual . We next add the fact that FSH drives the growth of follicles
To understand FSH as a function of time, recall that it is inhibited in the HPO axis by estrogen E (Figure 2A). Thus, we model
, as qualitatively observed in serum hormone measurements (Figure S2). Estrogen is assumed to be produced by each follicle in proportion to its size. Thus, total serum estrogen goes as the sum of follicle sizes:
Where xT is the total size of the follicles. We conclude that
Thus, the growth rate depends on the relative follicle size, . If this was the only control on follicular growth, each follicle would keep growing and never be removed. We next reasoned that since each follicle makes androgen, with a high intra-follicular concentration Ai, follicle growth rate also depends on the biphasic effect of androgen. Thus
Where ϕ is a biphasic function: it is negative at low and high levels of Ai, and positive in between (Figure 2B).
To understand the dynamics of local androgen Ai, we assume that each follicle produces androgen in proportion to its size,
where β(t) describes the time-dependent hormonal control by FSH and LH. Circulating androgen, secreted by the sum of all follicles, is observed to be nearly constant across the follicular phase except for the 2-3 days before ovulation (SI section 2, Figure S3). This adds a constraint to the model
where η is the fraction of androgen that leaves the ovary into the circulation and is assumed to be constant. The conclusion is that . The local androgen level is therefore, again, a function of the relative size of the follicle:
We end up with an equation in which the growth rate of a follicle depends in a biphasic manner on its size relative to the sum of all the other follicles:
We define the two zero crossing points of ϕ as and , as shown in Figure 2B M1 and M2 are parameters that will become important soon. For simplicity, in the simulations below we assume a parabolic form for ϕ , namely
Other biphasic functions of the relative follicle size lead to the same qualitative conclusions. The model is similar to Lacker's model with certain alterations that change its behavior. We relate this model to Lacker’s model in Methods.
The Model Solves ‘Choose M’ and Provides Linear Follicle Growth
This model can choose M follicles out of many, and provides growth at a constant velocity to the dominant follicle(s). Examples from simulations are shown in Figure 2 and Figure S4. The ovulation number is determined by the biphasic function ϕ, and in particular by its zero-crossing points and (Figure 2B), as shown below.
The simulations begin with follicles with random initial sizes. A simulation where a single follicle is chosen, M = 1, is shown in Figure 2C. The dominant follicle grows and the rest shrink. In this simulation the zero-crossing parameters are M1 = 1 and M2 = 3.5.
A simulation with parameters in which M = 3 is chosen is shown in Figure 2D, namely M1 = 2.9 and M2 = 7.5. Three dominant follicles emerge and converge onto identical growth trajectories while the rest shrink. Figures 2E, F show simulations with a larger number of initial follicles, and different parameters, in which M = 1 or M = 6 dominant follicles arise. Note that after an initial transient, the dominant follicles approach a constant velocity in all cases.
A mathematical analysis of the model is provided in the SI (SI section 3). An intuitive way to see the main features is to solve equation (1) for the case of M equal-sized follicles that grow while all others die (4). Such solutions are called symmetric solutions. The relative size of the Mgrowing follicles is . The other follicles have steady states
. For example, if we are interested in a symmetric solution with M = 3, each of the follicles is ⅓ of the total summed size. The other follicles have . Plugging this into Eq (1), we find that the growing symmetric follicles have a constant velocity
Thus
We note that a slightly more complicated model shown in the SI can provide polynomial growth to the follicle mass, as tq. This can resolve the relation between follicle diameter and mass. Cells in the follicle lie mainly on its surface in a layer that changes thickness with time, resulting in a relation between diameter and mass that lies somewhere between linear and quadratic (37) (SI section 4, Figure S8), so that follicle mass grows as tq where 1 < q < 2 when diameter grows linearly with time.
For the M follicles in the symmetric solution to grow in size, rather than shrink, the velocity must be positive. Thus ϕ must be positive. This requires that the relative size of the follicles,
, falls between the two zero points of ϕ so it can be in the region where the growth rate is positive. The condition for a growing solution is therefore:
Since M1 and M2 are parameters determined by the thresholds of local androgen for growth or death, they are assumed to be the same for all follicles, and depend on factors like androgen-receptor affinity. To obtain M = 3, for example, we need M1 < 3 and M2 > 3 (orange dot in Figure 3A). If M>M2 or M<M1, the growth rate is negative (purple dots in Figure 3A).
Figure 3 Ovulation number stability, twin ovulations, and the effect of excessive circulating androgen. (A) A stable solution of ovulating follicles should satisfy M1 < M < Mmax, thus is in the positive declining part of ϕ. (B) An intuitive but imprecise explanation for the criteria for the stable solution is that in the declining part of ϕ, deviations from the symmetric solution converge back to it, and thus the solution is stable. (C) Ovulation numbers as a function of M1 in the limit of M2 >> M1. When ovulation of one follicle is a possible ovulation number, it is the only stable solution. (D) A twin ovulation can occur even when M = 1 is the only stable solution, when two follicles start with a similar size and ovulate together near an unstable solution when high enough estrogen triggers the LH surge. Parameters: M1 = 0.5, M2 = 4, ovulation occurs when xT(t) = 4.6 xT(0) (38)(SI section 2). (E) Same simulation as in D continued to later times (without allowing ovulation), showing that the smaller dominant follicle eventually shrinks. (F) Excessive or exogenous androgen, Aex, shifts the parabola to the left, effectively changing the values of M1 and M2 and potentially disrupting ovulation.
In order for the symmetric solution to be stable, there is an additional condition. If we make one of the M follicles slightly bigger than the rest, we want that follicle to shrink back to be equal to the rest. To see this graphically we can note that the region of stability occurs when is to the right of the maximum of ϕ. That is, must lie in the declining phase of the biphasic function.
To see why, imagine that we make one of the M follicles slightly larger than the others (Figure 3B purple dot). If ϕ is declining, the follicle grows slightly slower than the rest, and thus shrinks in relative size. The solution is stable. Likewise, if one follicle is slightly smaller (Figure 3B yellow dot), it grows faster, and catches up, returning to the symmetric solution.
In contrast, in the rising part of ϕ, to the left of its maximum, the symmetric solution is unstable. A slightly larger follicle grows faster than the rest and keeps growing, breaking the symmetric solution.
Thus, the stability criterion is that
or equivalently M < Mmax. If we assume for simplicity that ϕ is a parabola, the maximum point is midway between the zeros, , and thus
. Another criterion for stability for a general form of ϕ is ; however, since we assume M1,M2 > 0 this criterion is redundant to the positive growth criterion. The criterion for stability and positive growth together is
If we assume for simplicity that M2 >> M1, we find
The SI shows that there are other fixed-point solutions for relative follicle sizes, but the symmetric solutions are the only stable ones. It also shows that a biphasic form of ϕ is essential for stable solutions to the choose-M problem (SI section 3).
This model explains how setting a physiological parameter like M1 can determine the ovulation number M. This parameter, is proportional to the intra-follicle androgen concentration that is toxic to the follicle. For ovulation number M = 3, for example, M1 needs to be in the open interval between 1.5 and 3. For the human case M = 1, M1 needs to be lower than 1. Figure 3C shows the relation between the androgen toxicity parameter M1 and the ovulation number, when M2 is assumed to be very large, which gives the maximal range of possible ovulation numbers for a given value of M1. At M1 = 3.4, for example one can have 4, 5 or 6 ovulating follicles in the model.
Twin Ovulations in the Model
The present model provides insight into the occurrence of twin ovulations. Such dizygotic (non-identical) twins occur in unassisted human ovulations with a frequency of about 1 out of 90 pregnancies (39, 40). Monozygotic twins, which originate after ovulation, are beyond the scope of this model.
In the model, dizygotic ovulations can occur even if the only stable solution is M = 1. If two follicles happen to start the race with very similar initial sizes, they grow together (Figure 3D). Given enough time, the smaller follicle would eventually die in the model because M = 2 is unstable (Figure 3E). But given the limited duration of the follicular phase due to the estradiol levels that trigger the LH surge, it sometimes happens that the second follicle is so close in size to the dominant follicle that it also makes it to the LH surge. In this case both follicles ovulate (Figure 3D).
The probability for triplets and quadruplets in this picture drops exponentially because it is increasingly unlikely to have three or four dominant follicles with such nearly identical initial sizes that can keep together throughout the process. This is similar to the observed drop in frequencies: Naturally, dizygotic twins occur in about one in 90 pregnancies, triplets in about one in 8000 pregnancies, and quadruplets in about one in 400,000 pregnancies (39, 41).
The Model Can Explain Qualitative Follicle Dynamics in PCOS
The present mechanism also explains how excessive chronic circulating androgen might interfere with ovulation, as occurs in many cases of PCOS. High levels of external androgen effectively shift the biphasic curve to the left (Figure 3F). High enough levels prevent growing solutions and no ovulations can occur.
The model also relates to the dynamics of follicles in PCOS. These dynamics have an interesting history with respect to modelling. In the 1980’s it was believed that in PCOS follicles become arrested and persist at an intermediate size. Lacker’s model (4) and later variants (28) showed such a steady-state of growth arrest for certain parameters. Mathematically, this growth arrest solution can be seen from the symmetric solution of Lacker’s model (Methods).
The assumption of persistent growth-arrest in PCOS recently changed when longitudinal ultrasound in humans was reported by Jarrett et al. (17). PCOS involves more follicles than normal ovulation. These follicles grow, reach a typical size of about 7.5mm, and then decline (Figure 4A). There are no persistent growth-arrested follicles in PCOS according to the available data.
Figure 4 The present model can provide dynamics that are qualitatively similar to ultrasound experiments in women with and without PCOS. (A) Data from Jarrett et al. on follicle diameter in a participant with PCOS shows continual entry of follicles, which stop growth at a typical size and then shrink. (B) Dynamics in simulations in which follicles grow independently until they reach a critical size of xc = 1 and then compete according to the androgen model. Parameters do not allow ovulation, with M1 = 0.9, M2 = 10, external androgen Aex = 5 and α~0.005. (C) Data from Jarrett et al. on follicle dynamics in a control participant without PCOS, showing a single dominant follicle that ovulates and gradual shrinkage of the other follicles. The dashed line indicates the start of the follicular phase. (D) A single dominant follicle is seen in simulations as in B, but with parameters that allow ovulation, M1 = 0.9, M2 = 10, Aex = 0, and α = 1. In the simulations, follicle size is in arbitrary units.
Notably, the present model differs from Lacker’s model in having no growth arrested steady-states. Instead, the symmetric solutions show either growth or decline with constant velocity, due to the scaling property of the model. Mathematically, provides a constant velocity for a symmetric solution , either positive (ovulation) or negative (decline).
We asked to what extent does the present model provide dynamics similar to those observed experimentally in humans. Due to its simplicity, it is unlikely that the model can give quantitatively precise agreement with experiments, and instead we aimed at a qualitative agreement. To compare the model to simulations, we adjusted the model to allow follicles to enter the cycle at different times as experimentally observed (Figure 4A). To do so, we assume that follicles begin to grow independently and exponentially at random times from a small initial size
When a follicle reaches a critical size, xc it undergoes a developmental switch. Thereafter, the follicle obeys the present model’s competition dynamics, even if it shrinks below xc,
where Aex is external androgen and xT = Σxi is the sum over all follicles regardless of size. In simulations, we used a parabolic form for ϕ given by Eq 2.
We simulated PCOS by adding high androgen levels that effectively shift M1and M2 to low numbers that prohibit ovulation. The simulations (Figure 4B) show follicles that grow to size xc and then shrink, similar to the experimental ultrasound data for women with PCOS (Figure 4A).
To simulate ovulatory conditions, we use M1=0.9, M2 = 10 with no external androgen. We also set a smaller number of follicles to enter the cycle per unit time, as observed (Figure 4C). The model displays a dominant follicle, typically one that starts early in the follicular phase of the cycle (Figure 4D). The dominant follicle grows with a constant velocity. The subdominant follicles grow, reach the critical size, and then gradually shrink, qualitatively similar to the experimental data (Figures 4C, D).
We were not able to obtain a closer fit to the experimental data with other model parameters. This indicates that a more complex model is required if precise fits to experiments are the goal.
The experimental data of Jarrett et al. also shows follicles that grow and then shrink during the luteal phase of the cycle, in both PCOS and non-PCOS participants. The present model concerns the follicular phase. The luteal phase involves a different hormonal profile, including progesterone produced by the corpus luteum and altered levels of FSH, LH and estrogen, which may require additional modelling.
We conclude that the proposed androgen-based mechanism leads to a mathematical model that can explain some of the qualitative features of follicle dynamics in normal and PCOS-like situations.
Discussion
We presented a physiological mechanism and mathematical model for the control of ovulation number in mammals. In this mechanism a specific number of follicles are chosen to ovulate whereas others die by means of a circuit in which follicles measure their relative sizes. The circuit has two signals: systemic control by FSH together with a new follicle-centric androgen signal. This local androgen signal inhibits a follicle’s growth if the relative size of the follicle is too large or too small. The model provides species-specific ovulation numbers regulated by the local androgen toxicity threshold. The model explains the observed constant velocity of the dominant follicle size, and explains the disruption of ovulation and decline of non-ovulating follicles in PCOS caused by high circulating androgen levels.
The mechanism is based on a biphasic (inverse U-shaped) effect of local androgen. Both high and low levels inhibit follicular growth. A biphasic function of follicle size was assumed in Lacker’s classic model of ovulation without a physiological basis; here we suggest that androgen provides a physiological basis.
Biphasic or non-monotonic behavior occurs in other physiological settings including endocrine, immune and neuronal systems (42, 43). Biphasic control of growth was suggested by Karin et al. (42, 43) to act as a mutant-resistance mechanism for mutants that mis-sense a feedback signal. In the present model, the biphasic effect plays a dynamical role, eliminating follicles that are too large or too small relative to the others.
The present model has a scaling property in which follicle growth rate depends only on its relative size compared to the sum of all other follicle sizes. Follicle sizes converge to symmetric solutions where all dominant follicles have equal sizes. This symmetry and scaling provides a constant growth velocity to the dominant follicles. Such a constant velocity is observed in longitudinal ultrasound studies of ovulation (1, 17). This differs from Lacker’s model and its variants, which do not scale each follicle to the total follicular mass, and as a result predict that the dominant follicles grow super-exponentially and reach infinite size at a finite time.
The present model has another key difference from Lacker’s model and its variants (28) in terms of the follicular growth profiles in PCOS. The earlier studies assumed that follicles in PCOS arrest and persist at a fixed size, and their models provide such a fixed-point solution. However, recent ultrasound measurements (17) indicate that in PCOS follicles stop growing prematurely and then shrink, rather than persisting at an intermediate size. The present model does not allow a persistent growth-arrest solution. Instead, it can provide turnover of follicles that grow and then shrink, similar to the experimental trajectories. This turnover is seen by simulating a process in which small follicles grow independently at first and then start competing according to the androgen mechanism at a certain developmental time or critical size.
Growth and shrinkage of follicles was also addressed in a recent mathematical model by Lange et al. (44) of anovulatory waves in cows. In these waves, follicles grow and then shrink without ovulating. To fit this behavior, Lange et al. introduced a term in which the follicular death rate increased linearly with time; this term lacks a known physiological basis, as stated by the authors, who intended to provide a minimal model to fit data rather than a physiological mechanism. The present androgen mechanism may be able to describe such anovulatory processes of follicular growth and shrinkage.
The proposed mechanism can be extended in several ways. One can consider competition between neighboring follicles, such as paracrine interactions and competition over blood vessels. It can also be extended by considering follicular developmental stages. This can be done by using the models of Clement, Monniaux and colleagues (5, 6, 14, 19, 45–50), which provide a continuous and deterministic description of follicle development based on the hormonally regulated partition of granulosa cells into different cell states: proliferation, differentiation and apoptosis. This detailed framework can help to explore the present mechanism by making the transition probabilities between cell states depend on intra-follicular androgen levels in a biphasic manner. Other extensions include changes that occur with age (51) and maturation and ageing of follicles (27, 52), as well as hormone-induced changes in pituitary gonadotroph cell mass, in analogy to recent models of the HPA axis (53–55).
The present model assumes that all follicles have the same response parameters to hormones and the same intrinsic growth rates. This can be extended to the case of heterogeneous follicle parameters using the methods of Chávez-Ross et al., which analyzed heterogeneity in Lacker’s model (28). Noise can also be added to the model. These extensions can break the strict ordering of the follicle sizes throughout their growth.
The present model suggests several experimental tests that can refute it. It would be important to explore in detail whether intra-follicular androgen has a biphasic effect on follicle growth rate. Such experiments can in principle be done in vitro. One can provide androgen at different concentrations and measure the effects on follicular-cell growth and death rates. These experiments can thus map the function ϕ.
In vivo, one can measure intra-follicular androgen levels from different follicles and relate these levels to single-cell gene expression from the same follicles. The model predicts that androgen concentration in a given follicle will have a biphasic effect on growth, differentiation and artesia expression programs in the theca and granulosa cells.
Another possible experiment is measurement of the intra-follicular androgen levels alongside the follicle diameter and its granulosa cell number and mass. Our theory predicts that intra-follicular androgen will be proportional to the relative size of each follicle, rather than to the absolute size.
One can use these approaches to compare ϕ from species with different ovulation numbers. The model predicts that the androgen thresholds for growth and toxicity, which are described by the zero points of ϕ denoted M1 and M2, will vary between species to provide the species-specific range of ovulation numbers. The model can also be tested by performing detailed longitudinal measurements of follicle sizes in species with different ovulation numbers using ultrasound (42, 55). The simplicity of the model suggests a qualitative comparison to data, and probably precludes precise quantitative fits. Future improvements can aim to provide such quantitative fits.
This study provides a perspective on PCOS in which the disorder is a fragility (exposed by modern conditions) of an essential androgen-based mechanism for ovulation number control. The biphasic androgen mechanism ensures that the ovulation number is in a specific range, which can be important for fitness of each species. The price of such a control mechanism is that excessive androgen can disrupt ovulation. Conditions of chronic excess ovarian androgen were probably rare historically before the rise of insulin resistance and obesity. We hope that the present framework will help to improve our understanding of ovulation, PCOS and infertility.
Methods
Simulations of the Follicular Growth Model
Initial conditions were N follicles with sizes given by random numbers from a uniform distribution between 0.05 and 0.15. Ovulation was simulated by the time when total follicle size xT exceeded a threshold. Based on experiments (38), we set the threshold at 4.6 times the initial xT (SI section 2). We used the solver scipy.integrate.ode (scipy version 1.7.0) of python 3.8.11. Code is available at github.com/michalshilo6/OvulationNumberControl.
To simulate continual entry of follicles in Figure 4, we used Eq 3-4. We generated a set of random times ti throughout the follicular phase, and started follicles at t = 0 with initial sizes xi=xce−γti so that they reach xc at times ti. Parameters were ′γ = 1 and xc = 1; times and sizes can be scaled by appropriate parameter changes. External androgen Aex shifts the biphasic function to the left. The larger M2, the less sensitive the dominant follicle dynamics to the entry of new follicles.
Relation to Lacker’s Model
Lacker’s model in its most commonly used form (4) is, upto parameters which can be rescaled away, dxi/dt = xi g(xi,xt) with g(xi,xt) = 1 – (xT – M1xi) (xT – M2xi) and xT = Σ xi. Symmetric solutions are obtained when M follicles have equal sizes x_i = x_T/M and the rest have xi =0. The symmetric solution obeys with . In ovulatory conditions, /mu is positive, and hence the dominant follicles diverge to infinite size at finite time. In anovulatory conditions /mu is negative and follicles reach growth arrest at a steady-state size .
The present model (Eq 1-2) can be obtained as follows. First delete the 1 from the function g, resulting in . This is the present model times . Time can be transformed to a new variable τ such that , so that with and ϕ(ui) = (1 – M1 ui) (M2 ui - 1). This is the present model with a transformed timescale. The transformed timescale changes the asymptotic growing solutions from singular with finite-time divergence in Lacker’s model to a constant velocity in the present model.
The constant factor 1 in Lacker’s function g provides exponential growth to small follicles. It is also responsible for the growth-arrested steady-states in Lacker’s model by balancing the parabolic part of g. In the present model of Figure 4, exponential growth occurs for small follicles but stops when follicles cross a developmental threshold at a critical size. Thereafter, there is no term to balance the biphasic function ϕ and hence no growth-arrested solutions.
Data Availability Statement
The original contributions presented in the study are publicly available. This data can be found here: https://github.com/michalshilo6/OvulationNumberControl.
Author Contributions
Conceptualization: MS and UA; Data curation: MS and UA; Formal analysis: all authors; Investigation: all authors; Methodology: MS and UA; Software: MS and AM; Supervision: UA; Visualization: MS and AM; Writing- original draft: MS and UA; Writing- review & editing: all authors. All authors contributed to the article and approved the submitted version.
Funding
This work was supported by funding from Cancer Research UK (C19767/A27145), and by internal funding of the Weizmann institute of science.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Acknowledgments
We thank Navot Silberstein, Charles Peskin, Nava Dekel, Philippa Melamed and our group members for discussions.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fendo.2022.816967/full#supplementary-material
References
1. Baerwald AR, Adams GP, Pierson RA. A New Model for Ovarian Follicular Development During the Human Menstrual Cycle. Fert Steril (2003) 80(1):116–22. doi: 10.1016/S0015-0282(03)00544-2
2. Paulini F, Silva RC, de P. Rôlo JLJ, Lucci CM. Ultrastructural Changes in Oocytes During Folliculogenesis in Domestic Mammals. J Ovarian Res 7(102):2014. doi: 10.1186/s13048-014-0102-6
3. Strauss JF, Williams CJ. “Chapter 8 - Ovarian Life Cycle,”. In: Strauss JF, Barbieri RL, editors. Yen and Jaffe’s Reproductive Endocrinology (Eighth Edition). Philadelphia: Elsevier (2019). p. 167–205.e9. doi: 10.1016/B978-0-323-47912-7.00008-1
4. Lacker HM. Regulation of Ovulation Number in Mammals. A Follicle Interaction Law That Controls Maturation. Biophys J (1981) 35(2):433–54. doi: 10.1016/S0006-3495(81)84800-X
5. Clément F, Monniaux D. Multiscale Modelling of Ovarian Follicular Selection. Prog Biophys Mol Biol (2013) 113(3):398–408. doi: 10.1016/j.pbiomolbio.2012.12.005
6. Echenim N, Monniaux D, Sorine M, Clément F. Multi-Scale Modeling of the Follicle Selection Process in the Ovary. Math Biosci (2005) 198(1):57–79. doi: 10.1016/j.mbs.2005.05.003
7. Finn CA. Reproductive Capacity and Litter Size in Mice: Effect of Age and Environment. Reproduction 6(2):205–214, 1963. doi: 10.1530/jrf.0.0060205
8. “Mathematical Physiology - II: Systems Physiology | James Keener | Springer.” . Available at: https://www-springer-com.ezproxy.weizmann.ac.il/gp/book/9780387793870 (Accessed Jul. 11, 2021).
9. Lack DL. The Natural Regulation of Animal Numbers (1954). Oxford: Clarendon Press. Available at: http://archive.org/details/naturalregulatio0000lack (Accessed Jul. 11, 2021).
10. Davies NB, Krebs JR, West SA. An Introduction to Behavioural Ecology. The Atrium, Southern Gate, Chichester West Sussex, PO19 8SQ, UK:John Wiley & Sons (2012).
11. Hall JE. “Chapter 7 - Neuroendocrine Control of the Menstrual Cycle,”. In: Strauss JF, Barbieri RL, editors. Yen and Jaffe’s Reproductive Endocrinology (Eighth Edition). Philadelphia: Elsevier (2019). p. 149–166.e5. doi: 10.1016/B978-0-323-47912-7.00007-X
12. McCartney CR, Marshall JC. Neuroendocrinology of Reproduction. 1600 John F. Kennedy Blvd. Ste 1800 Philadelphia, PA 19103-2899:Elsevier (2018). doi: 10.1016/B978-0-323-47912-7.00001-9
13. Bulun SE. Physiology and Pathology of the Female Reproductive Axis. Philadelphia:Elsevier Inc. (2011). doi: 10.1016/B978-1-4377-0324-5.00017-1
14. Clément F, Crépieux P, Yvinec R, Monniaux D. “Mathematical Modeling Approaches of Cellular Endocrinology Within the Hypothalamo-Pituitary-Gonadal Axis,”. In: Molecular and Cellular Endocrinology, vol. 518. 1600 John F. Kennedy Blvd. Ste 1800 Philadelphia, PA 19103-2899:Elsevier Inc. (2020). p. 110877. doi: 10.1016/j.mce.2020.110877
15. Escobar-Morreale HF. Polycystic Ovary Syndrome: Definition, Aetiology, Diagnosis and Treatment. Nat Rev Endocrinol (2018) 14(5):270–84. doi: 10.1038/nrendo.2018.24
16. Chang RJ, Dumesic DA. Polycystic Ovary Syndrome and Hyperandrogenic States Vol. 2. 1600 John F. Kennedy Blvd. Ste 1800 Philadelphia, PA 19103-2899:Elsevier Inc. (2019). doi: 10.1016/B978-0-323-47912-7.00021-4
17. Jarrett BY, Vanden Brink H, Oldfield AL, Lujan ME. Ultrasound Characterization of Disordered Antral Follicle Development in Women With Polycystic Ovary Syndrome. J Clin Endocrinol Metab (2020) 105(11):e3847–61. doi: 10.1210/clinem/dgaa515
18. Azziz R, Carmina E, Chen Z, Dunaif A, Laven JSE, Legro RS, et al. Polycystic Ovary Syndrome. Nat Rev Dis Primers (2016) 2:16057. doi: 10.1038/nrdp.2016.57
19. Echenim N, Clément F, Sorine M. Multiscale Modeling of Follicular Ovulation as a Reachability Problem. Multi Model Simul 6(3):895–912, 2007. doi: 10.1137/060664495
20. ScaramuzziABM RJ, BairdC DT, CampbellD BK, DriancourtE M-A, DupontA J, FortuneF JE, et al. Regulation of Folliculogenesis and the Determination of Ovulation Rate in Ruminants. Reprod Fertil. Dev 23(3):444–467, 2011. doi: 10.1071/RD09161
21. Biggers J, FINN C, McLAREN A. Long-Term Reproductive Performance of Female Mice. II. Variation of Litter Size With Parity. J Reprod Fertil (1962) 3:313–30. doi: 10.1530/jrf.0.0030313
22. Greenwald GS. Quantitative Study of Follicular Development in the Ovary of the Intact or Unilaterally Ovariectomized Hamster. Reproduction (1961) 2(3):351–61. doi: 10.1530/jrf.0.0020351
23. Akin E, Lacker HM. Ovulation Control: The Right Number or Nothing. J Math Biol (1984) 20(2):113–32. doi: 10.1007/BF00285341
24. Lacker HM, Akin E. How do the Ovaries Count? Math Biosci (1988) 90(1–2):305–32. doi: 10.1016/0025-5564(88)90072-7
25. Michael Lacker H, Percus A. How do Ovarian Follicles Interact? A Many-Body Problem With Unusual Symmetry and Symmetry-Breaking Properties. J Stat Phys (1991) 63(5–6):1133–61. doi: 10.1007/BF01030003
26. Gougeon A. Dynamics of Follicular Growth in the Human: A Model From Preliminary Results. Hum Reprod (1986) 1(2):81–7. doi: 10.1093/oxfordjournals.humrep.a136365
27. Chávez-Ross A. Analysis of Lacker’s Model With an Ageing Factor for the Control of Ovulation and Polycystic Ovary Syndrome. J Theor Med (2001) 3(4):247–70. doi: 10.1080/10273660108833079
28. Chávez-Ross A, Franks S, Mason HD, Hardy K, Stark J. Modelling the Control of Ovulation and Polycystic Ovary Syndrome. J Math Biol (1997) 36(1):95–118. doi: 10.1007/s002850050092
29. Young JM, McNeilly AS. Theca: The Forgotten Cell of the Ovarian Follicle. Reproduction (2010) 140(4):489–504. doi: 10.1530/REP-10-0094
30. Gervásio CG, Bernuci MP, Silva-de-Sá MF, de S. Rosa-E-Silva ACJ. The Role of Androgen Hormones in Early Follicular Development. ISRN Obstet Gynecol (2014) 2014:1–11. doi: 10.1155/2014/818010
31. Andersen CY, Lossl K. Increased Intrafollicular Androgen Levels Affect Human Granulosa Cell Secretion of Anti-Müllerian Hormone and Inhibin-B. Fert Steril (2008) 89(6):1760–5. doi: 10.1016/j.fertnstert.2007.05.003
32. Carmina E, Rosato F, Jannì A, Rizzo M, Longo RA. Relative Prevalence of Different Androgen Excess Disorders in 950 Women Referred Because of Clinical Hyperandrogenism. J Clin Endocrinol Metab (2006) 91(1):2–6. doi: 10.1210/jc.2005-1457
33. Pache T. D., Chadha S., Gooren L. J. G., Hop W. C. J., Jaarsma K. W., Dommerholt H. B. R., et al. Ovarian Morphology in Long-Term Androgen-Treated Female to Male Transsexuals. A Human Model for the Study of Polycystic Ovarian Syndrome? Histopathology (1991) 19(5):445–52. doi: 10.1111/j.1365-2559.1991.tb00235.x
34. Spinder T, Spijkstra J. J., van den Tweel J. G., Burger C. W., van Kessel H., Hompes P. G. A., et al. The Effects of Long Term Testosterone Administration on Pulsatile Luteinizing Hormone Secretion and on Ovarian Histology in Eugonadal Female to Male Transsexual Subjects. J Clin Endocrinol Metab (1989) 69(1):151–7. doi: 10.1210/jcem-69-1-151
35. Bulun SE. Aromatase and Estrogen Receptor α Deficiency. Fert Steril (2014) 101(2):323–9. doi: 10.1016/j.fertnstert.2013.12.022
36. Walters KA. Role of Androgens in Normal and Pathological Ovarian Function. Reproduction (2015) 149(4):R193–218. doi: 10.1530/REP-14-0517
37. Mc Natty KP. Hormonal Correlates of Follicular Development in the Human Ovary. Aust J Biol Sci (1981) 34(3):249–68. doi: 10.1071/BI9810249
38. Mclachlan RI, Cohen NL, Dahl KD, Bremner WJ, Souls MR. Serum Inhibin Levels During the Periovulatory Interval in Normal Women: Relationships With Sex Steroids and Gonadotrophin Levels. Clin Endocrinol (1990) 32(1):39–48. doi: 10.1111/j.1365-2265.1990.tb03748.x
39. Fellman J, Eriksson AW. On the History of Hellin’s Law. Twin Res Hum Genet (2009) 12(2):183–90. doi: 10.1375/twin.12.2.183
40. Zeleny C. The Relative Numbers of Twins and Triplets. Science (1921) 53(1368):262–3. doi: 10.1126/science.53.1368.262
41. “The Relative Numbers of Twins and Triplets | Science.” . Available at: https://science-sciencemag-org.ezproxy.weizmann.ac.il/content/53/1368/262/tab-pdf (Accessed Jul. 12, 2021).
42. Alon U. An Introduction to Systems Biology: Design Principles of Biological Circuits, 2nd Edition. Boca Raton, Fla: Chapman and Hall/CRC (2019).
43. Biphasic Response as a Mechanism Against Mutant Takeover in Tissue Homeostasis Circuits . Available at: https://scholar.google.co.il/citations?view_op=view_citaion&hl=en&user=dg7qorcAAAAJ&citation_for_view=dg7qorcAAAAJ:qjMakFHDy7sC (Accessed Jul. 12, 2021).
44. Lange A, Schwieger R, Plöntzke J, Schäfer S, Röblitz S. Follicular Competition in Cows: The Selection of Dominant Follicles as a Synergistic Effect. J Math Biol (2019) 78(3):579–606. doi: 10.1007/s00285-018-1284-0
45. Clément F. Optimal Control of the Cell Dynamics in the Granulosa of Ovulatory Follicles. Math Biosci (1998) 152(2):123–42. doi: 10.1016/S0025-5564(98)10026-3
46. Clément F, Coron J-M, Shang P. Optimal Control of Cell Mass and Maturity in a Model of Follicular Ovulation. SIAM J Cont Optim. (2013) 51(2):824–47. doi: 10.1137/120862247
47. Aymard B, Clément F, Monniaux D, Postel M. Cell-Kinetics Based Calibration of a Multiscale Model of Structured Cell Populations in Ovarian Follicles. SIAM J Appl Math (2016) 76(4):1471–91. doi: 10.1137/15M1030327
48. Clément F, Gruet MA, Monget P, Terqui M, Jolivet E, Monniaux D. Growth Kinetics of the Granulosa Cell Population in Ovarian Follicles: An Approach by Mathematical Modelling. Cell Prolif (2003) 30(6):255–70. doi: 10.1046/j.1365-2184.1997.00090.x
49. Clément F. D., Monniaux D., Stark J., Hardy K., Thalabard J.C., Franks S., et al. Mathematical Model of FSH-Induced cAMP Production in Ovarian Follicles. Am J Physiology-Endocrinol Metab (2001) 281(1):E35–53. doi: 10.1152/ajpendo.2001.281.1.E35
50. Clément F, Monniaux D, Thalabard J-C, Claude D. Contribution of a Mathematical Modelling Approach to the Understanding of the Ovarian Function. Comptes Rendus Biologies (2002) 325(4):473–85. doi: 10.1016/S1631-0691(02)01457-9
51. Faddy MJ, -Gosden R-G. “A Mathematical Model of Follicle Dynamics in the Human Ovary,” (1995). Available at: https://academic.oup.com/humrep/article-abstract/10/4/770/611171.
52. Mariana JC, Corpet F, Chevalet C. Lacker’s Model: Control of Follicular Growth and Ovulation in Domestic Species. Acta Biotheor (1994) 42(4):245–62. doi: 10.1007/BF00707391
53. Karin O, Raz M, Alon U. An Opponent Process for Alcohol Addiction Based on Changes in Endocrine Gland Mass. iScience (2021) 24(3):102127. doi: 10.1016/j.isci.2021.102127
54. “A New Model for the HPA Axis Explains Dysregulation of Stress Hormones on the Timescale of Weeks, in: Molecular Systems Biology . Available at: https://www-embopress-org.ezproxy.weizmann.ac.il/doi/full/10.15252/msb.20209510 (Accessed Dec. 16, 2021).
55. “Hormone Seasonality in Medical Records Suggests Circannual Endocrine Circuits, in: PNAS . Available at: https://www-pnas-org.ezproxy.weizmann.ac.il/content/118/7/e2003926118.short (Accessed Jul. 13, 2021).
Keywords: mathematical model for ovulation, ovulation number control, polycystic ovary syndrome, biphasic control, biphasic control of androgen, hyperanderogenism, mathematical model for follicular growth, Lacker’s model
Citation: Shilo M, Mayo A and Alon U (2022) A Mechanism for Ovulation Number Control. Front. Endocrinol. 13:816967. doi: 10.3389/fendo.2022.816967
Received: 17 November 2021; Accepted: 26 May 2022;
Published: 14 July 2022.
Edited by:
Frédérique Clément, Inria Saclay - Île-de-France Research Centre, FranceReviewed by:
James Sneyd, The University of Auckland, New ZealandGheorghe Craciun, University of Wisconsin-Madison, United States
Copyright © 2022 Shilo, Mayo and Alon. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Uri Alon, urialonw@gmail.com