A delay differential equation model of follicle waves in women

ABSTRACT

This article presents a mathematical model for hormonal regulation of the menstrual cycle which predicts the occurrence of follicle waves in normally cycling women. Several follicles of ovulatory size that develop sequentially during one menstrual cycle are referred to as follicle waves. The model consists of 13 nonlinear, delay differential equations with 51 parameters. Model simulations exhibit a unique stable periodic cycle and this menstrual cycle accurately approximates blood levels of ovarian and pituitary hormones found in the biological literature. Numerical experiments illustrate that the number of follicle waves corresponds to the number of rises in pituitary follicle stimulating hormone. Modifications of the model equations result in simulations which predict the possibility of two ovulations at different times during the same menstrual cycle and, hence, the occurrence of dizygotic twins via a phenomenon referred to as superfecundation. Sensitive parameters are identified and bifurcations in model behaviour with respect to parameter changes are discussed. Studying follicle waves may be helpful for improving female fertility and for understanding some aspects of female reproductive ageing.

KEYWORDS:

• Estradiol
• parameter
• pituitary gland
• superfecundation
• bifurcation

AMS 2000 MATHEMATICS SUBJECT CLASSIFICATION:

• 92D25
• 39A11

1. Introduction

Deterministic, mathematical models consisting of systems of differential equations have been used to describe hormonal regulation of the human menstrual cycle, for example, [4, 5, 8, 15, 29, 30, 32, 33]. This physiological system may be considered dual control because hormones secreted by the hypothalamus and the pituitary glands affect the ovaries and the ovarian hormones affect the brain, see Figure 1. Gonadotropin releasing hormone produced by the hypothalamus promotes the synthesis in the pituitary of follicle stimulating hormone (FSH) and luteinizing hormone (LH), which control ovarian activity and ovulation, see [16, 42, 45]. In turn, the ovaries produce estradiol (E2), progesterone (P4), inhibin A (InhA), and inhibin B (InhB) which influence the synthesis and release of FSH and LH, see [18, 20, 39].

Figure 1. The menstrual cycle is controlled by hormones secreted by the hypothalamus and the pituitary and by the ovaries. The cycle begins with menstruation and then follicles develop during the follicular phase. Ovulation occurs followed by the luteal phase during which hormones are secreted for pregnancy. The solid arrows indicate the effects of the gonadotropin hormones (red arrows) on the ovaries and the effects of the ovarian hormones (blue arrows) on the brain during the phases of the cycle.

The menstrual cycle consists of the follicular phase, ovulation, and the luteal phase (Figure 1) with an average duration of 28 days [37]. During the follicular phase, FSH initiates the development of 6–12 follicles and a dominant follicle is selected to grow to maturity [26, 27]. These follicles produce E2 which primes the pituitary to secrete LH in large amounts. At mid-cycle, a rise and fall of LH over a period of five days causes the dominant follicle to release its ovum and this LH surge is necessary for ovulation. The dominant follicle becomes the corpus luteum, which produces E2 and P4 in preparation for pregnancy. If fertilization does not occur, the corpus luteum atrophies which signals the end of one cycle and the beginning of the next.

In mammals such as cows, horses, and sheep, it has been known [11] that multiple follicles of ovulatory size develop during a single reproductive cycle but only one releases an ovum. Because these follicles grow sequentially (e.g. Figure 2, Panel a), they are described as follicle waves [12]. The growth of two or three follicle waves of this sort is stimulated by separate waves of FSH. Because luteal P4 suppresses a LH surge, the ovulatory follicle arises from the cohort that matures after the corpus luteum has atrophied and P4 has diminished. Baerwald et al. [1, 2] have observed follicle waves in normally cycling women. In Baerwald et al. [3], two or three follicle waves were illustrated during one menstrual cycle. These waves are initiated by rises in FSH levels. Specifically, for a two-wave cycle, the first wave is produced by the early follicular phase rise in FSH and the second wave, by the sharp FSH spike concurrent with the LH surge at mid-cycle (Figure 2, Panel b). The follicle wave corresponding to the mid-cycle FSH spike develops during the luteal phase of the cycle but does not ovulate because of the absence of a luteal LH surge (suppressed by luteal P4) and possibly because the follicle does not reach ovulatory size. When the third wave appears, it develops in response to an unusual late luteal/early follicular rise in FSH (Figure 3, Panel e) as reported by Baerwald et al. [3].

Figure 2. Graphs are taken from Baerwald et al. [3], based on results from Baerwald et al. [1, 2]. Panel (a) depicts two follicle waves per cycle and the corpus luteum (yellow body). Panels (b) and (c) graph corresponding LH, FSH, E2, and P4 concentrations. Two-thirds of the women studied had cycles exhibiting two follicle waves per cycle and their results are shown. Note that two follicle waves per cycle corresponds to two rises in FSH. From Baerwald et al. [3] by permission of Oxford University Press. Human Reproduction Update is published on behalf of the European Society of Human Reproduction and Embryology (ESHRE).

Figure 3. Graphs are taken from Baerwald et al. [3], based on results from Baerwald et al. [1, 2]. Panel (d) depicts three follicle waves per cycle and the corpus luteum (yellow body). Panels (e) and (f) graph corresponding LH, FSH, E2, and P4 concentrations. One-third of the women studied had cycles exhibiting three follicle waves per cycle and their results are shown. Note that three follicle waves per cycle corresponds to three rises in FSH. From Baerwald et al. [3] by permission of Oxford University Press. Human Reproduction Update is published on behalf of the ESHRE .

In this study, we illustrate follicle waves in normally cycling women by modifying the mathematical model introduced in Harris' thesis [14] and Harris-Clark et al. [15], which consists of 13 nonlinear, delay differential equations. We alter this system by including a single FSH threshold function to stimulate the growth of the first stage of follicular development as suggested by Zeleznik [44]. With a FSH threshold, the follicular phase rise in FSH and the mid-cycle surge in FSH will both promote the development of dominant follicles and the possibility of dual ovulations. To prevent a second dominant follicle, we add follicle atresia which depends on LH and this suppresses the development of a dominant follicle stage during the luteal phase. In her thesis, Margolskee [21] was the first to experiment with FSH thresholds and atresia and to demonstrate the possibility of follicle waves in a mathematical model.

The parameters of our system of differential equations are estimated here using the data of McLachlan et al. [23]. These data resemble the Baerwald hormone profiles in Figure 2 but also include inhibin A which is needed in our model because of its effect on FSH synthesis (Section 2.1). Simulations of our model approximate well the pituitary and ovarian hormone profiles in McLachlan et al. [23]. Also, simulations of the different stages of follicle development exhibit a follicle wave resulting in a dominant follicle during the follicular phase and a second wave of nonovulatory follicles during the early luteal phase similar to those depicted in Figure 2, Panel (a). Numerical experiments are carried out to show that this model may also produce three follicle waves per cycle. In particular, we input a FSH profile similar to that in Panel (e) of Figure 3 to illustrate that our model can produce three follicle waves like those in Panel (d) of Figure 3.

An interesting biological question is what conditions permit an ovulation from the second follicle wave and, hence, the possibility of fertilization of two ova at different times during the same menstrual cycle. The occurrence of dizygotic twins [13] in such a case is referred to as superfecundation [17, 40]. Our model may be used to illustrate the possibility of superfecundation. By removing the atresia term and making some other parameter modifications, we show that the second follicle wave can result in a second dominant follicle which will elicit a second LH surge and another ovulation.

2. Model development

Many mathematical models for hormonal control of the menstrual cycle track blood levels of hormones produced by the brain and the ovaries during the follicular and luteal phases of the cycle using systems of differential equations. Our model is based on a nonlinear system of 13 delay, differential equations first presented in Harris-Clark et al. [15] which includes the five reproductive hormones LH, FSH, E2, P4, and InhA. Model parameters are estimated using data for 33 normally cycling women in McLachlan et al. [23] which provides average daily blood levels of the 5 hormones over a 31-day month. Because the ovaries respond to daily averages [26], our model lumps the effects of the hypothalamus and the pituitary together and we consider ovarian hormone regulation of the synthesis and release of FSH and LH on a time scale of days (Section 2.1). In addition, because the clearance of the ovarian hormones from the blood is rapid compared to the clearance of the pituitary hormones and compared to the time scale for ovarian development, we assume that the blood levels of the ovarian hormones are at quasi-steady state [19]. Hence, their concentrations are expressed as linear combinations of ovarian stages of the cycle (Section 2.2). These simplifications permit the use of a daily time scale to track hormonal variation during the cycle. Finally, similar to the approach of Harris-Clark et al. [15], we construct our model in three steps: a pituitary system, an ovarian system, and the full system.

2.1. The pituitary system

Four differential equations are used to represent the synthesis, release, and clearance of LH and FSH from the pituitary into the blood under control of the ovarian hormones (Figure 4). These equations were derived in Schlosser and Selgrade [34] where the state variables, ${\mathrm{RP}}_{\mathrm{LH}}$ and ${\mathrm{RP}}_{\mathrm{FSH}}$, represent the amounts of LH and FSH in the pituitary reserve pool and the state variables, LH and FSH, represent concentrations in the blood. Assuming that V is constant blood volume, the four equations take the following general form: (1) $\frac{\mathrm{d}}{\mathrm{d}t}{\mathrm{RP}}_{\mathrm{LH}}={\mathrm{synthesis}}_{\mathrm{LH}}-{\mathrm{release}}_{\mathrm{LH}},$(1) (2) $\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{LH}=\frac{1}{V}{\mathrm{release}}_{\mathrm{LH}}-{\mathrm{clearance}}_{\mathrm{LH}},$(2) (3) $\frac{\mathrm{d}}{\mathrm{d}t}{\mathrm{RP}}_{\mathrm{FSH}}={\mathrm{synthesis}}_{\mathrm{FSH}}-{\mathrm{release}}_{\mathrm{FSH}},$(3) (4) $\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{FSH}=\frac{1}{V}{\mathrm{release}}_{\mathrm{FSH}}-{\mathrm{clearance}}_{\mathrm{FSH}}.$(4)

Figure 4. The brain and the blood compartments for the pituitary model are illustrated. E2, P4, and InhA stimulate (+) and/or inhibit (−) the synthesis and/or release of the pituitary hormones. Clearance from the blood is linear.

Synthesis and release rates for LH and FSH are described as rational functions of ovarian hormones in which stimulatory effects appear in the numerators and inhibitory effects appear in the denominators. Clinical research has indicated that P4 stimulates the release of LH and FSH into circulation [7] and that low levels of E2 inhibit the release of LH and FSH [38]. Also, it has been shown [20, 42] that the high levels of E2 produced by the dominant follicle elicit the mid-cycle LH surge, so the numerator of the LH synthesis term contains a Hill function of E2 to reflect this effect on LH. The inhibins inhibit FSH synthesis [42]. These features are discussed in detail in [34] and [15] where the following system expands Equations (1(1) $\frac{\mathrm{d}}{\mathrm{d}t}{\mathrm{RP}}_{\mathrm{LH}}={\mathrm{synthesis}}_{\mathrm{LH}}-{\mathrm{release}}_{\mathrm{LH}},$(1) )–(4(4) $\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{FSH}=\frac{1}{V}{\mathrm{release}}_{\mathrm{FSH}}-{\mathrm{clearance}}_{\mathrm{FSH}}.$(4) ): (5) $\frac{\mathrm{d}}{\mathrm{d}t}{\mathrm{RP}}_{\mathrm{LH}}=\frac{v_{0\mathrm{LH}}+\frac{v_{1\mathrm{LH}}{E}_2(t-d_E{)}^a}{K{m}_{\mathrm{LH}}^a+E_2(t-d_E{)}^a}}{1+\frac{P_4(t-d_P)}{K{i}_{\mathrm{LH},P}}}-\frac{k_{\mathrm{LH}}(1+c_{\mathrm{LH},P}{P}_4){\mathrm{RP}}_{\mathrm{LH}}}{1+c_{\mathrm{LH},E}{E}_2},$(5) (6) $\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{LH}=\frac{1}{V}\frac{k_{\mathrm{LH}}(1+c_{\mathrm{LH},P}{P}_4){\mathrm{RP}}_{\mathrm{LH}}}{1+c_{\mathrm{LH},E}{E}_2}-a_{\mathrm{LH}}\mathrm{LH},$(6) (7) $\frac{\mathrm{d}}{\mathrm{d}t}{\mathrm{RP}}_{\mathrm{FSH}}=\frac{v_{\mathrm{FSH}}}{1+\frac{\mathrm{InhA}(t-d_{\mathrm{InhA}})}{K{i}_{\mathrm{FSH},\mathrm{InhA}}}}-\frac{k_{\mathrm{FSH}}(1+c_{\mathrm{FSH},P}{P}_4){\mathrm{RP}}_{\mathrm{FSH}}}{1+c_{\mathrm{FSH},E}{E}_2^2},$(7) (8) $\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{FSH}=\frac{1}{V}\frac{k_{\mathrm{FSH}}(1+c_{\mathrm{FSH},P}{P}_4){\mathrm{RP}}_{\mathrm{FSH}}}{1+c_{\mathrm{FSH},E}{E}_2^2}-a_{\mathrm{FSH}}\mathrm{FSH}.$(8)

The data of McLachlan et al. [23] for the ovarian hormones are used to obtained statistical fits for input functions $E_2\left(t\right)$, $P_4\left(t\right)$, and $\mathrm{InhA}\left(t\right)$ in Equations (5(5) $\frac{\mathrm{d}}{\mathrm{d}t}{\mathrm{RP}}_{\mathrm{LH}}=\frac{v_{0\mathrm{LH}}+\frac{v_{1\mathrm{LH}}{E}_2(t-d_E{)}^a}{K{m}_{\mathrm{LH}}^a+E_2(t-d_E{)}^a}}{1+\frac{P_4(t-d_P)}{K{i}_{\mathrm{LH},P}}}-\frac{k_{\mathrm{LH}}(1+c_{\mathrm{LH},P}{P}_4){\mathrm{RP}}_{\mathrm{LH}}}{1+c_{\mathrm{LH},E}{E}_2},$(5) )–(8(8) $\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{FSH}=\frac{1}{V}\frac{k_{\mathrm{FSH}}(1+c_{\mathrm{FSH},P}{P}_4){\mathrm{RP}}_{\mathrm{FSH}}}{1+c_{\mathrm{FSH},E}{E}_2^2}-a_{\mathrm{FSH}}\mathrm{FSH}.$(8) ). Thus, in this step of the modelling process, Equations (5(5) $\frac{\mathrm{d}}{\mathrm{d}t}{\mathrm{RP}}_{\mathrm{LH}}=\frac{v_{0\mathrm{LH}}+\frac{v_{1\mathrm{LH}}{E}_2(t-d_E{)}^a}{K{m}_{\mathrm{LH}}^a+E_2(t-d_E{)}^a}}{1+\frac{P_4(t-d_P)}{K{i}_{\mathrm{LH},P}}}-\frac{k_{\mathrm{LH}}(1+c_{\mathrm{LH},P}{P}_4){\mathrm{RP}}_{\mathrm{LH}}}{1+c_{\mathrm{LH},E}{E}_2},$(5) )–(8(8) $\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{FSH}=\frac{1}{V}\frac{k_{\mathrm{FSH}}(1+c_{\mathrm{FSH},P}{P}_4){\mathrm{RP}}_{\mathrm{FSH}}}{1+c_{\mathrm{FSH},E}{E}_2^2}-a_{\mathrm{FSH}}\mathrm{FSH}.$(8) ) are time-dependent ordinary differential equations which are linear in the state variables, ${\mathrm{RP}}_{\mathrm{LH}}$, LH, ${\mathrm{RP}}_{\mathrm{FSH}}$, and FSH. The parameters $d_E$, $d_P$, and $d_{\mathrm{InhA}}$ denote the period of time between when changes in E2, P4, and InhA blood levels occur and changes in the synthesis rates of LH and FSH result. These time delays are used only in the synthesis terms because we assume that hormone synthesis is a more complicated biological process than hormone release. Then, the data of McLachlan et al. [23] for LH and FSH are used to estimate the 16 unknown parameters of Equations (5(5) $\frac{\mathrm{d}}{\mathrm{d}t}{\mathrm{RP}}_{\mathrm{LH}}=\frac{v_{0\mathrm{LH}}+\frac{v_{1\mathrm{LH}}{E}_2(t-d_E{)}^a}{K{m}_{\mathrm{LH}}^a+E_2(t-d_E{)}^a}}{1+\frac{P_4(t-d_P)}{K{i}_{\mathrm{LH},P}}}-\frac{k_{\mathrm{LH}}(1+c_{\mathrm{LH},P}{P}_4){\mathrm{RP}}_{\mathrm{LH}}}{1+c_{\mathrm{LH},E}{E}_2},$(5) )–(8(8) $\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{FSH}=\frac{1}{V}\frac{k_{\mathrm{FSH}}(1+c_{\mathrm{FSH},P}{P}_4){\mathrm{RP}}_{\mathrm{FSH}}}{1+c_{\mathrm{FSH},E}{E}_2^2}-a_{\mathrm{FSH}}\mathrm{FSH}.$(8) ) as was done in Harris-Clark et al. [15].

2.2. The ovarian system

The ovarian model, first developed by Selgrade and Schlosser [35], describes nine stages in the monthly growth and development of the ovaries. FSH and LH promote the growth of follicular tissue within a stage and the transfer of follicular tissue from one stage to the next. The capacity to produce hormones at each stage of the cycle is assumed to be proportional to the mass of the ovarian follicles or corpus lutea at that stage. The nine stages are as follows:

• (RcF) Recruited Stage: FSH-induced recruitment from the preantral pool,

• (GrF) Growing Stage: FSH sustaining growth of cohort,

• (DomF) Dominant Stage: selection of a dominant follicle,

• (${\mathrm{OvF}}_1$) First Ovulatory Stage: follicle ovulation,

• (${\mathrm{OvF}}_2$) Second Ovulatory Stage: follicle luteinization, and

• (${\mathrm{Lut}}_{1-4}$) Luteal Stages: four stages of corpus luteum decay.

The increase in mass and transfer of mass for each stage is described by a differential equation which is linear in the state variables with coefficients depending on the pituitary hormones. A transfer term in one equation appears as a growth term in the next equation as follows: (9) $\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{RcF}=b\cdot \mathrm{FSH}+\frac{c_1{\mathrm{FSH}}^p}{K{m}_F^p+{\mathrm{FSH}}^p}\mathrm{RcF}-c_2{\mathrm{LH}}^{\alpha }\mathrm{RcF},$(9) (10) $\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{GrF}=c_2{\mathrm{LH}}^{\alpha }\mathrm{RcF}+(c_3{\mathrm{LH}}^{\beta }-c_4\mathrm{LH})\mathrm{GrF},$(10) (11) $\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{DomF}=c_4\mathrm{LH}\cdot \mathrm{GrF}-c_5{\mathrm{LH}}^{\gamma }\mathrm{DomF}-\frac{c_{\mathrm{atr}}\mathrm{DomF}}{1+(\mathrm{LH}/K{i}_{\mathrm{atr}}{)}^q},$(11) (12) $\frac{\mathrm{d}}{\mathrm{d}t}{\mathrm{OvF}}_1=c_5{\mathrm{LH}}^{\gamma }\mathrm{DomF}-d_1{\mathrm{OvF}}_1,$(12) (13) $\frac{\mathrm{d}}{\mathrm{d}t}{\mathrm{OvF}}_2=d_1{\mathrm{OvF}}_1-d_2{\mathrm{OvF}}_2,$(13) (14) $\frac{\mathrm{d}}{\mathrm{d}t}{\mathrm{Lut}}_1=d_2{\mathrm{OvF}}_2-k_1{\mathrm{Lut}}_1,$(14) (15) $\frac{\mathrm{d}}{\mathrm{d}t}{\mathrm{Lut}}_2=k_1{\mathrm{Lut}}_1-k_2{\mathrm{Lut}}_2,$(15) (16) $\frac{\mathrm{d}}{\mathrm{d}t}{\mathrm{Lut}}_3=k_2{\mathrm{Lut}}_2-k_3{\mathrm{Lut}}_3,$(16) (17) $\frac{\mathrm{d}}{\mathrm{d}t}{\mathrm{Lut}}_4=k_3{\mathrm{Lut}}_3-k_4{\mathrm{Lut}}_4.$(17)

Applying the quasi-steady state assumption as did Bogumil et al. [4, 5], we write ovarian hormone concentrations as linear combinations of the stages that secrete these hormones: (18) $E_2=e_0+e_1\mathrm{GrF}+e_2\mathrm{DomF}+e_3{\mathrm{Lut}}_4,$(18) (19) $P_4=p_0+p_1{\mathrm{Lut}}_3+p_2{\mathrm{Lut}}_4,$(19) (20) $InhA=h_0+h_1\mathrm{DomF}+h_2{\mathrm{Lut}}_2+h_3{\mathrm{Lut}}_3+h_4{\mathrm{Lut}}_4.$(20)

The ovarian system of previous models [15, 22, 35, 36] used first-order FSH terms to initiate follicular growth. However, Zeleznik [44] suggests that FSH levels must exceed a threshold to initiate follicular development. Such a threshold effect may be captured mathematically using the Hill function of FSH in Equation (9(9) $\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{RcF}=b\cdot \mathrm{FSH}+\frac{c_1{\mathrm{FSH}}^p}{K{m}_F^p+{\mathrm{FSH}}^p}\mathrm{RcF}-c_2{\mathrm{LH}}^{\alpha }\mathrm{RcF},$(9) ): $\frac{c_1{\mathrm{FSH}}^p}{K{m}_F^p+{\mathrm{FSH}}^p}.$

This function of FSH has an increasing sigmoidal shaped graph. The Hill coefficient p determines the steepness of the graph, the parameter $c_1$ is the limiting value of the function, and the half-saturation parameter $K{m}_F$ is the FSH value giving a function value of $0.5c_1$. Here we show that parameters may be found so that this single Hill function as it appears in the differential Equation (9(9) $\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{RcF}=b\cdot \mathrm{FSH}+\frac{c_1{\mathrm{FSH}}^p}{K{m}_F^p+{\mathrm{FSH}}^p}\mathrm{RcF}-c_2{\mathrm{LH}}^{\alpha }\mathrm{RcF},$(9) ) results in the appearance of two follicle waves.

For a normally cycling woman, only the wave occurring in the folliclular phase will result in an ovulating follicle. To prevent a second dominant follicle, we add follicle atresia in Equation (11(11) $\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{DomF}=c_4\mathrm{LH}\cdot \mathrm{GrF}-c_5{\mathrm{LH}}^{\gamma }\mathrm{DomF}-\frac{c_{\mathrm{atr}}\mathrm{DomF}}{1+(\mathrm{LH}/K{i}_{\mathrm{atr}}{)}^q},$(11) ) which depends inversely on LH: $\frac{c_{\mathrm{atr}}\mathrm{DomF}}{1+(\mathrm{LH}/K{i}_{\mathrm{atr}}{)}^q}.$

Because luteal P4 inhibits LH synthesis, LH is low during the luteal phase resulting in a larger atresia term. Hence, a luteal phase follicle wave does not produce a dominant follicle and no second ovulation occurs. In Section 5, we show that the removal of this atresia term does result in a second ovulation and the possibility of fraternal twins.

In Equations (9(9) $\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{RcF}=b\cdot \mathrm{FSH}+\frac{c_1{\mathrm{FSH}}^p}{K{m}_F^p+{\mathrm{FSH}}^p}\mathrm{RcF}-c_2{\mathrm{LH}}^{\alpha }\mathrm{RcF},$(9) )–(17(17) $\frac{\mathrm{d}}{\mathrm{d}t}{\mathrm{Lut}}_4=k_3{\mathrm{Lut}}_3-k_4{\mathrm{Lut}}_4.$(17) ), $\mathrm{FSH}\left(t\right)$ and $\mathrm{LH}\left(t\right)$ are explicit functions of time t obtained statistically from the pituitary hormone data of [23]. Then, the McLachlan data for E2, P4, and InhA are used to estimate the 32 parameters in Equations (9(9) $\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{RcF}=b\cdot \mathrm{FSH}+\frac{c_1{\mathrm{FSH}}^p}{K{m}_F^p+{\mathrm{FSH}}^p}\mathrm{RcF}-c_2{\mathrm{LH}}^{\alpha }\mathrm{RcF},$(9) )–(20(20) $InhA=h_0+h_1\mathrm{DomF}+h_2{\mathrm{Lut}}_2+h_3{\mathrm{Lut}}_3+h_4{\mathrm{Lut}}_4.$(20) ) as was done in Section 2.1.

2.3. The merged system

The third and final step of model construction is to merge the pituitary system and ovarian system together to create a single 13-dimensional system of differential Equations (5(5) $\frac{\mathrm{d}}{\mathrm{d}t}{\mathrm{RP}}_{\mathrm{LH}}=\frac{v_{0\mathrm{LH}}+\frac{v_{1\mathrm{LH}}{E}_2(t-d_E{)}^a}{K{m}_{\mathrm{LH}}^a+E_2(t-d_E{)}^a}}{1+\frac{P_4(t-d_P)}{K{i}_{\mathrm{LH},P}}}-\frac{k_{\mathrm{LH}}(1+c_{\mathrm{LH},P}{P}_4){\mathrm{RP}}_{\mathrm{LH}}}{1+c_{\mathrm{LH},E}{E}_2},$(5) )–(17(17) $\frac{\mathrm{d}}{\mathrm{d}t}{\mathrm{Lut}}_4=k_3{\mathrm{Lut}}_3-k_4{\mathrm{Lut}}_4.$(17) ) with three auxiliary equations  (18(18) $E_2=e_0+e_1\mathrm{GrF}+e_2\mathrm{DomF}+e_3{\mathrm{Lut}}_4,$(18) )–(20(20) $InhA=h_0+h_1\mathrm{DomF}+h_2{\mathrm{Lut}}_2+h_3{\mathrm{Lut}}_3+h_4{\mathrm{Lut}}_4.$(20) ). The nonlinear coefficient functions of the inputs in Sections 2.1 and 2.2 now become nonlinear functions of state variables so the system is highly nonlinear. The parameters $d_E$, $d_P$, and $d_{\mathrm{InhA}}$ are discrete time delays so Equations (5(5) $\frac{\mathrm{d}}{\mathrm{d}t}{\mathrm{RP}}_{\mathrm{LH}}=\frac{v_{0\mathrm{LH}}+\frac{v_{1\mathrm{LH}}{E}_2(t-d_E{)}^a}{K{m}_{\mathrm{LH}}^a+E_2(t-d_E{)}^a}}{1+\frac{P_4(t-d_P)}{K{i}_{\mathrm{LH},P}}}-\frac{k_{\mathrm{LH}}(1+c_{\mathrm{LH},P}{P}_4){\mathrm{RP}}_{\mathrm{LH}}}{1+c_{\mathrm{LH},E}{E}_2},$(5) )–(17(17) $\frac{\mathrm{d}}{\mathrm{d}t}{\mathrm{Lut}}_4=k_3{\mathrm{Lut}}_3-k_4{\mathrm{Lut}}_4.$(17) ) are delay differential equations. The 51 parameters of Equations (5(5) $\frac{\mathrm{d}}{\mathrm{d}t}{\mathrm{RP}}_{\mathrm{LH}}=\frac{v_{0\mathrm{LH}}+\frac{v_{1\mathrm{LH}}{E}_2(t-d_E{)}^a}{K{m}_{\mathrm{LH}}^a+E_2(t-d_E{)}^a}}{1+\frac{P_4(t-d_P)}{K{i}_{\mathrm{LH},P}}}-\frac{k_{\mathrm{LH}}(1+c_{\mathrm{LH},P}{P}_4){\mathrm{RP}}_{\mathrm{LH}}}{1+c_{\mathrm{LH},E}{E}_2},$(5) )–(20(20) $InhA=h_0+h_1\mathrm{DomF}+h_2{\mathrm{Lut}}_2+h_3{\mathrm{Lut}}_3+h_4{\mathrm{Lut}}_4.$(20) ) must be re-estimated using the data of McLachlan et al. [23] but the parameters found in Sections 2.1 and 2.2 are good starting values. Optimized parameters are obtained using least-squares and the Nelder–Mead Simplex method implemented in MATLAB [25] (see [28] for details). These parameters are listed in Table 2 .

3. Model simulations – two follicle waves

Numerical solutions to Equations (5(5) $\frac{\mathrm{d}}{\mathrm{d}t}{\mathrm{RP}}_{\mathrm{LH}}=\frac{v_{0\mathrm{LH}}+\frac{v_{1\mathrm{LH}}{E}_2(t-d_E{)}^a}{K{m}_{\mathrm{LH}}^a+E_2(t-d_E{)}^a}}{1+\frac{P_4(t-d_P)}{K{i}_{\mathrm{LH},P}}}-\frac{k_{\mathrm{LH}}(1+c_{\mathrm{LH},P}{P}_4){\mathrm{RP}}_{\mathrm{LH}}}{1+c_{\mathrm{LH},E}{E}_2},$(5) )–(20(20) $InhA=h_0+h_1\mathrm{DomF}+h_2{\mathrm{Lut}}_2+h_3{\mathrm{Lut}}_3+h_4{\mathrm{Lut}}_4.$(20) ) are computed using the delay equations solver DDE23 in MATLAB. Testing a variety of initial conditions indicates that solutions approach a unique stable solution of period approximately 31.5 days. Initial conditions for this stable cycle are listed in Table 3 and are used to plot the follicle stages for two complete cycles in Figure 5 and the four hormones LH, FSH, E2, and P4 against the data of [23] in Figures 6 and 7. The left panel of Figure 5 displays the first three ovarian follicle stages, the recruited stage RcF, the growing stage GrF, and the dominant stage DomF. RcF and GrF exhibit two follicle waves – a large follicular phase wave which becomes the dominant, ovulatory follicle, and a smaller luteal phase wave. The atresia term in Equation (11(11) $\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{DomF}=c_4\mathrm{LH}\cdot \mathrm{GrF}-c_5{\mathrm{LH}}^{\gamma }\mathrm{DomF}-\frac{c_{\mathrm{atr}}\mathrm{DomF}}{1+(\mathrm{LH}/K{i}_{\mathrm{atr}}{)}^q},$(11) ) eliminates a luteal phase dominant follicle. In fact, the atresia is so successful at reducing the dominant stage after mid-cycle that all subsequent stages, ${\mathrm{OvF}}_1$ through ${\mathrm{Lut}}_4$, are quite small (see the right panel in Figure 5).

Figure 5. All ovarian follicle stages for the merged model Equations (5(5) $\frac{\mathrm{d}}{\mathrm{d}t}{\mathrm{RP}}_{\mathrm{LH}}=\frac{v_{0\mathrm{LH}}+\frac{v_{1\mathrm{LH}}{E}_2(t-d_E{)}^a}{K{m}_{\mathrm{LH}}^a+E_2(t-d_E{)}^a}}{1+\frac{P_4(t-d_P)}{K{i}_{\mathrm{LH},P}}}-\frac{k_{\mathrm{LH}}(1+c_{\mathrm{LH},P}{P}_4){\mathrm{RP}}_{\mathrm{LH}}}{1+c_{\mathrm{LH},E}{E}_2},$(5) )–(20(20) $InhA=h_0+h_1\mathrm{DomF}+h_2{\mathrm{Lut}}_2+h_3{\mathrm{Lut}}_3+h_4{\mathrm{Lut}}_4.$(20) ) with parameters of Table 2 are shown for two cycles. The model maintains two follicle waves for the recruited RcF and growing GrF stages (left panel). The atresia term reduces the dominant stage DomF after mid-cycle and effectively results in small subsequent ovarian stages (right panel).

Figure 6. Ovarian hormones E2 and P4 for the merged model Equations (5(5) $\frac{\mathrm{d}}{\mathrm{d}t}{\mathrm{RP}}_{\mathrm{LH}}=\frac{v_{0\mathrm{LH}}+\frac{v_{1\mathrm{LH}}{E}_2(t-d_E{)}^a}{K{m}_{\mathrm{LH}}^a+E_2(t-d_E{)}^a}}{1+\frac{P_4(t-d_P)}{K{i}_{\mathrm{LH},P}}}-\frac{k_{\mathrm{LH}}(1+c_{\mathrm{LH},P}{P}_4){\mathrm{RP}}_{\mathrm{LH}}}{1+c_{\mathrm{LH},E}{E}_2},$(5) )–(20(20) $InhA=h_0+h_1\mathrm{DomF}+h_2{\mathrm{Lut}}_2+h_3{\mathrm{Lut}}_3+h_4{\mathrm{Lut}}_4.$(20) ) with parameters of Table 2 are plotted against the McLachlan data [23] for two cycles.

Figure 7. Pituitary hormones LH and FSH for the merged model Equations (5(5) $\frac{\mathrm{d}}{\mathrm{d}t}{\mathrm{RP}}_{\mathrm{LH}}=\frac{v_{0\mathrm{LH}}+\frac{v_{1\mathrm{LH}}{E}_2(t-d_E{)}^a}{K{m}_{\mathrm{LH}}^a+E_2(t-d_E{)}^a}}{1+\frac{P_4(t-d_P)}{K{i}_{\mathrm{LH},P}}}-\frac{k_{\mathrm{LH}}(1+c_{\mathrm{LH},P}{P}_4){\mathrm{RP}}_{\mathrm{LH}}}{1+c_{\mathrm{LH},E}{E}_2},$(5) )–(20(20) $InhA=h_0+h_1\mathrm{DomF}+h_2{\mathrm{Lut}}_2+h_3{\mathrm{Lut}}_3+h_4{\mathrm{Lut}}_4.$(20) ) with parameters of Table 2 are plotted against the McLachlan data [23] for two cycles.

From Equations (18(18) $E_2=e_0+e_1\mathrm{GrF}+e_2\mathrm{DomF}+e_3{\mathrm{Lut}}_4,$(18) )–(20(20) $InhA=h_0+h_1\mathrm{DomF}+h_2{\mathrm{Lut}}_2+h_3{\mathrm{Lut}}_3+h_4{\mathrm{Lut}}_4.$(20) ) for the ovarian hormones, note that the second follicle wave affects only E2 and this is because of the GrF term in Equation (18(18) $E_2=e_0+e_1\mathrm{GrF}+e_2\mathrm{DomF}+e_3{\mathrm{Lut}}_4,$(18) ). As a result, the luteal E2 peak is higher than the data and almost 80% of the follicular E2 peak (Figure 6, left panel). Because E2 promotes LH synthesis, Equation (5(5) $\frac{\mathrm{d}}{\mathrm{d}t}{\mathrm{RP}}_{\mathrm{LH}}=\frac{v_{0\mathrm{LH}}+\frac{v_{1\mathrm{LH}}{E}_2(t-d_E{)}^a}{K{m}_{\mathrm{LH}}^a+E_2(t-d_E{)}^a}}{1+\frac{P_4(t-d_P)}{K{i}_{\mathrm{LH},P}}}-\frac{k_{\mathrm{LH}}(1+c_{\mathrm{LH},P}{P}_4){\mathrm{RP}}_{\mathrm{LH}}}{1+c_{\mathrm{LH},E}{E}_2},$(5) ), there is additional LH in the luteal phase appearing as a small bump in the LH profile (Figure 7, left panel). P4, which depends on the luteal stages, does not reach the luteal data peak (Figure 6, right panel) because these stages are small.

4. Three follicle waves

Baerwald et al. [3] illustrated three follicle waves per cycle depicted in Figure 3, Panel (d). Each wave corresponds to a different peak in FSH as observed in Figure 3, Panel (e), where a late luteal rise in FSH is interrupted by an early follicular dip in the FSH graph followed by its normal follicular phase rise. Daily FSH data for such a FSH profile are not available for us to use in a full model simulation. However, in order to test if our model can generate three follicle waves from this type of FSH profile, we use the ovarian model component discussed and implemented in Section 2.2. An input function $\mathrm{LH}\left(t\right)$ is obtained from the data of [23] as done in Section 2.2 but a new input for $\mathrm{FSH}\left(t\right)$ is defined using negative exponential functions of t to position three appropriately spaced rises as follows: (21) $\mathrm{FSH}=82.5+190{\mathrm{e}}^{-\frac{(t-8{)}^2}{10}}+180{\mathrm{e}}^{-\frac{(t-0.3{)}^2}{0.4}}+400{\mathrm{e}}^{-\frac{(t-15.2{)}^2}{3}}-20{\mathrm{e}}^{-\frac{(t-24{)}^2}{15}}.$(21)

This FSH profile is scaled to be consistent with the McLachlan data and is plotted in Figure 8.

Figure 8. The FSH input function (21(21) $\mathrm{FSH}=82.5+190{\mathrm{e}}^{-\frac{(t-8{)}^2}{10}}+180{\mathrm{e}}^{-\frac{(t-0.3{)}^2}{0.4}}+400{\mathrm{e}}^{-\frac{(t-15.2{)}^2}{3}}-20{\mathrm{e}}^{-\frac{(t-24{)}^2}{15}}.$(21) ) has three distinct rises during one cycle and is used to produce three follicle waves.

These input functions for $\mathrm{LH}\left(t\right)$ and $\mathrm{FSH}\left(t\right)$ are inserted into the system of differential Equations (9(9) $\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{RcF}=b\cdot \mathrm{FSH}+\frac{c_1{\mathrm{FSH}}^p}{K{m}_F^p+{\mathrm{FSH}}^p}\mathrm{RcF}-c_2{\mathrm{LH}}^{\alpha }\mathrm{RcF},$(9) )–(17(17) $\frac{\mathrm{d}}{\mathrm{d}t}{\mathrm{Lut}}_4=k_3{\mathrm{Lut}}_3-k_4{\mathrm{Lut}}_4.$(17) ) with auxiliary Equations (18(18) $E_2=e_0+e_1\mathrm{GrF}+e_2\mathrm{DomF}+e_3{\mathrm{Lut}}_4,$(18) )–(20(20) $InhA=h_0+h_1\mathrm{DomF}+h_2{\mathrm{Lut}}_2+h_3{\mathrm{Lut}}_3+h_4{\mathrm{Lut}}_4.$(20) ) to produce three follicle waves per cycle (Figure 9). The data for E2, P4, and InhA from McLachlan et al. [23] are used to optimize the parameters (see [28] for details) and these parameters are listed in Table 4. Figure 10 plots model simulations of the ovarian hormones against the data of [23]. Notice that three waves occur for the first two ovarian stages RcF and GrF but the atresia term eliminates the first and third waves in the dominate stage (Figure 9, left panel). To prevent an early follicular rise in E2 (Figure 10), we reduce the parameter $e_1$ in Equation (18(18) $E_2=e_0+e_1\mathrm{GrF}+e_2\mathrm{DomF}+e_3{\mathrm{Lut}}_4,$(18) ) to one-half of its value for the two-wave model (compare Tables 2 and 4).

Figure 9. Follicle stages are plotted for Equations (9(9) $\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{RcF}=b\cdot \mathrm{FSH}+\frac{c_1{\mathrm{FSH}}^p}{K{m}_F^p+{\mathrm{FSH}}^p}\mathrm{RcF}-c_2{\mathrm{LH}}^{\alpha }\mathrm{RcF},$(9) )–(20(20) $InhA=h_0+h_1\mathrm{DomF}+h_2{\mathrm{Lut}}_2+h_3{\mathrm{Lut}}_3+h_4{\mathrm{Lut}}_4.$(20) ) where the FSH input (21(21) $\mathrm{FSH}=82.5+190{\mathrm{e}}^{-\frac{(t-8{)}^2}{10}}+180{\mathrm{e}}^{-\frac{(t-0.3{)}^2}{0.4}}+400{\mathrm{e}}^{-\frac{(t-15.2{)}^2}{3}}-20{\mathrm{e}}^{-\frac{(t-24{)}^2}{15}}.$(21) ) has three distinct rises. The parameters are given in Table 4. The left panel depicts the first three stages of follicle growth and isolates waves in the recruited RcF (−) and the growing GrF (−) stages. Note that DomF (−) exhibits only one wave. The right panel gives the complete set of nine follicle stages.

Figure 10. The data of McLachlan et al. [23] for the ovarian hormones are compared to simulations of Equations (9(9) $\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{RcF}=b\cdot \mathrm{FSH}+\frac{c_1{\mathrm{FSH}}^p}{K{m}_F^p+{\mathrm{FSH}}^p}\mathrm{RcF}-c_2{\mathrm{LH}}^{\alpha }\mathrm{RcF},$(9) )–(20(20) $InhA=h_0+h_1\mathrm{DomF}+h_2{\mathrm{Lut}}_2+h_3{\mathrm{Lut}}_3+h_4{\mathrm{Lut}}_4.$(20) ) with parameters of Table 4. The FSH input function (21(21) $\mathrm{FSH}=82.5+190{\mathrm{e}}^{-\frac{(t-8{)}^2}{10}}+180{\mathrm{e}}^{-\frac{(t-0.3{)}^2}{0.4}}+400{\mathrm{e}}^{-\frac{(t-15.2{)}^2}{3}}-20{\mathrm{e}}^{-\frac{(t-24{)}^2}{15}}.$(21) ) with three distinct rises is used.

5. Superfecundation

Dizygotic twinning may happen when multiple ovulations occur during a single or consecutive menstrual cycles and two oocytes are fertilized [13]. In recent years, the incidence of this has increased partially due to more women receiving fertility treatments and more older women becoming pregnant [5, 6, 41]. Though rare, dizygotic twins via superfecundation could be a direct consequence of multiple follicle waves per cycle. Superfecundation is the ovulation and fertilization of two oocytes at different times during the same menstrual cycle [17]. This may occur if two follicle waves result in two dominant, ovulatory follicles during the same cycle. Most cases of superfecundation are noticed only when the twins have different fathers because they look drastically different or by chance a paternity test is conducted. A recent child support case was deemed ‘groundbreaking’ because the paternity test on a woman's twins showed that they had different fathers and the father who was being sued only had to pay support for one of the twins [24].

Here we show that removing the atresia term from Equation (11(11) $\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{DomF}=c_4\mathrm{LH}\cdot \mathrm{GrF}-c_5{\mathrm{LH}}^{\gamma }\mathrm{DomF}-\frac{c_{\mathrm{atr}}\mathrm{DomF}}{1+(\mathrm{LH}/K{i}_{\mathrm{atr}}{)}^q},$(11) ) and changing some parameters cause a luteal phase dominant follicle which secretes enough E2 to elicit a second LH surge and a second ovulation. If both oocytes are fertilized, then dizygotic twins may result. First, setting the parameter $c_{\mathrm{atr}}$ to zero removes the atresia term in Equation (11(11) $\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{DomF}=c_4\mathrm{LH}\cdot \mathrm{GrF}-c_5{\mathrm{LH}}^{\gamma }\mathrm{DomF}-\frac{c_{\mathrm{atr}}\mathrm{DomF}}{1+(\mathrm{LH}/K{i}_{\mathrm{atr}}{)}^q},$(11) ). This permits the growth of a dominant follicle stage during the luteal phase (Figure 11) and the resulting secretion of large amounts of luteal E2 (Figure 12). E2 promotes the synthesis of LH but P4 inhibits LH (see the first term in Equation 5(5) $\frac{\mathrm{d}}{\mathrm{d}t}{\mathrm{RP}}_{\mathrm{LH}}=\frac{v_{0\mathrm{LH}}+\frac{v_{1\mathrm{LH}}{E}_2(t-d_E{)}^a}{K{m}_{\mathrm{LH}}^a+E_2(t-d_E{)}^a}}{1+\frac{P_4(t-d_P)}{K{i}_{\mathrm{LH},P}}}-\frac{k_{\mathrm{LH}}(1+c_{\mathrm{LH},P}{P}_4){\mathrm{RP}}_{\mathrm{LH}}}{1+c_{\mathrm{LH},E}{E}_2},$(5) ). To guarantee a luteal LH surge (Figure 12), we diminish luteal P4 via Equation (19(19) $P_4=p_0+p_1{\mathrm{Lut}}_3+p_2{\mathrm{Lut}}_4,$(19) ) by decreasing parameter $p_1$ slightly and parameter $p_2$ by a factor of 10 (compare Table 2 and Table 5). Other ad hoc parameter adjustments are made to insure two LH surges of comparable strength, for example, $v_{1\mathrm{LH}}$ is increased and $v_{0\mathrm{LH}}$ is decreased. Then, the data for FSH, E2, and InhA from [23] are used to optimize remaining parameters (see [28] for details).

Figure 11. The first three stages of ovarian development for the model of follicle wave superfecundation are plotted for two cycles. The parameters are given in Table 5. The luteal phase DomF (−) stage produces large amounts of luteal E2 which causes a second LH surge, Figure 12.

Figure 12. E2 and LH from the model of follicle wave superfecundation is compared to the McLachlan data [28] for two cycles. During one cycle, a second LH surge is present due to increased luteal E2 levels produced by a second dominant follicle. A second LH surge allows for ovulation of the second dominant follicle for potential fertilization resulting in superfecundation.

6. Sensitivity and bifurcation analyses

Sensitivity analysis measures changes in system behaviour with respect to changes in parameters [43]. The purpose is to identify the parameters that have the largest impact on the dynamical behaviour of the system. The parameters that greatly affect the model output when changed slightly are known as sensitive parameters, while the parameters that have little effect on the model output are known as insensitive parameters [31]. Local sensitivity analysis examines the effects that the system undergoes when each parameter is changed one at a time. Here model outputs of crucial importance are the follicular phase E2 peak because it affects the LH surge and the FSH profile because it determines the follicle waves. So Panza [27] computes a normalized sensitivity coefficient for each model parameter measuring changes in E2 for four days centred at the E2 peak (ranked in Table 1, column 2) and another coefficient for changes in FSH for one complete cycle (ranked in Table 1, column 3). The five most sensitive parameters are ranked in Table 1 with $v_{\mathrm{FSH}}$ being the most sensitive for both model outputs. Notice that the five most sensitive parameters are the same for both outputs. The parameters $v_{\mathrm{FSH}}$, $v_{0\mathrm{LH}}$, and $K{m}_{\mathrm{LH}}$ affect the synthesis of the pituitary hormones, which control ovarian development. $K{m}_F$ is the half-saturation for the Hill function which promotes the growth of the first follicle stage RcF and α affects the transfer of mass from RcF to GrF (Equation (9(9) $\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{RcF}=b\cdot \mathrm{FSH}+\frac{c_1{\mathrm{FSH}}^p}{K{m}_F^p+{\mathrm{FSH}}^p}\mathrm{RcF}-c_2{\mathrm{LH}}^{\alpha }\mathrm{RcF},$(9) )). These two are the stages which exhibit follicle wave behaviour.

Table 1. The five most sensitive parameters ranked in order for model Equations (5(5) $\frac{\mathrm{d}}{\mathrm{d}t}{\mathrm{RP}}_{\mathrm{LH}}=\frac{v_{0\mathrm{LH}}+\frac{v_{1\mathrm{LH}}{E}_2(t-d_E{)}^a}{K{m}_{\mathrm{LH}}^a+E_2(t-d_E{)}^a}}{1+\frac{P_4(t-d_P)}{K{i}_{\mathrm{LH},P}}}-\frac{k_{\mathrm{LH}}(1+c_{\mathrm{LH},P}{P}_4){\mathrm{RP}}_{\mathrm{LH}}}{1+c_{\mathrm{LH},E}{E}_2},$(5) )–(20(20) $InhA=h_0+h_1\mathrm{DomF}+h_2{\mathrm{Lut}}_2+h_3{\mathrm{Lut}}_3+h_4{\mathrm{Lut}}_4.$(20) ) with parameters from Table 2. The model output for the second column is E2 peak and for the third column is FSH cycle profile. The five most sensitive parameters are the same for both outputs.

Rank	E2 peak	FSH cycle
1	$v_{\mathrm{FSH}}$	$v_{\mathrm{FSH}}$
2	$K{m}_F$	$v_{0\mathrm{LH}}$
3	$v_{0\mathrm{LH}}$	$K{m}_F$
4	α	$K{m}_{\mathrm{LH}}$
5	$K{m}_{\mathrm{LH}}$	α

Numerical simulations of Equations (5(5) $\frac{\mathrm{d}}{\mathrm{d}t}{\mathrm{RP}}_{\mathrm{LH}}=\frac{v_{0\mathrm{LH}}+\frac{v_{1\mathrm{LH}}{E}_2(t-d_E{)}^a}{K{m}_{\mathrm{LH}}^a+E_2(t-d_E{)}^a}}{1+\frac{P_4(t-d_P)}{K{i}_{\mathrm{LH},P}}}-\frac{k_{\mathrm{LH}}(1+c_{\mathrm{LH},P}{P}_4){\mathrm{RP}}_{\mathrm{LH}}}{1+c_{\mathrm{LH},E}{E}_2},$(5) )–(20(20) $InhA=h_0+h_1\mathrm{DomF}+h_2{\mathrm{Lut}}_2+h_3{\mathrm{Lut}}_3+h_4{\mathrm{Lut}}_4.$(20) ) with parameters from Table 2, which were optimized using the data of [23], indicate that the stable cycle is unique (Figures 6 and 7). If the sensitive parameters in Table 1 are varied near their best-fit values, the unique stable cycle appears to persist. Using the bifurcation software for delay differential equations, DDEBIFTOOL [10] confirms this suspicion by displaying a simple continuation curve through the best-fit value and no bifurcations nearby (see Figure 13 for the parameter $K{m}_F$). Hence, a woman whose reproductive cycle is determined by this parameter set would cycle normally and small perturbations would not significantly disturb her cycle.

Figure 13. This bifurcation diagram plots the maximal LH value along a periodic solution of Equations (5(5) $\frac{\mathrm{d}}{\mathrm{d}t}{\mathrm{RP}}_{\mathrm{LH}}=\frac{v_{0\mathrm{LH}}+\frac{v_{1\mathrm{LH}}{E}_2(t-d_E{)}^a}{K{m}_{\mathrm{LH}}^a+E_2(t-d_E{)}^a}}{1+\frac{P_4(t-d_P)}{K{i}_{\mathrm{LH},P}}}-\frac{k_{\mathrm{LH}}(1+c_{\mathrm{LH},P}{P}_4){\mathrm{RP}}_{\mathrm{LH}}}{1+c_{\mathrm{LH},E}{E}_2},$(5) )–(20(20) $InhA=h_0+h_1\mathrm{DomF}+h_2{\mathrm{Lut}}_2+h_3{\mathrm{Lut}}_3+h_4{\mathrm{Lut}}_4.$(20) ) against $K{m}_F$ with the remaining parameters from Table 2. SN denotes a saddle-node bifurcation and the * indicates the position of the cycle for best-fit $K{m}_F=115.18$. The solid blue curve represents stable cycles and the dashed red curve, unstable cycles.

The superfecundation system, Equations (5(5) $\frac{\mathrm{d}}{\mathrm{d}t}{\mathrm{RP}}_{\mathrm{LH}}=\frac{v_{0\mathrm{LH}}+\frac{v_{1\mathrm{LH}}{E}_2(t-d_E{)}^a}{K{m}_{\mathrm{LH}}^a+E_2(t-d_E{)}^a}}{1+\frac{P_4(t-d_P)}{K{i}_{\mathrm{LH},P}}}-\frac{k_{\mathrm{LH}}(1+c_{\mathrm{LH},P}{P}_4){\mathrm{RP}}_{\mathrm{LH}}}{1+c_{\mathrm{LH},E}{E}_2},$(5) )–(20(20) $InhA=h_0+h_1\mathrm{DomF}+h_2{\mathrm{Lut}}_2+h_3{\mathrm{Lut}}_3+h_4{\mathrm{Lut}}_4.$(20) ) with parameters from Table 5, produces dominant follicles DomF in both the follicular and luteal phases of the cycle primarily because of the absence of atresia. Introducing atresia by permitting the atresia coefficient $c_{\mathrm{atr}}$ to be nonzero and by taking $q=2.36$ and $K{i}_{\mathrm{atr}}=4.69$ results in interesting dynamical behaviour. Increasing $c_{\mathrm{atr}}$ diminishes the dominant follicle and the LH surge during the follicular phase but not the luteal phase (Figure 14). For $c_{\mathrm{atr}}$ near 20, the follicular phase LH surge is clearly too small to cause ovulation so ovulation occurs only during the luteal phase and the possibility of superfecundation is lost. For a woman whose cycle is described by this situation, her cycle length is normal (here about 30 days) but ovulation occurs much later than mid-cycle. Observe the LH profile in Figure 15, left panel, where the LH surge occurs at day 19 instead of day 15.

Figure 14. DomF and LH for the superfecundation model for three values of $c_{\mathrm{atr}}$ show a diminishing dominant follicle and LH surge during the follicular phase. If $c_{\mathrm{atr}}\approx 20$, there is just one ovulation per cycle.

Figure 15. LH profiles for $c_{\mathrm{atr}}=22.6517$ and $c_{\mathrm{atr}}=22.7841$ illustrating that the period of the stable cycle doubles between these values, that is, a period-doubling bifurcation occurs.

As $c_{\mathrm{atr}}$ increases above 20, a cascade of period-doubling bifurcations occurs. At $c_{\mathrm{atr}}\approx 22.654$, the cycle of period 30 days becomes unstable and a stable cycle of about 60 days appears as observed by comparing the LH profiles for $c_{\mathrm{atr}}=22.6517$ and $c_{\mathrm{atr}}=22.7841$ (Figure 15). Then at $c_{\mathrm{atr}}\approx 22.785$, a second period-doubling bifurcation occurs and a stable cycle of approximately 120 days appears. A portion of the bifurcation diagram showing these two period-doubling bifurcations is drawn in Figure 16.

Figure 16. The maximal LH value along a periodic solution of Equations (5(5) $\frac{\mathrm{d}}{\mathrm{d}t}{\mathrm{RP}}_{\mathrm{LH}}=\frac{v_{0\mathrm{LH}}+\frac{v_{1\mathrm{LH}}{E}_2(t-d_E{)}^a}{K{m}_{\mathrm{LH}}^a+E_2(t-d_E{)}^a}}{1+\frac{P_4(t-d_P)}{K{i}_{\mathrm{LH},P}}}-\frac{k_{\mathrm{LH}}(1+c_{\mathrm{LH},P}{P}_4){\mathrm{RP}}_{\mathrm{LH}}}{1+c_{\mathrm{LH},E}{E}_2},$(5) )–(20(20) $InhA=h_0+h_1\mathrm{DomF}+h_2{\mathrm{Lut}}_2+h_3{\mathrm{Lut}}_3+h_4{\mathrm{Lut}}_4.$(20) ) is plotted against $c_{\mathrm{atr}}$ with $q=2.36$, $K{i}_{\mathrm{atr}}=4.69$, and the remaining parameters from Table 5. PD denotes period-doubling bifurcations which occur at $c_{\mathrm{atr}}\approx 22.654$ and $c_{\mathrm{atr}}\approx 22.785$. The solid blue curves represent stable cycles and the dashed red curves, unstable cycles.

The sequence of bifurcations as $c_{\mathrm{atr}}$ increases from 22.5 to 23 resembles the standard period-doubling route to chaos [9] with a chaotic cycle at $c_{\mathrm{atr}}=22.85$ and a period-6 window around $c_{\mathrm{atr}}=22.894$ (Figure 17). The left panel of Figure 17 plots a menstrual cycle which is not exactly periodic but there are 12 LH surges per year. The right panel of Figure 17 plots a stable cycle of period about 180 days which is roughly six times the period of the cycle at $c_{\mathrm{atr}}=22.6517$. Such a interval in a bifurcation diagram is referred to as a period-6 window. A woman whose reproductive cycle is represented by this parameter regime could experience unpredictable cycles. This indicates that the two models corresponding to the parameters of Table 2 or of Table 5 with atresia are quite different. The model corresponding to Table 2 depicts a normally cycling woman with hormone levels consistent with data in the literature [23] and exhibiting two follicle waves per cycle. The model corresponding to Table 5 with atresia added represents a woman who may experience superfecundation but also whose cycle is easily disturbed by small variations in follicle atresia.

Figure 17. The left panel illustrates the LH profile of a stable, chaotic cycle at $c_{\mathrm{atr}}=22.85$. There are 12 LH surges per year but they do not have the same magnitude. The right panel plots the LH profile of a stable cycle of period approximately 180 days at $c_{\mathrm{atr}}=22.894$. This indicates the presence of a period-6 window in the $c_{\mathrm{atr}}$ bifurcation diagram.

7. Summary and discussion

This study describes a system of 13 nonlinear, delay differential equations with 51 parameters which models the phenomenon of follicle waves in normally cycling women first observed by Baerwald et al. [1–3]. Using model Equations (5(5) $\frac{\mathrm{d}}{\mathrm{d}t}{\mathrm{RP}}_{\mathrm{LH}}=\frac{v_{0\mathrm{LH}}+\frac{v_{1\mathrm{LH}}{E}_2(t-d_E{)}^a}{K{m}_{\mathrm{LH}}^a+E_2(t-d_E{)}^a}}{1+\frac{P_4(t-d_P)}{K{i}_{\mathrm{LH},P}}}-\frac{k_{\mathrm{LH}}(1+c_{\mathrm{LH},P}{P}_4){\mathrm{RP}}_{\mathrm{LH}}}{1+c_{\mathrm{LH},E}{E}_2},$(5) )–(20(20) $InhA=h_0+h_1\mathrm{DomF}+h_2{\mathrm{Lut}}_2+h_3{\mathrm{Lut}}_3+h_4{\mathrm{Lut}}_4.$(20) ), we exhibit menstrual cycles with two or three follicle waves corresponding to the number of rises in blood levels of FSH. Simulations with the parameters of Table 2 depict two follicle waves per cycle and resulting hormone profiles accurately approximate the data of [23] (Figures 6 and 7). Data with three distinctive rises in FSH are not available but we use the input function plotted in Figure 8 with three FSH hills timed to correspond to the hills of Panel (e) in Figure 3 from [3]. This input function in Equations (9(9) $\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{RcF}=b\cdot \mathrm{FSH}+\frac{c_1{\mathrm{FSH}}^p}{K{m}_F^p+{\mathrm{FSH}}^p}\mathrm{RcF}-c_2{\mathrm{LH}}^{\alpha }\mathrm{RcF},$(9) )–(20(20) $InhA=h_0+h_1\mathrm{DomF}+h_2{\mathrm{Lut}}_2+h_3{\mathrm{Lut}}_3+h_4{\mathrm{Lut}}_4.$(20) ) with the parameters of Table 4 causes the formation of three follicle waves as depicted in Figure 9 and the resulting ovarian hormones are plotted against the data of [23] in Figure 10. We conjecture that variations in InhB levels during the late luteal to early follicular transition might cause a FSH profile like Figure 3, Panel (e), but we do not have data to substantiate this conjecture. Since InhB inhibits FSH synthesis and is prominent during the late luteal to early follicular transition, decreases in InhB would result in increases in FSH.

The inclusion of an atresia term in Equation (11(11) $\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{DomF}=c_4\mathrm{LH}\cdot \mathrm{GrF}-c_5{\mathrm{LH}}^{\gamma }\mathrm{DomF}-\frac{c_{\mathrm{atr}}\mathrm{DomF}}{1+(\mathrm{LH}/K{i}_{\mathrm{atr}}{)}^q},$(11) ) prevents more than one ovulation per cycle. We show that removing this term and making some parameter adjustments permit two ovulations per cycle and, hence, the possibility of dizygotic twins via superfecundation. Simulations of Equations (5(5) $\frac{\mathrm{d}}{\mathrm{d}t}{\mathrm{RP}}_{\mathrm{LH}}=\frac{v_{0\mathrm{LH}}+\frac{v_{1\mathrm{LH}}{E}_2(t-d_E{)}^a}{K{m}_{\mathrm{LH}}^a+E_2(t-d_E{)}^a}}{1+\frac{P_4(t-d_P)}{K{i}_{\mathrm{LH},P}}}-\frac{k_{\mathrm{LH}}(1+c_{\mathrm{LH},P}{P}_4){\mathrm{RP}}_{\mathrm{LH}}}{1+c_{\mathrm{LH},E}{E}_2},$(5) )–(20(20) $InhA=h_0+h_1\mathrm{DomF}+h_2{\mathrm{Lut}}_2+h_3{\mathrm{Lut}}_3+h_4{\mathrm{Lut}}_4.$(20) ) with parameters of Table 5 exhibit two dominant follicles per cycle (Figure 11) and two LH surges (Figure 12, right panel). This situation results in two ovulations per cycle. Introducing atresia back into the superfecundation system by taking $q=2.36$ and $K{i}_{\mathrm{atr}}=4.69$ and increasing $c_{\mathrm{atr}}$ from zero, causes the loss of superfecundation and a cascade of period-doubling bifurcations. This system may produce erratic menstrual cycle behaviour. Hence, although model equations are the same, the two models corresponding to the parameters of Table 2 or of Table 5 with atresia have different dynamical behaviours and represent women with noticeably different reproductive hormone profiles.

Acknowledgements

The authors thank Dr Angela R. Baerwald of University of Saskatchewan for helpful communications.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported by National Science Foundation grant [DMS-1225607].

Appendix. Parameters and initial conditions

Table A.1. Parameters for the merged, two-wave model Equations (5(5) $\frac{\mathrm{d}}{\mathrm{d}t}{\mathrm{RP}}_{\mathrm{LH}}=\frac{v_{0\mathrm{LH}}+\frac{v_{1\mathrm{LH}}{E}_2(t-d_E{)}^a}{K{m}_{\mathrm{LH}}^a+E_2(t-d_E{)}^a}}{1+\frac{P_4(t-d_P)}{K{i}_{\mathrm{LH},P}}}-\frac{k_{\mathrm{LH}}(1+c_{\mathrm{LH},P}{P}_4){\mathrm{RP}}_{\mathrm{LH}}}{1+c_{\mathrm{LH},E}{E}_2},$(5) )–(20(20) $InhA=h_0+h_1\mathrm{DomF}+h_2{\mathrm{Lut}}_2+h_3{\mathrm{Lut}}_3+h_4{\mathrm{Lut}}_4.$(20) ) fitting the McLachlan data [23].

Param.	Value	Units
a	11	Dimensionless
V	2.5	L
$a_{\mathrm{LH}}$	14	1/day
$a_{\mathrm{FSH}}$	8.21	1/day
$v_{0\mathrm{LH}}$	1577.78	µ g/day
$v_{1\mathrm{LH}}$	25,288.06	µ g/day
$K{m}_{\mathrm{LH}}$	228.65	ng/L
$K{i}_{\mathrm{LH},P}$	3.01	µ g/L
$k_{\mathrm{LH}}$	12.26	1/day
$c_{\mathrm{LH},P}$	0.085	L/ µ g
$c_{\mathrm{LH},E}$	9.22E−08	L/ng
$d_E$	0.47	days
$d_P$	0.80	days
$v_{\mathrm{FSH}}$	5190.68	µ g/day
$K{i}_{\mathrm{FSH},\mathrm{InhA}}$	654.11	U/L
$k_{\mathrm{FSH}}$	12.83	1/day
$c_{\mathrm{FSH},P}$	2.46	L/ $\mu g$
$c_{\mathrm{FSH},E}$	0.0019	(L/ng)${}^2$
$d_{\mathrm{InhA}}$	1.66	days
α	0.31	Dimensionless
β	0.022	Dimensionless
γ	0.011	Dimensionless
b	0.60	L/day
$c_1$	1.59	1/day
$c_2$	0.51	(L/µ g)${}^{\alpha }$(1/day)
$c_3$	0.057	(L/µ g)${}^{\beta }$(1/day)
$c_4$	0.025	(L/µ g)(1/day)
$c_5$	0.30	(L/µ g)${}^{\gamma }$(1/day)
$d_1$	0.70	1/day
$d_2$	0.59	1/day
$k_1$	0.66	1/day
$k_2$	0.63	1/day
$k_3$	0.60	1/day
$k_4$	0.80	1/day
$K{m}_F$	115.18	µ g/L
p	19.80	Dimensionless
$c_{\mathrm{atr}}$	615.89	1/day
$K{i}_{\mathrm{atr}}$	4.69	µ g/L
q	2.36	Dimensionless
$e_0$	0.00	ng/L
$e_1$	0.023	1/kL
$e_2$	0.025	1/kL
$e_3$	0.14	1/kL
$p_0$	0.00	µ g/L
$p_1$	0.0063	1/L
$p_2$	0.0025	1/L
$h_0$	208.76	U/L
$h_1$	0.0040	(U/L)/µ g
$h_2$	0.00	(U/L)/µ g
$h_3$	0.66	(U/L)/µ g
$h_4$	0.00098	(U/L)/µ g

Note: All of the following optimized parameters were truncated to two significant decimal digits.

Table A.2. Initial conditions for the two-wave merged model Equations (5(5) $\frac{\mathrm{d}}{\mathrm{d}t}{\mathrm{RP}}_{\mathrm{LH}}=\frac{v_{0\mathrm{LH}}+\frac{v_{1\mathrm{LH}}{E}_2(t-d_E{)}^a}{K{m}_{\mathrm{LH}}^a+E_2(t-d_E{)}^a}}{1+\frac{P_4(t-d_P)}{K{i}_{\mathrm{LH},P}}}-\frac{k_{\mathrm{LH}}(1+c_{\mathrm{LH},P}{P}_4){\mathrm{RP}}_{\mathrm{LH}}}{1+c_{\mathrm{LH},E}{E}_2},$(5) )–(20(20) $InhA=h_0+h_1\mathrm{DomF}+h_2{\mathrm{Lut}}_2+h_3{\mathrm{Lut}}_3+h_4{\mathrm{Lut}}_4.$(20) ) used for simulations plotted in Figures 5– 7.

State	Value	Units
${\mathrm{RP}}_{\mathrm{LH}}\left(0\right)$	83.43	µ g
$\mathrm{LH}\left(0\right)$	25.34	µ g/L
${\mathrm{RP}}_{\mathrm{FSH}}\left(0\right)$	442.29	µ g
$\mathrm{FSH}\left(0\right)$	142.50	µ g/L
$\mathrm{RcF}\left(0\right)$	790.68	µ g
$\mathrm{GrF}\left(0\right)$	1160.43	µ g
$\mathrm{DomF}\left(0\right)$	113.01	µ g
${\mathrm{OvF}}_1\left(0\right)$	24.77	µ g
${\mathrm{OvF}}_2\left(0\right)$	15.19	µ g
${\mathrm{Lut}}_1\left(0\right)$	14.94	µ g
${\mathrm{Lut}}_2\left(0\right)$	37.76	µ g
${\mathrm{Lut}}_3\left(0\right)$	104.52	µ g
${\mathrm{Lut}}_4\left(0\right)$	140.12	µ g

Table A.3. Parameters for the three-wave ovarian component model Equations (9(9) $\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{RcF}=b\cdot \mathrm{FSH}+\frac{c_1{\mathrm{FSH}}^p}{K{m}_F^p+{\mathrm{FSH}}^p}\mathrm{RcF}-c_2{\mathrm{LH}}^{\alpha }\mathrm{RcF},$(9) )–(20(20) $InhA=h_0+h_1\mathrm{DomF}+h_2{\mathrm{Lut}}_2+h_3{\mathrm{Lut}}_3+h_4{\mathrm{Lut}}_4.$(20) ) fitting the McLachlan data [23].

Parameter	Value	Units
α	0.31	Dimensionless
β	0.20	Dimensionless
γ	0.035	Dimensionless
b	1.86	L/day
$c_1$	1.31	1/day
$c_2$	0.42	(L/µ g)${}^{\alpha }$(1/day)
$c_3$	0.14	(L/µ g)${}^{\beta }$(1/day)
$c_4$	0.035	(L/µ g)(1/day)
$c_5$	0.39	(L/µ g)${}^{\gamma }$(1/day)
$d_1$	0.59	1/day
$d_2$	0.56	1/day
$k_1$	0.56	1/day
$k_2$	0.63	1/day
$k_3$	0.60	1/day
$k_4$	0.58	1/day
$K{m}_F$	122.10	µ g/L
p	19.09	Dimensionless
$c_{\mathrm{atr}}$	285.45	1/day
$K{i}_{\mathrm{atr}}$	8.58	µ g/L
q	1.83	Dimensionless
$e_0$	16.00	ng/L
$e_1$	0.01	1/kL
$e_2$	0.045	1/kL
$e_3$	0.13	1/kL
$p_0$	0.00	µ g/L
$p_1$	0.011	1/L
$p_2$	0.001	1/L
$h_0$	188.05	U/L
$h_1$	0.087	(U/L)/(µ g)
$h_2$	0.01	(U/L)/(µ g)
$h_3$	0.95	(U/L)/(µ g)
$h_4$	0.18	(U/L)/(µ g)

Note: All of the following optimized parameters were truncated to two significant decimal digits.

Table A.4. Superfecundation parameters for the merged, two-wave model fitting the McLachlan data [23].

Param.	Value	Units
a	11	Dimensionless
V	2.5	L
$a_{\mathrm{LH}}$	14	1/day
$a_{\mathrm{FSH}}$	8.21	1/day
$v_{0\mathrm{LH}}$	212.68	µ g/day
$v_{1\mathrm{LH}}$	35,656.90	µ g/day
$K{m}_{\mathrm{LH}}$	170.00	ng/L
$K{i}_{\mathrm{LH},P}$	2.81	µ g/L
$k_{\mathrm{LH}}$	11.41	1/day
$c_{\mathrm{LH},P}$	0.084	L/µ g
$c_{\mathrm{LH},E}$	0.0020	L/ng
$d_E$	0.50	days
$d_P$	1.00	days
$v_{\mathrm{FSH}}$	4500.00	µ g/day
$K{i}_{\mathrm{FSH},\mathrm{InhA}}$	646.30	U/L
$k_{\mathrm{FSH}}$	11.55	1/day
$c_{\mathrm{FSH},P}$	3.16	L/µ g
$c_{\mathrm{FSH},E}$	0.0023	(L/ng)${}^2$
$d_{\mathrm{InhA}}$	1.74	days
α	0.30	Dimensionless
β	0.020	Dimensionless
γ	0.010	Dimensionless
b	0.60	L/day
$c_1$	1.56	1/day
$c_2$	0.55	(L/µ g)${}^{\alpha }$(1/day)
$c_3$	0.051	(L/µ g)${}^{\beta }$(1/day)
$c_4$	0.027	(L/µ g)(1/day)
$c_5$	0.60	(L/µ g)${}^{\gamma }$(1/day)
$d_1$	0.70	1/day
$d_2$	0.70	1/day
$k_1$	0.60	1/day
$k_2$	0.65	1/day
$k_3$	0.70	1/day
$k_4$	0.80	1/day
$K{m}_F$	115.00	$\mu g$/L
p	15.00	Dimensionless
$e_0$	0.00	ng/L
$e_1$	0.025	1/kL
$e_2$	0.047	1/kL
$e_3$	0.060	1/kL
$p_0$	0.00	µ g/L
$p_1$	0.0047	1/L
$p_2$	0.00016	1/L
$h_0$	140.00	U/L
$h_1$	0.020	(U/L)/$\mu \mathrm{g}$
$h_2$	0.00	(U/L)/$\mu \mathrm{g}$
$h_3$	0.63	(U/L)/$\mu \mathrm{g}$
$h_4$	0.020	(U/L)/$\mu \mathrm{g}$

Note: All of the following optimized parameters were truncated to two significant decimal digits.