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.

S(t) ≥ 0 denotes the influencing substance. The threshold T > 0 represents the amount of S at which the Hill function has the value 1/2, i.e., where 50% of the maximum effect is reached. n ≥ 1 is called the Hill coefficient and regulates the rate of switching. i denotes the influencing species and j the influenced species.
Equations S1-4 deal with LH dynamics in two compartments. LH synthesis is described by a basal synthesis rate constant b LH Syn and is modulated by two Hill equations (stimulation by E2 and inhibition by P4).
Syn LH (t) = (b LH Syn + k LH E2 · H + (E2(t), T LH E2 ; n LH E2 )) · H − (P 4(t), T LH P 4 ; n LH P 4 ) (S1) LH is released from the pituitary with a basal release rate constant b LH Rel , which depends on the amount of LH in the pituitary and which is stimulated by the GnRH-receptor complex.
The amount of LH in the pituitary can be described by the difference between the synthesized and the released amount of LH.
The released amount of LH is diluted in the blood volume V blood . Through the blood stream, LH reaches its receptor and forms complexes. LH is removed from the systems with the clearance rate constant k LH cl .
FSH synthesis is inhibited by progesterone and the frequency of pulsatile GnRH release.
The fundamentals of the equations S6-8 are analogue to S2-4.
FSH reaches the ovaries through the blood stream and binds to its receptors on the follicular surfaces.
The dynamics of the FSH receptors is described in three steps (S10-12). FSH binds to its receptor R F SH , resulting in the formation of an active complex F SH-R. Desensitised receptors R F SH,des are caused by the dissociation of the complex. Each step has a characteristic reaction rate constant.
f req and mass characterize the GnRH pulse generator needed to describe the GnRH dynamics in the system. The basal frequency f 0 is modulated by the stimulatory effect of E2 and the inhibitory effect of P4.
The basal released mass a 0 is regulated by E2 stimulation as well as inhibition depending on the E2 concentration.
The GnRH receptor is present in four states in the model (S16-19). Free active receptors R G,a are available for binding GnRH with a forward rate k G on and a reverse rate k G of f . The receptors and receptor complexes can switch between active and inactive states.
The concept to simulate follicular maturation is derived from Lange et al. (2018). The Hill term in equation (S20) represents the stimulation of follicular growth and competition as soon as the individual FSH sensitivity threshold T F SH-R (i) is exceeded.
The growth rate γ is inhibited by P4 and stimulated by the FSH receptor complex, whereas κ is inhibited by the FSH receptor complex, The P4 concentration in the model is described by a Gaussian with three parameters c 1 , v 1,2 that characterize its height, mean and variance and whose values were determined by fitting the curve to P4 data, see Fig. S1. The mean is located with respect to the last time point of ovulation, T ovu .
We assume that all follicles with a radius larger than T F S contribute to E2 production, and that the amount of E2 produced depends linearly on their surface F S.
The second summand represents the amount of E2 produced by the corpus luteum.
Simulating ovarian stimulation with r-FSH, which is structurally identical to FSH, is based on the idea that the administration of a drug which is identical to a species in the system can be modelled as an additional time dependent source term. The following algebraic equation is used to calculate the drug concentration c(t) in the system.
The three parameters D, β and c L were calculated uniquely for each administered drug from the corresponding values of three pharmacokinetic parameters given in (Kompendium, 2020). For this work, equation (S25) is used to calculate the concentration of r-FSH, which is then added to the concentration of FSH to form the total FSH concentration in the system.

LIST OF PARAMETERS
The following table includes all model parameters, their values and short descriptions. Hill terms and coefficients are not described in further detail, since they are represented in a consistent notation which is introduced at the beginning of the first section (1 EQUATIONS) of the supplementary material. Parameters marked with * are either newly introduced or modified compared to Röblitz et al. (2013 Figure S2. Representation of the cycle length of 42 simulated menstrual cycles.