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 ).

To model ovulation and follicle growth, we developed an equation for the size of follicle i, denoted x_i . We next describe the reasoning step by step. Since follicle cells replicate to make more of themselves, we begin with the usual d x i d t ∼   x i . 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

F S H ∼ 1 E ′ , 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 x_T is the total size of the follicles. We conclude that

Thus, the growth rate depends on the relative follicle size, x i   x T . 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 A_i , 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 A_i , and positive in between ( Figure 2B ).

To understand the dynamics of local androgen A_i , 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 β ( t ) = A T η x T ( t ) . 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 1 M 1 and 1 M 2 , as shown in Figure 2B M ₁ and M ₂ 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.

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 1 M 1 and 1 M 2 ( 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 M ₁ = 1 and M ₂ = 3.5.

A simulation with parameters in which M = 3 is chosen is shown in Figure 2D , namely M ₁ = 2.9 and M ₂ = 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 x i x T = 1 M . The other follicles have steady states

x i x T = 0 . 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 x i x T = 0 . 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 t^q . 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 t^q 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,

1 M 1 , 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 M ₁ and M ₂ 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 M ₁ < 3 and M ₂ > 3 (orange dot in Figure 3A ). If M>M ₂ or M<M ₁, the growth rate is negative (purple dots in Figure 3A ).

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 1 M 1 is to the right of the maximum of ϕ. That is, 1 M 1 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

1 M 1 > 1 M max or equivalently M < M_max . If we assume for simplicity that ϕ is a parabola, the maximum point is midway between the zeros, 1 M m a x = 1 2 ( 1 M 1 + 1 M 2 ) = 1 2 M 1 + M 2 M 1 M 2 , and thus

M m a x     = 2 M 1 M 2 M 1 + M 2 . Another criterion for stability for a general form of ϕ is ϕ ( 1 M ) > ϕ ( 0 ) ; however, since we assume M ₁,M ₂ > 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 M ₂ >> M ₁, 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 M ₁ can determine the ovulation number M. This parameter, 1 M 1 is proportional to the intra-follicle androgen concentration that is toxic to the follicle. For ovulation number M = 3, for example, M ₁ needs to be in the open interval between 1.5 and 3. For the human case M = 1, M ₁ needs to be lower than 1. Figure 3C shows the relation between the androgen toxicity parameter M ₁ and the ovulation number, when M ₂ is assumed to be very large, which gives the maximal range of possible ovulation numbers for a given value of M ₁. At M ₁ = 3.4, for example one can have 4, 5 or 6 ovulating follicles 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 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.

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, d x i d t = x i x T ϕ ( x i x T ) provides a constant velocity for a symmetric solution d x i d t = 1 M ϕ ( 1 M ) , 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, x_c it undergoes a developmental switch. Thereafter, the follicle obeys the present model’s competition dynamics, even if it shrinks below x_c ,

where A_ex is external androgen and x_T = Σx_i 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 M ₁and M ₂ to low numbers that prohibit ovulation. The simulations ( Figure 4B ) show follicles that grow to size x_c and then shrink, similar to the experimental ultrasound data for women with PCOS ( Figure 4A ).

To simulate ovulatory conditions, we use M₁=0.9, M₂ = 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.

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 x_T exceeded a threshold. Based on experiments (38), we set the threshold at 4.6 times the initial x_T (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 t_i throughout the follicular phase, and started follicles at t = 0 with initial sizes x _i =x _c e ^−γt_i so that they reach x_c at times t_i . Parameters were ′γ = 1 and x_c = 1; times and sizes can be scaled by appropriate parameter changes. External androgen A_ex shifts the biphasic function to the left. The larger M2, the less sensitive the dominant follicle dynamics to the entry of new follicles.

Lacker’s model in its most commonly used form (4) is, upto parameters which can be rescaled away, dx_i /dt = x_i g(x_i ,x_t ) with g(x_i ,x_t ) = 1 – (x_T – M ₁ x_i ) (x_T – M ₂ x_i ) and x_T = Σ x_i . Symmetric solutions are obtained when M follicles have equal sizes x__i = x__T /M and the rest have x_i =0. The symmetric solution obeys d x T d t = x T + μ x T 3 with μ = − ( 1 − M 1 M ) ( 1 − M 2 M ) . 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 x i = 1 | μ | .

The present model (Eq 1-2) can be obtained as follows. First delete the 1 from the function g, resulting in d x i d t = x i x T 2 ( 1 − M 1 x i x T ) ( M 2 x i x T − 1 ) . This is the present model times x T 3 . Time can be transformed to a new variable τ such that d τ / d t = x T 3 , so that d x i d τ = u i ϕ ( u i ) with u i = x i x T and ϕ(u_i ) = (1 – M ₁ u_i ) (M ₂ u_i - 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.