Mathematical Modeling and Simulation Provides Evidence for New Strategies of Ovarian Stimulation.

New approaches to ovarian stimulation protocols, such as luteal start, random start or double stimulation, allow for flexibility in ovarian stimulation at different phases of the menstrual cycle. It has been proposed that the success of these methods is based on the continuous growth of multiple cohorts ("waves") of follicles throughout the menstrual cycle which leads to the availability of ovarian follicles for ovarian controlled stimulation at several time points. Though several preliminary studies have been published, their scientific evidence has not been considered as being strong enough to integrate these results into routine clinical practice. This work aims at adding further scientific evidence about the efficiency of variable-start protocols and underpinning the theory of follicular waves by using mathematical modeling and numerical simulations. For this purpose, we have modified and coupled two previously published models, one describing the time course of hormones and one describing competitive follicular growth in a normal menstrual cycle. The coupled model is used to test ovarian stimulation protocols in silico. Simulation results show the occurrence of follicles in a wave-like manner during a normal menstrual cycle and qualitatively predict the outcome of ovarian stimulation initiated at different time points of the menstrual cycle.

Introduction

Infertility is a worldwide problem. According to the World Health Organization, about 48.5 million couples worldwide were affected by unwanted childlessness in 2010, and the number continues to grow (1). Men and women are just as likely to contribute to the couple’s infertility (2). Infertility as a disease of the female reproductive system affects approximately 10% of women of reproductive age worldwide (3). Unbalanced hormone levels are one cause, in a wide range of conditions, leading to infertility. For many couples, unwanted childlessness is a burden. Assisted reproductive technologies (ART) provide strategies to deal with infertility. Both unwanted childlessness and ART increase the risk for negative psycho-social functioning, such as depression and anxiety disorders (4–6), whereby the treatment burden has fallen mainly on women (2). Therefore, new ART approaches deserve to be highlighted. We want to add further scientific evidence for the efficiency of those new approaches by using mathematical modeling and numerical simulations.

Female reproduction is essentially enabled by a feedback mechanism between ovarian hormones, mainly progesterone (P4) and estradiol (E2), and the pituitary hormones luteinizing hormone (LH) and follicular stimulating hormone (FSH), see . The hormone interaction network is important for regulating folliculogenesis. While the initial recruitment of follicles does not depend on gonadotropins (7, 8), the growth of cohorts of larger follicles relies on a stimulatory effect of FSH. FSH signaling is mediated by the expression of FSH receptors on granulosa cells (9, 10). The gonadodropins LH and FSH are responsible for follicular estradiol production. LH stimulates androstenedione production, which is the substrate for the FSH stimulated aromatase reaction producing estradiol (8, 11, 12). Around mid-cycle, usually one dominant follicle ovulates and releases an oocyte. The remaining parts of the dominant follicle transform into the corpus luteum, which has a key role in preparing the body for a possible pregnancy. If the oocyte is not fertilized, the corpus luteum decays and a new cycle starts (13–15). Interruptions in the feedback system are one reason for infertility.

Figure 1

Flowchart illustrating the interactions included in the given model. This is a simplified feedback interactions network for the hormonal control of the female menstrual cycle. Green arrows indicate positive feedback effects, while red arrows express negative feedback interactions. Gray arrows show other types of interactions. The pulsatile release of GnRH stimulates the release of the pituitary hormones LH and FSH. These hormones effect follicular maturation. Growing follicles produce E2 which has a positive feedback effect on the LH concentration. A high LH concentration triggers the ovulation of one selected follicle (light gray arrow) followed by the formation of the corpus luteum (dark gray dashed arrow). The simultaneous release of E2 and P4 by the corpus luteum (dark gray arrows) inhibits the release of GnRH. Additionally, P4 has an inhibitory effect on LH and FSH. While P4 only has an inhibitory effect on GnRH, E2 has either a stimulatory or an inhibitory effect on GnRH, depending on the E2 level.

Modern assisted reproductive technologies like in vitro fertilization (IVF) or intracytoplasmic sperm injection (ICSI) have increased the chance for pregnancy. Ovarian stimulation, which aims at obtaining multiple fertilizable oocytes, is a critical step in ART (16). Since the 1980s, the long gonadotropin-releasing hormone (GnRH) agonist protocol has been commonly used to prepare for oocyte retrieval and in-vitro fertilization (17, 18). This protocol starts around mid-luteal phase with GnRH agonist administration for about 14 days. Right after the beginning of GnRH agonist administration, a short period of gonadotropin (FSH and LH) hypersecretion is observable. The treatment leads to GnRH-receptor down-regulation in the pituitary (19, 20). In the next step, the growth of multiple follicles is stimulated by FSH administration alone, e.g. with recombinant FSH (rFSH), or by a combination of FSH and LH, e.g. with human menopausal gonadotropin (hMG). Continuation of GnRH agonist administration during the stimulation phase prevents an LH surge and hence ovulation. In the final step, ovulation is induced by injecting human chorionic gonadotropin (hCG) (18). Patient-specific and clinic-dependent modifications of these general procedures are common. The two most common alternatives are the short GnRH agonist protocol and the antagonist protocol. Both protocols work without downregulation, though some clinics perform a pre-treatment phase for 10 to 25 days with a P4 antagonist that inhibits ovulation.

The stimulation phase in the short GnRH agonist protocol is the same as in the long protocol. It includes the stimulation with hMG or rFSH and the concurrent administration of a GnRH agonist. The antagonist protocol also includes the stimulation with hMG or rFSH but, in contrast to the agonist protocols, a GnRH antagonist is administered from day 5 of the stimulation period. The final step in all protocols is the induction of ovulation by hCG.

In general, infertility treatment is a long-term and expensive therapy with high dropout rates (21), mainly because it imposes physical, mental, and emotional burdens (22). Often, life has to be subordinated to medical procedures. Therefore, treatment alternatives are of interest. Both the short and the antagonist protocol are less time-consuming than the long protocol. However, the stimulation phase in these protocols conventionally starts in the early follicular phase. This constraint could cause too long waiting times, e.g. for women requiring emergency fertility preservation. Hence, the advancement of a new class of ovarian stimulation approaches called random - and luteal phase-start ovarian stimulation protocol (23) has progressed. In recent years, several studies investigating ovarian stimulation protocols with various starting points have been published (24–26). Originally, these protocols were invented for fertility preservation in cancer patients, where time is a limiting factor (27). However, they might be beneficial for patients outside an oncological setting (23), though there is an ongoing debate whether the oocyte quality differs between protocols. Other approaches like the double ovarian stimulation, where two waves within one cycle are stimulated, might help to increase the number of accumulated oocytes within one treatment cycle (28). That strategy could be of particular interest for the therapy of poor ovarian response patients (29, 30).

One possible explanation for the success of stimulation initiated in different phases of the cycle is the “wave” theory. The use of high-resolution transvaginal ultrasonography has underpinned the hypothesis that, similar to ruminants, follicular growth and development in human is characterized by waves (31, 32), whereby each wave involves the recruitment of a cohort of follicles and the possible selection of a dominant follicle. Given that multiple waves of follicles appear each cycle, there are multiple time points during one cycle that are suitable to start ovarian stimulation.

The mathematical model underlying this study simulates the time-evolution of key hormones and growth behavior of multiple follicles. In particular, we test the hypothesis that random recruitment of follicles leads to the emergence of follicular waves. Based on the occurrence of follicular waves that we observe in our simulation results, we study variable-start ovarian stimulation protocols in silico. We demonstrate simulation results for two protocols, namely (i) stimulation initiated in the late follicular phase and (ii) stimulation initiated in the luteal phase. We analyze statistics of treatment duration and numbers of follicles in our simulation results and compare them with the literature.

Materials And Methods

Mathematical Modeling of the Female Menstrual Cycle

Mathematical modeling is a useful tool to better understand the human menstrual cycle by validating or testing hypothesis in silico, and predicting possible dynamics. A first mathematical model for the human menstrual cycle was introduced in a series of articles by Schlosser, Selgrade, and Harris-Clark (33). Their model allows to simulate the time course of hormones and follicular maturation stages over several cycles and is able to display multiple follicular waves (34). This model was extended by pharmacokinetic sub-models to simulate the administration of drugs, including ovarian contraceptive pills (35, 36) and GnRH analogs (37). These pharmacokinetic-pharmacodynamic (PKPD) models allow to study the influence of dose and time point of administration of various drugs on the cycle dynamics.

All those models are based on ordinary or delay differential equations since they allow to simulate the time evolution of hormone concentrations and follicles. Hill functions have been used to characterize stimulatory and inhibitory effects, as it is common practice for modeling regulatory networks. The model by Röblitz et al. (37) consists of 33 ordinary differential equations that describe the feedback mechanisms between the hormones that are of particular importance for the female menstrual cycle (GnRH, FSH, LH, E2, P4, inhibin A, inhibin B) and the development of follicles and corpus luteum throughout consecutive cycles. Compared to previous models, it does not use delay differential equations and consists of fewer equations and parameters. However, all those models have in common that follicular growth is described in terms of activity levels of different follicular maturation stages, but not in terms of follicle numbers and sizes. Thus, the simulation results cannot be compared with ultrasound data.

A mathematical model that quantifies the time evolution of the sizes of multiple follicles comparable to observations by ultrasound measurements in mono-ovulatory species was presented by (38). This model contains a separate differential equation for each follicle, whereby the structure of this equation is the same for all follicles, but the initial follicle sizes are different. The equations are coupled via a term that accounts for competitive interactions between follicles. Together with the model by (37) a previous version of the model by (38) formed the basis for the development of computational tools to enable in silico clinical trials in reproductive endocrinology (39, 40). In particular, by introducing variability into model parameters (41–43), the authors could analyze inter-individual variability in the cycle and automatically synthesize, by means of artificial intelligence guided by patient digital twins, optimal personalized treatments for the patients at hand (44). However, the tools could only be applied to the downregulation phase before follicular stimulation, because the feedback mechanisms from the ovaries to the pituitary were not implemented in the modified model. This drawback motivated the development of the fully coupled model as presented in this work. To our knowledge, this is the first mathematical model that allows for the simulation of stimulation protocols that start at different time points in the cycle.

Model Construction and Assumptions

The mathematical model underlying this work is the result of modifying and coupling the two previously published models by Röblitz et al. (37) and Lange et al. (38). In a first step, the model by Röblitz et al. (37) was reduced by removing the equations for the development of follicles and the corpus luteum and the hormones produced by them (inhibin A, inhibin B, E2, P4). In addition, we removed the equations for LH receptor binding mechanisms, since they were not needed for our purpose. The remaining equations were kept exactly as in (37), except for the FSH synthesis rate. In the new model, this rate is inhibited by P4 instead of inhibin A and B [ ], since P4 reaches its peak in the mid-luteal phase exactly as inhibin A. The influence of inhibin B could be neglected without any consequences for the qualitative behavior of the model. In addition, we have introduced a new equation for the amount of FSH that reaches the follicles [ ] to account for delays caused by transportation and for changes in concentration caused by different volumes. In contrast to (37), the equations for FSH receptor binding now take into account FSH in the ovaries instead of the FSH blood concentration [ ].

Instead of re-introducing a corpus luteum into the model equations, we decided to use algebraic equations to directly model the amounts of E2 and P4 produced in the luteal phase of the cycle [ ]. The model describes E2 and P4 levels in the luteal phase by Gaussian-shaped curves with fixed parameters based on fits to experimental data (for P4 see ). This simplification is based on the observation that the variability in the length of the luteal phase is significantly lower than the variability in the length of the follicular phase (45).

We modified the follicle equation introduced by (38) in a way that the hormone dynamics in the system have a direct effect on the follicular growth behavior [ ]. The maturation of each follicle is modeled by a separate ODE. All ODEs have the same structure and include both shared and follicle specific parameters. Each follicle carries two random properties that are follicle specific, hence there are two follicle specific parameters: the time point at which a follicle is recruited, and its FSH sensitivity. The following assumptions are made about these two parameters:

The time point at which a follicle is recruited and starts growing is follicle-specific and follows a Poisson process. The overall number of follicles that are recruited within a specific time interval is a Poisson random variable. The parameter of this distribution, in the following named Poisson parameter, corresponds to the probability that a given number of follicles is recruited in a fixed time interval. In the model, the Poisson parameter is modulated by the FSH concentration: if the FSH concentration is above a certain threshold, more follicles are recruited.

The second property is a follicle specific FSH value, referred to as FSH sensitivity threshold value, which has to be exceeded in order to stimulate the follicle’s growth. This refers to the biological finding that follicle growth does not occur below a certain level of FSH (46), and that any two follicles might respond differently to FSH, even if the two have the same size, because they differ in the FSH receptor density. The distribution of the FSH sensitivity threshold values in the population of follicles is assumed to follow a normal distribution. Follicles that are more sensitive to FSH, i.e. which require less FSH to start growing, have a competitive advantage for being selected as the dominant follicle. Whether a follicle becomes dominant, however, depends on both its FSH sensitivity and its recruitment time point.

The competition between follicles, which is represented by a common parameter [ ], is inhibited by FSH concentrations above a certain threshold, taking into account the “FSH window concept” (47–49). This concept is based on the observation that the period of time during which FSH is above a certain threshold effects the number of follicles reaching the dominant follicle’s size (50, 51). Moreover, we assume that the follicular growth rate is inhibited by P4 and stimulated by the FSH receptor complex level in a threshold dependent way [ ] (52).

Growing follicles are the main source of E2 in the female body and the dominant follicle produces the most E2 (12, 53, 54). Estradiol is produced by granulosa cells, which proliferate and form a multilayered structure. This is included in the model by an additional term in E2 production which is dependent on the follicular size [ ].

To sum up, the coupling between the hormone dynamics model and the follicular growth model is realized as follows (compare ). The levels of FSH in the blood and of the FSH receptor complex enter into the equations for the follicles in a threshold dependent way. In addition, the LH level plays a role in determining the time point of ovulation. Ovulation of a follicle that exceeds the size threshold occurs 12 h after the LH level is higher than a certain threshold. The levels of E2 and P4 in the luteal phase depend on the time point of the last ovulation. E2 and P4 levels enter into the equations for LH and FSH synthesis and for the frequency and mass of GnRH, in the same way as in (37). The coupled model contains in total 72 parameters, i.e. less than the two original models taken together (114 parameters in (37) and 5 parameters in (38). We adopted 44 parameters from (37) and only changed the values of three of them. A detailed parameter list can be found in the Supplement. The model has been implemented in MATLAB and numerical simulations were performed using the ODE solver ode15s. The code is available at https://github.com/SoFiwork/GynCycle.

Ovarian Stimulation Protocols

Stimulation protocols are introduced to the model by a pharmacokinetic approach. The dosing concentrations of the administered drug, as used in the ovarian stimulation protocols, are calculated during the simulation based on three drug specific pharmacokinetic parameters using the information given by (55) [ ]. In order to study treatment outcomes, two different stimulation protocols were implemented. The two studies were selected based on the accessibility of results, the size of study cohorts and the physiological stage of patients. Each study includes data from more than 100 women. Patients were at the age of 18 to 40 years with a body mass index of 18 to 30 kg/m3. All women showed spontaneous ovulation.

Stimulation Initiated in the Late Follicular Phase

Our simulated treatment protocol for ovarian stimulation during the late follicular phase follows the description in Zhu and Fu (24). As a simplification, we did not vary the administered hMG dose during the first days of stimulation. The stimulation starts with a daily administration of 150 IU hMG when at least one follicle measures 14 mm in diameter. After 6 days, the daily dose is increased to 225 IU per day. We chose day 6 to change the hMG concentration because re-examination and dose adjustment in the clinical trial took place after 5 - 7 days. The stimulation stops whenever at least 3 follicles reach a diameter of at least 18 mm. The ovulation of a dominant follicle during the stimulation phase is characteristic for this protocol.

Stimulation Initiated in the Luteal Phase

The protocol described in (26) served as a reference to simulate the stimulation of multiple follicular growth during the luteal phase. In this clinical trial, the drug administration in the simulation starts between day 1 and 3 after ovulation under the condition that there exist follicles smaller than 8 mm. Follicular growth is stimulated by the daily administration of 225 IU hMG. The stimulation terminates if at least three follicles have reached a diameter of 18 mm.

Results

Unstimulated Cycle

As indicated in , the model generates quasi-periodic solutions for all four hormones. Due to the individual growth behavior of follicles implemented in the model, variations in cycle length and number of follicles per cycle occur. Simulations for a normal cycle were performed for more than 1000 time steps in order to get an idea of the variability in the model outcome. In total, 42 simulated menstrual cycles (here, one menstrual cycle is defined from one ovulation to the next one) were used for a statistical analysis. In the simulations, the average cycle length was 30.56 days, with a standard deviation of 7.00 days ( ). On average, 16.19 follicles greater than 4 mm were detected during one cycle, with a standard deviation of 3.08 follicles. The results were tested for normality using the Shapiro-Wilk test with a 95 confidence interval. A correlation between the cycle length and the follicular count was not observed.

Figure 2

Simulation results of the female menstrual cycle model are displayed. The left column illustrates the simulation outcome for two menstrual cycles and the right column zooms into details. Here, one cycle is defined from one ovulation to the next one. Sub-figures (A–D) represent the simulated hormone concentration profiles for LH, FSH, E2 and P4. (E) portrays growth trajectories of follicles >4 mm. The ovulation of a dominant follicle is indicated by terminating trajectories, as seen for example around day 80 of the simulation. (F) illustrates competition between follicles indicated by crossing growth trajectories. (G) Points out that the ovulation of a dominant follicle occurs 12 h after the LH peak concentration as a result of the way the ovulation process is implemented in the model.

The simulated hormone curves are supposed to be comparable to serum hormone concentration profiles in terms of shape and peak values. display consecutive menstrual cycles in the time period between day 50 and day 130 from one simulation run. The time evolution of all four hormone profiles is illustrated, and the described interplay between hormones and follicles is apparent.

The wave-like growth behavior of the follicles ( ) is generated by the model itself and is not enforced by the implementation. shows an example of the ovulation of a dominant follicle that occurs 12 h after LH reached its peak concentration. This 12-h gap is accomplished by the way the ovulation event is defined in our model (see Discussion). Once ovulation is detected during the run time of the simulation, the ovulated follicle is taken out from the cohort of follicles (indicated by the terminating trajectory in ). This follicle no longer contributes to steroid production. Keeping it in the simulation would needlessly increase computational time. The growth behavior of follicles causes variation in the length of the follicular phase. In contrast to that, the luteal phase has a constant length of 14 days due to its implementation.

The follicular growth equation, as introduced by (38) and modified for the given model, includes a term addressing the competition for dominance between follicles. In the simulation results, its effect is visible by crossing growth trajectories ( ). This crossing only is possible because each follicle has its specific parameters. As it can be seen in , competition is stronger during the early follicular phase before a dominant follicle emerges.

Ovarian Stimulation

The simulations of ovarian stimulation initiated in the luteal phase or the late follicular phase are characterized by the growth of multiple follicles. Additionally, the ovulation of a dominant follicle during a stimulation protocol occurs only during stimulation in the late follicular phase. In the model, the competition term is inhibited by high FSH concentrations, enabling the growth of multiple follicles under stimulatory treatment.

exemplarily displays hormone concentration profiles and follicle development for one simulation of each treatment approach. Additionally, error bars at four characteristic time points (one day before treatment, one day after first drug administration, six days after first drug administration, last day of drug administration) indicate the variability in the hormone levels between 20 simulations using the same treatment conditions. The characteristic time points where chosen in a way that the results are easily comparable to the clinical data. In both cases, the FSH concentration rises with each day of the treatment. Due to the growth of multiple large follicles, which are the main source of E2, the E2 level increases significantly during ovarian stimulation. The levels are almost ten times higher compared to the normal cycle ( ).

Figure 3

Simulation results for two different ovarian stimulation protocols. The growth of multiple large follicles, caused by the stimulation treatment, is characteristic for both strategies. The left column represents simulation results from a luteal phase stimulation protocol, while the right column shows the effect of a stimulation during the late follicular phase. Sub-figures (A–D, F–I) exemplary represent hormone profiles originating from one simulation in red. Purple dots and error bars represent mean values and variances, respectively, from 20 simulations at four characteristic time points: 1 day before the stimulation treatment starts, 1 day after starting the treatment, 6 days after starting the treatment, and the last day of treatment. Sub-figures (E, J) illustrate the growth trajectories of the follicles.

Simulations of an ovarian stimulation during the luteal phase are dominated by high P4 levels during the stimulation with hMG. The high P4 concentration prevents the ovulation of follicles (through the negative feedback mechanisms of P4 on LH). The concentrations of LH, FSH, P4 and E2 in are comparable to observations by (26).

illustrates the follicular growth behavior under stimulation in the late follicular phase, initiated after the occurrence of a dominant follicle. The ovulation of the dominant follicle is followed by an increase in P4 concentration comparable to non-treated conditions. The E2 level decreases after the ovulation of the dominant follicle but starts to increase again. This increase is caused by multiple large follicles as a result of the stimulation.

represents the individual outcomes (treatment duration and follicular count) of 20 simulations per treatment protocol. The mean and standard deviation of these results are given in . The simulation results for ovarian stimulation initiated in the luteal phase match the observations from Kuang et al. (26). The simulated treatment duration for the late follicular phase stimulation approach is noticeably lower than the clinical observations, which goes along with comparably low counts of follicles >14 mm. convincingly shows that simulations differ among each other even if non-follicular parameters are the same in all simulations. Hence, the individual growth behaviors of the follicles have a major effect on treatment simulations and outcomes.

Figure 4

Simulation outcomes of 20 independent cycles for each treatment: ovarian stimulation induced either during the luteal phase (top, A, B) or the late follicular phase (bottom, C, D). In the upper row (A, B), follicular counts and treatment duration for the luteal phase stimulation approach are displayed (red: follicles 10–14 mm; purple: follicles >14 mm). On average, 11.1 ± 3.5 follicles with a diameter of 10 to 14 mm and 8.9 ± 3.7 follicles with a diameter >14 mm are observed. The average treatment duration is 9.4 ± 0.7 days. The lower row (C, D) shows follicular counts and treatment durations for simulated stimulations in the late follicular phase. A treatment cycle takes about 6.0 ± 0.7 days. The average count of follicles with diameters 10 - 14 mm is 6.3 ± 2.2 and the one for follicles >14 mm is 8.0 ± 2.2. (Numbers refer to mean ± standard deviation.)

Table 1

Comparison between simulation results and clinical observations.

	Luteal phase ovarian stimulation		Late follicular phase ovarian stimulation
	Kuang et al. (26)	Simulation	Zhu and Fu (24)	Simulation
Num. of follicles with	13.9 ± 7.8	11.1 ± 3.5		6.3 ± 2.2
diameter 10 - 14 mm
Num. of follicles with	11.1 ± 5.5	8.9 ± 3.7	11.7 ± 6.2	8.0 ± 2.2
diameter > 14 mm
Duration of treatment	10.2 ± 1.6	9.4 ± 0.7	10.93 ± 1.66	6.0 ± 0.7
with hMG

Ovarian stimulation is induced either during the luteal phase or the late follicular phase. Each of the two studies includes data from more than 100 woman. Patients were at the age of 18 – 40 years with a body mass index of 18 – 30 kg/m3. All woman showed spontaneous ovulation.

Discussion

The mathematical model developed in this work addresses the interplay between pituitary hormones, ovarian hormones and follicular growth. Simulation results for the unstimulated cycle agree qualitatively and quantitatively with observations reported in literature. In particular:

The time evolution of the four hormone profiles for LH, FSH, P4 and E2 is consistent with the scientific literature (56).

An average cycle length of around 29 days, ranging from cycles with a duration of 22–25 up to 36 days, is reported in experimental studies (56–58). The simulation results are in line with these observations.

In the literature, it is described that the variability in the length of the follicular phase is significantly higher than for the luteal phase (58, 59). The given simulation results fulfill the same property.

The observed intra-cycle variability of 7 days is comparable to experimental results by (58).

(32) observed the emergence of two to three waves carrying 4 to 14 follicles greater than 4 mm. The given simulation results of 16.19 ± 3.08 follicles in two waves per cycle match their experimental investigations.

The discontinuity in the profile of the E2 curve ( at day 85 of the simulation) is related to the growth behavior of the follicles and is caused by atresia of larger sub-dominant follicles.

By comparing the results in Table 4, it is visible that variations in the experimental data are higher than in the simulation results. That indicates the fact that the inter-individual variability in human is higher than the variability between simulations sharing one set of non-follicular parameters. The stochastic growth behavior of follicles is the only source of variability between simulations. According to (24), the LH concentration under stimulatory treatment in the late follicular phase is not supposed to increase after the ovulation of the dominant follicle due to the inhibitory effect of P4. However, this effect is not visible in the simulation results ( ). This might be due to the comparably lower P4 concentrations in the simulation results. Here, the P4 concentration at day 6 is about 0.99 ± 0.6 ng/mL, whereas the figures published by (24) indicate P4 concentrations up to more than five times as high. In the present model, the P4 concentration is linked to the formation of the corpus luteum as the only source of P4. Minor P4 sources such as the adrenal cortex are neglected. However, the equations for the P4 concentration matches experimental measurements quite well ( in Supplementary Material). A relation between the high LH concentrations, the low P4 concentrations and the follicular growth behavior are conceivable as well. Since the simulated treatment duration is several days shorter than those in the clinical observations, it appears that follicles are growing too fast during the simulation of ovarian stimulation. If this is the reason for the mismatch between the simulation results and the observations by (24), two explanations are credible: (i) the model parameters should have other values, or (ii) at least one mechanism is missing. However, at this point it was not possible to compare the simulated follicular growth under treatment to detailed experimental investigations since ultrasound measurement data were not available from literature.

Another reason for the mismatch could be that we could not simulate the clinical treatment procedures in full detail. In a clinical setting the dose is adjusted according to the treatment response, which is based on an evaluation of follicular growth during the stimulation procedure. Since the criteria for dose adjustment were not described in the available publications, we did not implement adjustments in our model.

We have not yet simulated double ovarian stimulation due to technical difficulties with the model implementation. However, we will do this in future work in order to address some of the problems that are still unsolved (60), for example the choice of the best day to start the second stimulation or the necessity of using a GnRH antagonist during the second stimulation.

Finally, we want to point out that clinical data are mainly reported as summary statistics, usually in terms of means and standard deviations, and for very few indicators, e.g. treatment duration or number and sizes of follicles on certain treatment days. However, with our model-based approach we could go beyond a simple comparison of moments. Since the model simulations generate distributions, we could compare them with data from literature if the publications about clinical trial outcomes reported the complete data distributions.

Conclusion

This study demonstrates how mathematical modeling and simulations can contribute to enhance our mechanistic understanding of ovarian stimulation protocols. In particular, our approach allows to study the extend of variability in both treated and untreated cycles. The model simulations confirm that follicular size is not a reliable parameter for determining treatment outcome since the receptor status of each individual follicle (modeled by the FSH sensitivity threshold) and the timing of growth matter. However, we cannot (yet) make use of that knowledge in a clinical setting as long as the receptor status cannot be inferred from measurements. Making predictions on the level of individuals, either in-vivo or in-silico, will therefore remain notoriously difficult. However, models that include random effects can be used to quantify uncertainties in the predictions. Even though these uncertainties might be large, being aware of what could happen as well as identifying outliers can assist in making decisions. Moreover, the model presented here could be used to compare the outcome of different treatment strategies in terms of specific success criteria (e.g. average number of follicles larger than a threshold size at the end of the stimulation), similar to the approach in (39). This requires to first validate the model with data from other stimulation protocols. For example, in order to compare the two protocols simulated here with the three currently most often used protocols (long, short, and antagonist), we would need data on each protocol from cohorts that are comparable in terms of size and physiological stage (e.g. race, age, BMI). We therefore invite clinicians to share their data and to join interdisciplinary research projects with the ultimate goal to develop model-based clinical decision support systems.

Data Availability Statement

The original contributions presented in the study are included in the article/ . Further inquiries can be directed to the corresponding author.

Author Contributions

SF and SR conceived the study. ET, TM, ME, FI, TK, BL, and SR obtained the funding. BL collected the data. SF, RE, SS, BL, and SR analyzed and interpreted the data. SF, RE, SS, and SR developed the mathematical model. SF and SS implemented the model and performed the simulations. All authors contributed to the article and approved the submitted version.

Funding

The work of SF and SR was supported by the Trond Mohn Foundation (BSF, https://www.mohnfoundation.no/), grant no. BFS2017TMT01. The funder had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

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.