H M T Boer

C Stötzel

S Röblitz

P Deuflhard

R F Veerkamp

H Woelders

A simple mathematical model of the bovine estrous cycle: Follicle development and endocrine interactions

Keywords: cow, reproduction, hormone patterns, differential equations, systems biology 2000 MSC: 92C42, 92C30, 90C31, 65L09

Bovine fertility is the subject of extensive research in animal sciences, especially because fertility of dairy cows has declined during the last decades.

The regulation of estrus is controlled by the complex interplay of various organs and hormones. Mathematical modeling of the bovine estrous cycle could help in understanding the dynamics of this complex biological system.

In this paper we present a mechanistic mathematical model of the bovine estrous cycle that includes the processes of follicle and corpus luteum development and the key hormones that interact to control these processes. The model generates successive estrous cycles of 21 days, with three waves of follicle growth per cycle. The model contains 12 differential equations and 54 parameters. Focus in this paper is on development of the model, but also some simulation results are presented, showing that a set of equations and

Introduction

Systems biology is a relatively new research area in the field of animal sciences. It aims at understanding how the various components of a biological system function together, rather than investigating only individual parts.

One approach is the translation of a conceptual biological model into a set of mathematical equations that represent the dynamic relations between system components. The purpose of building such mathematical models is to interpret and predict the dynamics of complex biological systems, and to identify new research questions.

One example of a dynamic biological system is the bovine estrous cycle, the hormonally controlled recurrent periods when the cow is preparing for reproduction by producing a fertilizable oocyte. Concurrent with selection for increased milk yield, a decrease in dairy cow fertility has been observed during the last decades (for reviews see [START_REF] Pryce | Fertility in the high-producing dairy cow[END_REF][START_REF] Veerkamp | Effects of genetic selection for milk yield on energy balance, levels of hormones, and metabolites in lactating cattle, and possible links to reduced fertility[END_REF]). This decline in fertility is shown by e.g. alterations in hormone patterns during the estrous cycle, reduced expression of estrous behavior and lower conception rates [START_REF] Wiltbank | Changes in reproductive physiology of lactating dairy cows due to elevated steroid metabolism[END_REF]. However, it is hard to understand which underlying mechanisms cause this decline in fertility. The regulation of estrus is controlled by the interplay of various organs and hormones. Mathematical modeling of the involved mechanisms is expected to improve insight in the biological processes underlying the bovine estrous cycle, and could thereby help to find causes of declined fertility in dairy cows [START_REF] Boer | Estrous behavior in dairy cows: identification of underlying mechanisms and gene functions[END_REF].

Although the endocrine and physiologic regulation of the bovine estrous cycle is studied extensively, mathematical models of cycle regulation are scarce and of limited scope [START_REF] Meier | A compartmental model desribing changes in progesterone concentrations during the estrous cycle[END_REF][START_REF] Soboleva | A model of follicular development and ovulation in sheep and cattle[END_REF]. A number of models have been developed for other ruminant species, especially ewes [START_REF] Clément | Contribution of a mathematical modelling approach to the understanding of the ovarian function[END_REF][START_REF] Heinze | A mathematical model of luteinizing hormone release from ovine pituitary cells in perifusion[END_REF], but these models do not contain all the key players that are required to simulate follicle development and the accompanying hormone levels throughout consecutive cycles.

A model that integrates the major tissues and hormones involved, and that is able to simulate the dynamics of follicular development, has been developed for the human menstrual cycle by Reinecke [START_REF] Reinecke | A complex mathematical model of the human menstrual cycle[END_REF]. This model, which is based on previous work by Selgrade and colleagues [START_REF] Selgrade | A model for the production of ovarian hormones during the menstrual cycle[END_REF][START_REF] Schlosser | A model of gonadotropin regulation during the menstrual cycle in women: Qualitative features[END_REF][START_REF] Clark | Multiple stable periodic solutions in a model for hormonal control of the menstrual cycle[END_REF], describes the dynamics of hormones, enzymes, receptors, and follicular phases throughout the cycle in a set of differential equations.

The objective of the work described in this paper was to develop a mathematical model of the dynamics of the bovine estrous cycle on individual cow level, that is able to simulate follicle development and the accompanying fluctuations in hormone concentrations. Physiologic and endocrine mechanisms that regulate the cycle are very similar between human and cows. There-fore, some mechanisms of the human model in [START_REF] Reinecke | A complex mathematical model of the human menstrual cycle[END_REF] could be used (although sometimes with simplifications), and extended with other mechanisms like follicular wave emergence and corpus luteum regression.

Focus in this paper is on the model development. Section 2 describes the biological mechanisms of the bovine estrous cycle and how these mechanisms are incorporated in the model. In Section 3, the mathematical description and all model equations and parameters are given. Simulation results are presented in Section 4, showing that a set of equations and parameters is obtained that describes the system consistent with biological data for cows.

In Section 5, it is discussed how the current model could be applied and extended.

Biological background

Follicles

Two different patterns of follicle development are identified in mammals.

In humans (and rats and pigs), the development of follicles to ovulatory size occurs only during the follicular phase, while in cattle (and sheep and horses), development of follicles to ovulatory or near-ovulatory size occurs throughout the cycle [START_REF] Fortune | Ovarian follicular growth and development in mammals[END_REF]. A normal cycle includes two or three wave-like patterns of follicle development, in which a cohort of follicles start to grow.

The average duration of the bovine estrous cycle is 20 days for 2-wave and 22 days for 3-wave cycles (reviewed in [START_REF] Adams | Progress in understanding ovarian follicular dynamics in cattle[END_REF]). Each follicular wave is initiated by an increase of follicle stimulating hormone (FSH) release from the anterior pituitary [START_REF] Ginther | Role of low circulating FSH concentrations in controlling the interval to emergence of the subsequent follicular wave in cattle[END_REF]. The growing follicles produce estradiol (E2) and inhibin (Inh), which are released into peripheral blood. In the first one or two waves, a dominant follicle deviates from the cohort of growing follicles that does not ovulate, but undergoes regression under influence of progesterone (P4) produced by the corpus luteum (CL). When the CL is regressed under influence of PGF2α, the concentration of P4 decreases [START_REF] Niswender | Mechanisms controlling the function and life span of the corpus luteum[END_REF]. The dominant follicle present at that moment develops and matures, and ovulation can then take place because the inhibiting effect of P4 on the surge of luteinizing hormone (LH) is removed [START_REF] Mann | Progesterone inhibition of the development of the luteolytic signal in cows[END_REF]. Elevated E2 levels increase the secretion of gonadotropin releasing hormone (GnRH), which triggers the LH surge and thereby induces ovulation. Once an oocyte is successfully ovulated, the remains of the follicle form a new P4-producing CL. If conception has failed, the CL regresses, P4 levels decrease, and the cycle restarts (reviewed in [START_REF] Boer | Estrous behavior in dairy cows: identification of underlying mechanisms and gene functions[END_REF]).

The ovaries contain a pool of small follicles with immature oocytes. Under influence of FSH, a cohort of 8-41 growing follicles emerge [START_REF] Adams | Progress in understanding ovarian follicular dynamics in cattle[END_REF]. Approximately two days after cohort recruitment, one follicle is selected to become the dominant follicle, and continues to grow [START_REF] Bao | Expression of steroidogenic enzyme and gonadotropin receptor genes in bovine follicles during ovarian follicular waves: A review[END_REF]. This deviation of the dominant follicle is associated with increased FSH and LH receptor binding, activating the enzymes that catalyze steroidogenesis, resulting in increased E2 production and higher E2 serum levels [START_REF] Bao | Expression of steroidogenic enzyme and gonadotropin receptor genes in bovine follicles during ovarian follicular waves: A review[END_REF]. The dominant follicle expresses more FSH receptors, and it can therefore continue to grow even when FSH serum levels are low [START_REF] Beg | Follicle selection in cattle: Dynamics of follicular fluid factors during development of follicle dominance[END_REF]. In the model, the emergence of a follicular wave is induced when FSH exceeds a threshold which becomes lower when follicles become larger, representing that larger follicles are more sensitive to FSH. Dominant follicles also secrete increasing amounts of inhibin (Inh). Inh suppresses FSH and, hence, suppresses the growth of subordinate follicles.

Ovulation or regression of the dominant follicle eliminates this suppression, allowing the onset of the next follicular wave [START_REF] Bleach | Plasma inhibin A in heifers: Relationship with follicle dynamics, gonadotropins, and steroids during the estrous cycle and after treatment with bovine follicular fluid[END_REF][START_REF] Ginther | Follicle selection in monovular species[END_REF].

Small follicles of an emerging cohort may release very small amounts of E2 and Inh per follicle, but taken together, this amount is not negligible. Furthermore, there is always a medium-size or large follicle present [START_REF] Ireland | Historical perspective of turnover of dominant follicles during the bovine estrous cycle: Key concepts, studies, advancements, and terms[END_REF][START_REF] Wise | Biochemical analysis of bovine follicular fluid: albumin, total protein, lysosomal enzymes, ions, steroids and ascorbic acid content in relation to follicular size, rank, atresia classification and day of estrous cycle[END_REF][START_REF] Wolfenson | Follicular dynamics and concentrations of steroids and gonadotropins in lactating cows and nulliparous heifers[END_REF], which results in a basal hormone production throughout the cycle. Different follicles are recruited, growing, and regressing in each cycle and in each wave. However, total E2 and Inh production capacity is modeled as a continuous function throughout subsequent waves and cycles, representing the total amount of hormone production of the follicles present at any moment. Follicle regression is promoted by high P4 levels and by the LH surge (Equation 7). The capacity of follicles to produce E2 and Inh is denoted as "follicular function" in the rest of this paper.

Corpus luteum

The CL develops within 2-3 days after ovulation, starting the synthesis and release of P4, which maintains the readiness of the endometrium for receiving the embryo. In absence of a conceptus, the CL will regress at day 17-18 of the cycle [START_REF] Miyamoto | Local regulation of corpus luteum development and regression in the cow: Impact of angiogenic and vasoactive factors[END_REF][START_REF] Taylor | Follicular dynamics and corpus luteum growth and function in pregnant versus nonpregnant cows[END_REF]. In each cycle a new CL develops, but CL development is modeled as a continuous function of P4 producing tissue, denoted as "CL function" in the rest of this paper. In the model, CL development is induced by the LH surge. A threshold and delay are incorporated in the effect of LH on the CL, to account for the time required for the process of transition from follicle to CL [START_REF] Niswender | Mechanisms controlling the function and life span of the corpus luteum[END_REF] and the shift from E2 to P4 production [START_REF] Dieleman | Peripheral plasma concentrations of oestradiol, progesterone, cortisol, LH and prolactin during the oestrous cycle in the cow, with emphasis on the peri-oestrous period[END_REF][START_REF] Dieleman | Progesterone-synthesizing ability of preovulatory follicles of cows relative to the peak of LH[END_REF]. If the CL reaches a certain size, it continues to grow without further stimulation by LH [START_REF] Skarzynski | Luteotropic mechanisms in the bovine corpus luteum: Role of oxytocin, prostaglandin F2α, progesterone and noradrenaline[END_REF]. CL regression is induced by PGF2α secretion from the uterus (described in Section 2.4). Growth and regression of CL function are described by Equation 9.

Estradiol and inhibin

E2 affects LH synthesis and release [START_REF] Glidewell-Kenney | Nonclassical estrogen receptor α signaling mediates negative feedback in the female mouse reproductive axis[END_REF] and FSH release [START_REF] Beg | Follicle selection in cattle: Dynamics of follicular fluid factors during development of follicle dominance[END_REF][START_REF] Lane | Alterations in the ability of the bovine pituitary gland to secrete gonadotropins in vitro during the first follicle-stimulating hormone increase of the estrous cycle and in response to exogenous steroids[END_REF]. E2 serum levels are higher in ovulatory than in non-ovulatory waves [START_REF] Bleach | Plasma inhibin A in heifers: Relationship with follicle dynamics, gonadotropins, and steroids during the estrous cycle and after treatment with bovine follicular fluid[END_REF][START_REF] Echternkamp | Concurrent changes in bovine plasma hormone levels prior to and during the first postpartum estrous cycle[END_REF] and reach peak levels around estrus [START_REF] Bleach | Plasma inhibin A in heifers: Relationship with follicle dynamics, gonadotropins, and steroids during the estrous cycle and after treatment with bovine follicular fluid[END_REF][START_REF] Echternkamp | Concurrent changes in bovine plasma hormone levels prior to and during the first postpartum estrous cycle[END_REF][START_REF] Evans | Changes in androgen secretion and luteinizing hormone pulse amplitude are associated with the recruitment and growth of ovarian follicles during the luteal phase of the bovine estrous cycle[END_REF][START_REF] Glencross | Relationships between the incidence of pre-ovulatory behaviour and the concentrations of oestradiol-17β and progesterone in bovine plasma[END_REF][START_REF] Kanchev | Concentration of steroids in bovine peripheral plasma during the oestrous cycle and the effect of betamethasone treatment[END_REF][START_REF] Wise | Biochemical analysis of bovine follicular fluid: albumin, total protein, lysosomal enzymes, ions, steroids and ascorbic acid content in relation to follicular size, rank, atresia classification and day of estrous cycle[END_REF]. This suggests that the preovulatory follicle has the largest capacity to produce and release E2, although its maximum size is not significantly different from the maximum size of non-ovulatory dominant follicles. Considering the results in [START_REF] Acosta | Studies of follicular vascularity associated with follicle selection and ovulation in cattle[END_REF][START_REF] Acosta | Periovulatory changes in the local release of vasoactive peptides, prostaglandin F2α, and steroid hormones from bovine mature follicles in vivo[END_REF],

where a better vascularity of the ovulatory follicle is reported, it is reasonable that the ovulatory follicle can secrete more E2 than non-ovulatory follicles and, consequently, E2 serum levels are highest at estrus. In the model, the rate of E2 production and release to the blood is taken as proportional to follicular function (Equation 11).

Inh inhibits FSH synthesis and thus reduces FSH release [START_REF] Ginther | Follicle selection in monovular species[END_REF]. Compared to basal Inh serum levels, peak levels are almost doubled in non-ovulatory waves and increase further in ovulatory waves [START_REF] Parker | Plasma concentrations of inhibin A and follicle-stimulating hormone differ between cows with two or three waves of ovarian follicular development in a single estrous cycle[END_REF]. There are different forms of inhibin, but only inhibin A is considered in the model, as it is the predominant form in bovine follicular fluid [START_REF] Beg | Follicle selection in cattle: Dynamics of follicular fluid factors during development of follicle dominance[END_REF]. In the model, Inh production rate is taken as proportional to follicular function (Equation 12).

Progesterone and prostaglandin F2α

The CL is the main source of P4. Serum P4 concentration is near to zero around estrus and high during the luteal phase [START_REF] Adams | Effect of progesterone on ovarian follicles, emergence of follicular waves and circulating follicle-stimulating hormone in heifers[END_REF][START_REF] Díaz | Plasma progesterone levels during the estrous cycle of Holstein and Brahman cows, Carora type and cross-bred heifers[END_REF][START_REF] Echternkamp | Concurrent changes in bovine plasma hormone levels prior to and during the first postpartum estrous cycle[END_REF][START_REF] Kaneko | Changes in plasma concentrations of immunoreactive inhibin, estradiol and FSH associated with follicular waves during the estrous cycle of the cow[END_REF][START_REF] Stabenfeldt | Peripheral plasma progesterone levels during the bovine oestrous cycle[END_REF]. A high correlation between CL diameter and P4 output was reported in [START_REF] Perry | Endocrine changes and ultrasonography of ovaries in suckled beef cows during resumption of postpartum estrous cycles[END_REF][START_REF] Savio | Pattern of growth of dominant follicles during the oestrous cycle of heifers[END_REF][START_REF] Wolfenson | Follicular dynamics and concentrations of steroids and gonadotropins in lactating cows and nulliparous heifers[END_REF].

In the model, the rate of P4 release into the blood is taken as proportional to CL function (Equation 10). Pulsatile PGF2α release from the uterus induces CL regression. The rise of P4 early in the cycle initiates a series of events or mechanisms that eventually lead to the rise of PGF2α, followed by a decline of PGF2α a few days later. It was shown that administration of P4 prior to its natural rise resulted in an equally earlier onset of CL regression [START_REF] Mann | Role of early luteal phase progesterone in control of the timing of the luteolytic signal in cows[END_REF]. Exposure to effective amounts of P4 must last for 10-13 days to induce PGF2α pulses [START_REF] Mann | Role of early luteal phase progesterone in control of the timing of the luteolytic signal in cows[END_REF][START_REF] Poyser | The control of prostaglandin production by the endometrium in relation to luteolysis and menstruation[END_REF][START_REF] Mccracken | Luteolysis: A neuroendocrine-mediated event[END_REF][START_REF] Goff | Steroid hormone modulation of prostaglandin secretion in the ruminant endometrium during the estrous cycle[END_REF]. Peak PGF2α levels are 3-4 times higher than basal levels [START_REF] Santos | Elevated progesterone concentrations enhance prostaglandin F2α synthesis in dairy cows[END_REF][START_REF] Mann | Relationship between maternal endocrine environment, early embryo development and inhibition of the luteolytic mechanism in cows[END_REF][START_REF] Shaham-Albalancy | Delayed effect of low progesterone concentrations on bovine uterine PGF2α secretion in the subsequent oestrous cycle[END_REF][START_REF] Araujo | Role of follicular estradiol-17beta in timing of luteolysis in heifers[END_REF].

PGF2α is regulated by oxytocin (OT), P4 and E2 [START_REF] Silvia | Review: Hormonal regulation of uterine secretion of prostaglandin F2α during luteolysis in ruminants[END_REF]. P4 first prevents a too early release of PGF2α pulses, but simultaneously stimulates synthesis of enzymes required for PGF2α production. In the later luteal phase, changed expression of P4 and OT receptors results in a gradual decrease in the suppression of PGF2α [START_REF] Santos | Elevated progesterone concentrations enhance prostaglandin F2α synthesis in dairy cows[END_REF], leading to an OT induced pulsatile release of PGF2α [START_REF] Araujo | Role of follicular estradiol-17beta in timing of luteolysis in heifers[END_REF][START_REF] Poyser | The control of prostaglandin production by the endometrium in relation to luteolysis and menstruation[END_REF]. How these mechanisms are regulating each other is quite complex and not understood in full detail.

What is clear is that the rise in P4 levels and the continued presence of P4 above an effective level sets in motion a series of events that lead to CL regression. Hence, we incorporated these series of events as a black box using time delays to obtain the right timing of PGF2α signaling. In the model, PGF2α increases a specific number of days (delay τ P 4,1 ) after P4 levels reach a threshold. Similarly, PGF2α declines another (larger) number of days (delay τ P 4,2 ) after P4 levels reached a threshold (Equation 8).

Gonadotropin releasing hormone, luteinizing hormone and follicle stimulating hormone

Pulsatile signaling of GnRH regulates LH and FSH secretion [START_REF] Pawson | The pituitary effects of GnRH[END_REF]. Because GnRH induces the LH surge, it indirectly induces ovulation [START_REF] Troxel | The effect of progestin and GnRH treatments on ovarian function and reproductive hormone secretions of anestrous postpartum suckled beef cows[END_REF]. The GnRH pulse generator is located in the hypothalamus and is modulated by P4 and E2 [START_REF] Goodman | The physiology of reproduction[END_REF]. During the luteal phase, both P4 and E2 suppress the activity of the GnRH pulse generator. During pro-estrus however, elevated E2 levels change estrogen receptor signaling, which induces a GnRH surge [START_REF] Glidewell-Kenney | Nonclassical estrogen receptor α signaling mediates negative feedback in the female mouse reproductive axis[END_REF][START_REF] Goodman | The physiology of reproduction[END_REF]. GnRH is released into the portal circulation of the pituitary and binds to GnRH receptors of the anterior pituitary [START_REF] Vizcarra | Effect of gonadotropin-releasing hormone (GnRH) pulse frequency on serum and pituitary concentrations of luteinizing hormone and folliclestimulating hormone, GnRH receptors, and messenger ribonucleic acid for gonadotropin subunits in cows[END_REF]. In the model, GnRH stimulates LH release, resulting in an LH surge concurrently with the GnRH surge.

GnRH synthesis is taken constant as long as the amount of GnRH in the hypothalamus is below a threshold (Equation 1). GnRH release is inhibited when P4 levels are above a threshold and when both P4 and E2 levels are above a threshold. GnRH release is stimulated when P4 levels are low and E2 reaches a threshold (Equation 1b), resulting in a surge of GnRH. GnRH concentration in the pituitary depends on GnRH amount released from the hypothalamus, and is further increased by high E2 levels, representing that E2 up-regulates expression of GnRH receptors [START_REF] Goodman | The physiology of reproduction[END_REF][START_REF] Vizcarra | Effect of gonadotropin-releasing hormone (GnRH) pulse frequency on serum and pituitary concentrations of luteinizing hormone and folliclestimulating hormone, GnRH receptors, and messenger ribonucleic acid for gonadotropin subunits in cows[END_REF] (Equation 2).

The LH surge at the day before ovulation induces ovulation of the ovulatory follicle and formation of the CL. The LH surge will shut down E2 and Inh production capacity of the ovulatory follicle [START_REF] Chenault | Transitory changes in plasma progestins, estradiol, and luteinizing hormone approaching ovulation in the bovine[END_REF][START_REF] Wolfenson | Follicular dynamics and concentrations of steroids and gonadotropins in lactating cows and nulliparous heifers[END_REF]. High P4 levels suppress the release of LH via the inhibition of the GnRH pulse generator [START_REF] Bergfeld | Changing dose of progesterone results in sudden changes in frequency of luteinizing hormone pulses and secretion of 17β-estradiol in bovine females[END_REF]. Additionally, high P4 levels decrease pituitary sensitivity to E2, thereby increasing the amount of E2 required to induce an LH surge above physio-logical levels [START_REF] Goodman | The physiology of reproduction[END_REF]. Peak LH levels are about five times as high as basal levels or higher [START_REF] Bleach | Plasma inhibin A in heifers: Relationship with follicle dynamics, gonadotropins, and steroids during the estrous cycle and after treatment with bovine follicular fluid[END_REF][START_REF] Echternkamp | Concurrent changes in bovine plasma hormone levels prior to and during the first postpartum estrous cycle[END_REF][START_REF] Kotwica | Relationship of plasma testosterone concentrations to pituitary-ovarian hormone secretion during the bovine estrous cycle and the effects of testosterone propionate administered during luteal regression[END_REF][START_REF] Dieleman | Peripheral plasma concentrations of oestradiol, progesterone, cortisol, LH and prolactin during the oestrous cycle in the cow, with emphasis on the peri-oestrous period[END_REF]]. In the model, LH synthesis is stimulated by E2 and inhibited by P4 (Equation 5a). Besides a small basal LH release, there is a surge of LH when GnRH in the pituitary reaches a threshold (Equation 5b).

FSH synthesis is inhibited by Inh [START_REF] Beg | Follicle selection in cattle: Dynamics of follicular fluid factors during development of follicle dominance[END_REF]. P4 and E2 modulate FSH release via effects on the anterior pituitary and on the GnRH pulse generator in the hypothalamus. Peak FSH serum levels are about three times higher than basal levels [START_REF] Bleach | Plasma inhibin A in heifers: Relationship with follicle dynamics, gonadotropins, and steroids during the estrous cycle and after treatment with bovine follicular fluid[END_REF][START_REF] Evans | Changes in androgen secretion and luteinizing hormone pulse amplitude are associated with the recruitment and growth of ovarian follicles during the luteal phase of the bovine estrous cycle[END_REF]. In the model, FSH synthesis in the pituitary is increased when Inh levels are below a threshold (Equation 3a). FSH release from the pituitary to the blood is stimulated by P4 and GnRH, and inhibited by E2 (Equation 3b).

Mathematical formulation

The mathematical approach used for the bovine model is comparable to the approach used for the model of the human menstrual cycle, which originally has been developed at North Carolina State University by Selgrade and colleagues [START_REF] Selgrade | A model for the production of ovarian hormones during the menstrual cycle[END_REF][START_REF] Schlosser | A model of gonadotropin regulation during the menstrual cycle in women: Qualitative features[END_REF][START_REF] Harris | Differential equation models for the hormonal regulation of the menstrual cycle[END_REF][START_REF] Clark | Multiple stable periodic solutions in a model for hormonal control of the menstrual cycle[END_REF], and has been extended at the Zuse Institute [START_REF] Reinecke | A complex mathematical model of the human menstrual cycle[END_REF][START_REF] Reinecke | Mathematical modeling and simulation of the female menstrual cycle[END_REF].

The system is considered in four compartments: hypothalamus, anterior pituitary, ovaries and uterus, connected through peripheral and portal blood (Figure 1). The model includes the processes of follicle and CL development and the key hormones that interact to control these processes as described in Section 2. The gonadotropin equations are based on synthesis-releaseclearance relations. This structure was first introduced in [START_REF] Schlosser | A model of gonadotropin regulation during the menstrual cycle in women: Qualitative features[END_REF] Based on these mechanisms, 12 ordinary differential equations (ODEs) with 54 parameters are formulated. If necessary, time delays are incorporated to model the time between events and their effects, representing the duration of intermediate steps in biological processes. In this case, the ODE is turned into a delay differential equation (DDE). To solve the system of differential equations, we use the solver RADAR5 [START_REF] Guglielmi | [END_REF], which has been designed for the solution of stiff delay differential equations.

Hill functions

Because the exact mechanisms are often not known or more specific than necessary, Hill functions are used to model stimulatory and inhibitory effects of the hormones. They are used whenever there is a nonlinear relation between two substances. A Hill function is a sigmoidal function between zero and one, which switches at a specified threshold from one level to the other with a specified steepness. Positive Hill functions are used for stimulating effects and are defined as

h + (S(t); T, n) := S(t) n T n + S(t) n .
S(t) represents the effector, T the threshold for change of behavior, and n controls the steepness of the curve. Negative Hill functions are used for inhibitory effects and are defined as

h -(S(t); T, n) := T n T n + S(t) n .
Here, the value of the function has its maximum at the lowest value of the initiating substrate S(t), and switches to zero if this substrate passes the threshold T . Whenever a Hill function is used, it is provided with another parameter m that controls the height of the switch. This parameter serves as maximum stimulatory respectively inhibitory effect. For abbreviation of notation, we use H + (S) instead of m • h + (S(t); T, n). We usually choose the steepness coefficient n = 2, but, when appropriate, we set n = 1, 5, or 10 to capture

Model equations

The amount of GnRH in the hypothalamus is a result of synthesis in the hypothalamus and release into the pituitary,

d dt GnRH Hypo (t) = Syn GnRH (t) -Rel GnRH (t). (1)
GnRH synthesis depends on its current level in the hypothalamus. If this level approaches a specified threshold, synthesis decreases to zero. This effect is modeled as

Syn GnRH (t) = c GnRH ,1 • 1 - GnRH Hypo (t) GnRH max Hypo . ( 1a
)
As long as GnRH is far below its maximum, the factor 1 -

GnRH Hypo (t) GnRH max
Hypo has only a small impact. The release of GnRH from the hypothalamus to the pituitary is dependent on its current level in the hypothalamus. E2 inhibits GnRH release during the luteal phase, i.e. if P4 and E2 are high at the same time, described by H - 1 (P4&E2). H - 1 (P4&E2) denotes the sum of two Hill functions minus their product, and inhibits GnRH release only if both substrates are above their threshold. Additionally, the release of GnRH is inhibited by P4 only,

Rel GnRH (t) = (H - 1 (P4&E2) + H - 2 (P4)) • GnRH Hypo (t). ( 1b
)
Changes in GnRH amount in the pituitary are dependent on the released amount from the hypothalamus, but also on the presence of E2. E2 increases the number of GnRH receptors in the pituitary. This effect is included in the equation as a positive Hill function. GnRH clearance from pituitary portal blood is proportional to the GnRH level in the pituitary, i.e. GnRH clearance

is represented by c GnRH,2 • GnRH P it (t), in which c GnRH,2 is a constant, d dt GnRH Pit (t) = Rel GnRH (t) • H + 3 (E2) -c GnRH ,2 • GnRH Pit (t). (2)
FSH is synthesized in the pituitary and released into the blood,

d dt FSH Pit (t) = Syn FSH (t) -Rel FSH (t). (3)
FSH synthesis rate in the pituitary is only dependent on delayed Inh, as in [START_REF] Harris | Differential equation models for the hormonal regulation of the menstrual cycle[END_REF]. FSH is synthesized when the Inh level is low, i.e. high Inh levels inhibit FSH synthesis, which is included as a negative Hill function,

Syn FSH = H - 4 (Inh τ ). ( 3a
)
The index τ stands for a delayed effect of Inh, i.e. Inh is considered at time tτ . FSH release from the pituitary to the blood is stimulated by P4 and GnRH, and inhibited by E2,

Rel FSH = (H + 5 (P4) + H - 6 (E2) + H + 7 (GnRH Pit )) • FSH Pit (t). (3b)
Concluding, FSH serum level is a result of the difference between the released amount from the pituitary and clearance in the blood,

d dt FSH Blood (t) = Rel FSH (t) -c FSH • FSH Blood (t), (4)
where c FSH is the FSH clearance rate constant.

Like FSH, the LH serum level depends on synthesis in the pituitary, release into the blood and clearance thereof,

d dt LH Pit (t) = Syn LH (t) -Rel LH (t). ( 5
)
LH synthesis in the pituitary is stimulated by E2 and inhibited by P4,

Syn LH (t) = H + 8 (E2) + H - 9 (P4). (5a)
We assume a low constant basal LH release b LH from the pituitary into the blood. On top of that, LH release is stimulated by GnRH,

Rel LH (t) = (b LH + H + 10 (GnRH Pit )) • LH Pit (t). ( 5b
)
Summarizing, LH in the blood is obtained as

d dt LH Blood (t) = Rel LH (t) -c LH • LH Blood (t), (6)
where c LH is the LH clearance rate constant.

Follicular function is stimulated by FSH, whereas its decrease is promoted by P4 and the LH surge,

d dt Foll (t) = H + 11 (FSH ) -(H + 12 (P4) + H + 13 (LH Blood )) • Foll (t). ( 7
)
The sensitivity of the follicles to respond to FSH grows with their size. In the model, the threshold of FSH to stimulate the follicular function decreases with increasing follicular function. For this effect of a rising FSH sensitivity, a negative Hill function is included to control the threshold of FSH,

T Foll FSH (t) := T Foll FSH • h -(Foll (t); T FSH Foll , 1),
and the Hill function for the effect of FSH on follicular function becomes

H + 11 (FSH ) := m Foll FSH • h + (FSH Blood (t); T Foll FSH (t), 2). ( 7a
)
PGF2α initiates the functional regression of the CL, and thereby the decrease in P4 levels. After a large time delay, PGF2α synthesis is stimulated by elevated P4 levels above a specified threshold value. The PGF2α level declines a couple of days after its rise, which is included as a delayed positive effect of P4 on the decay of PGF2α,

d dt PGF2 α(t) = H + 14 (P4 τ 1 ) -H + 15 (P4 τ 2 ) • PGF2 α(t). ( 8
)
The LH peak initiates growth of the CL with a specified delay. After reaching a certain size, the CL continues to grow on its own as long as PGF2α is low. The CL starts to regress when PGF2α levels rise above a threshold,

d dt CL(t) = H + 16 (LH τ ) + H + 17 (CL) -H + 18 (PGF2 α) • CL(t). ( 9
)
The production of P4 in the ovary is assumed to be proportional to CL function, and the production of E2 and Inh is assumed to be proportional to follicular function. Therefore, the equations for P4, E2, and Inh do not contain any Hill functions,

d dt P4(t) = c P4 CL • CL(t) -c P4 • P4(t), ( 10
)
d dt E2(t) = c E2 Foll • Foll (t) -c E2 • E2(t), ( 11
)
d dt Inh(t) = c Inh Foll • Foll (t) -c Inh • Inh(t). ( 12
)
The parameters c P4 , c E2 and c Inh denote the respective clearance rate constants.

Figure 2 gives an overview of all mechanisms described by the model equations. Detailed notations for the Hill functions, parameters, and equations are given in Appendix A, Appendix B, and Appendix C respectively.

Parameter identification and sensitivity analysis

The main difficulty is not to simulate the system, i.e. to solve the differential equations, but to identify the unknown parameters. Unfortunately, many of the parameters are not measurable. Sometimes the range of values is known, but some parameters are completely unknown. The techniques for parameter estimation that are used in this model are implemented in the software packages PARKIN [START_REF] Nowak | Numerical identification of selected rate constants in large chemical reaction systems[END_REF][START_REF] Deuflhard | Efficient numerical simulation and identification of large chemical reaction systems[END_REF] and NLSCON [START_REF] Nowak | Nonlinear Least Squares with nonlinear equality CONstraints[END_REF], which have been developed at the Zuse Institute for many years. These programs take into account parameter sensitivities and linear dependencies, and include a number of optimization methods such as, for example, affine covariant Gauss-Newton methods [START_REF] Deuflhard | Newton Methods for Nonlinear Problems: Affine Invariance and Adaptive Algorithms[END_REF]. A renewed version of this software, especially adapted to parameter identification in ordinary differential equation models, has been used throughout the paper. The mathematical background is described in [START_REF] Deuflhard | Newton Methods for Nonlinear Problems: Affine Invariance and Adaptive Algorithms[END_REF].

To obtain a good initial guess for the parameter optimization procedure, we use a model decomposition approach and successively enlarge the set of estimated parameters. The first step is to define input curves representing the development of Inh, P4, and E2 levels in the blood over time. This use of explicit functions, which simplifies parameter identification, was already suggested by Schlosser [START_REF] Schlosser | A model of gonadotropin regulation during the menstrual cycle in women: Qualitative features[END_REF]. Composition of these input curves is based on published data for endocrine profiles of cows with a normal estrous cycle, see for example [START_REF] Perry | The bovine estrous cycle[END_REF]. Following the approach in [START_REF] Harris | Differential equation models for the hormonal regulation of the menstrual cycle[END_REF], we use the input curves to successively fit the profiles of the other components. The detailed procedure can be found in [START_REF] Boer | A simple mathematical model of the bovine estrous cycle: follicle development and endocrine interactions[END_REF]. In the last step, the input curves for P4, E2, and Inh are replaced by their original ODE/DDE description to obtain a closed network. The final parameter values are listed in Table B.1, and the corresponding simulation results are illustrated in Figure 4.

A sensitivity analysis has been performed with the techniques desribed in [START_REF] Deuflhard | Newton Methods for Nonlinear Problems: Affine Invariance and Adaptive Algorithms[END_REF]. A more detailed description including column norms of the sensitivity matrix and subconditions, which provide information about the sensitivities and the dependencies of the parameters, can be found in [START_REF] Boer | A simple mathematical model of the bovine estrous cycle: follicle development and endocrine interactions[END_REF]. It turns out that among the most sensitive and best predictable parameters are p 36 = τ P 4,1 , p 11 = τ Inh , p 20 = c F SH , and p 39 = τ P 4,2 .

Simulation results

The figures in this section show the computed dynamics of follicle and CL development and accompanying fluctuations in hormone levels over consecutive cycles. The simulation results show that the current set of model parameters generates curves consistent with empirical knowledge for cows with a normal estrous cycle with three follicular waves. Notice that the model generates consecutive cycles that are not entirely identical (quasi-periodic behavior), but that vary slightly in patterns and peak heights between cycles.

Small differences in model output at the end of a cycle result in a different

starting point of the next cycle, which leads to variation between the curves. This variation in hormone levels between cycles could well resemble variation within a cow over consecutive cycles. However, a different parameterization can be used to produce a stable limit cycle.

Each estrous cycle contains three waves of follicular growth (Figure 5). The CL starts to grow a few days after ovulation and is large during the first two follicular waves, which suppresses follicle growth. As the larger follicles become more sensitive to P4, at a certain size the effect of P4 becomes so large that it induces follicle regression. After regression of the CL, the dominant follicle of the third follicular wave can continue to grow, leading to ovulation, which causes a sharp decline in follicular function.

The pattern of serum E2 levels is a result of follicular function (Figures 5 and6). The third wave of follicular growth takes place when P4 levels Increased FSH levels induce the growth of a follicular wave and thereby the start of Inh increase, but FSH is suppressed when Inh levels are above a certain level (Figure 7). Notice that FSH peak levels in the third wave of the cycle differ in consecutive cycles because of corresponding differences in height of the GnRH surge (Figures 6 and7). When Inh has declined due to follicular regression, FSH increases again and induces the next follicular wave. Because follicular growth is modeled in three waves, also Inh levels rise in three waves in a cycle. P4 serum levels are proportional to CL function. P4 concentration is small during the first days of the cycle and rises when the CL starts to grow (Figure 5). Notice that a lower LH peak height results in a less steep P4 increase and lower levels of P4 in the following cycle (Figures 6 and8).

Increased P4 levels induce a rise in PGF2α after a couple of days, which causes CL regression and declining P4.

Discussion and Outlook

The current mathematical model describes the interaction between a number of key physiological processes of the bovine estrous cycle. The model is able to simulate the dynamics of follicle and CL growth and development, as well as the associated hormone level changes in consecutive cycles. The current model comprises 12 equations and 54 parameters. The estrous cycles generated by the model are not entirely identical and could well resemble variations within a cow over consecutive cycles.

The above simulations show a quasi-periodic behavior, but a different parameterization (not listed in this paper) could be used to produce a stable limit cycle. This shows that the variations between simulated cycles are not an intrinsic characteristic of the model, but depend on the parameterization. However, the cycles of a real cow are usually quite irregular, and we think this is not due to changes in external factors for that cow but rather arises from the fact that each cycle presents slightly new and somewhat different 'starting values' for the next cycle, which we think that our model mimics.

Alternatively, one could add a stochastic component to the regular system (representing small variations in external factors) to induce variations in consecutive cycles, but this was not in the scope of our work. The modeling method with ODEs/DDEs as used for the presented model of the bovine estrous cycle was also used for the model of the human menstrual cycle [START_REF] Selgrade | A model for the production of ovarian hormones during the menstrual cycle[END_REF][START_REF] Schlosser | A model of gonadotropin regulation during the menstrual cycle in women: Qualitative features[END_REF][START_REF] Harris | Differential equation models for the hormonal regulation of the menstrual cycle[END_REF][START_REF] Reinecke | A complex mathematical model of the human menstrual cycle[END_REF]. As we aimed at the development of a model for the dynamical changes of a biological system, including the information about how components influence the rates of change of other components, our approach to model the system with differential equations appears to be the most reasonable. Maybe qualitative results could have been obtained

with other methods such as, for example, boolean networks, but differential equations allow for a simulation of quantitative profiles of the involved components. To our knowledge, no comparable models of the bovine estrous cycle are available.

The current model describes the mechanisms of an idealized cow, based on average numbers obtained from several data sources. It would in principle be possible to fit the model to measurement data of an individual cow that would show small deviations of the cycle, or even a pathological abnormal cycle due to certain disorders. This would represent the next step in the modeling approach. Because empirical data are usually noisy, parameter optimization would then also have to take measurement errors into account.

Although the current model could thus offer possibilities to simulate fertility disorders, its predictive ability may be limited in those parts and for those aspects in which the model is not entirely mechanistic but rather descriptive. One example thereof is the modeling of PGF2α. Because the detailed biological mechanisms that induce the rise of PGF2α are very complex and not completely understood, we chose to restrict the number of state variables for this part of the model, and to include time delays. This mimics the situation in cows that the rise of P4 early in the cycle starts a series of events or mechanisms that eventually lead to the rise of PGF2α, followed by a decline of PGF2α several days later. The time delays are thus a 'black box' where the intermediate events that regulate PGF2α levels are not described.

In this way, we were able to obtain the right time point of CL regression even though we don't know the biological mechanisms exactly. By reducing the delays, the duration of the luteal phase can be reduced. This could mean that P4 serum levels already decline during the second wave of follicle development, which could then become the ovulatory wave. The shorter delays could thus result in a shorter cycle with only two follicular waves. However, the consequence of the chosen approach is that the predictive abilities for this part of the model are limited. Model improvement and refinement of this sub-model will play an important role in future work.

Apart from fitting of the model to individual cow data, mentioned above, we plan to use this model to determine the level of control exerted by various system components on the functioning of the system. Examples of such model applications are to explore the mechanisms that influence the pattern of follicular waves, or to study hormone patterns associated with subfertility. Also, the model can serve as a basis for more elaborate models and simulations, with the ability to study effects of external manipulations and genetic differences. Possible extensions of the model could be in the field of energy metabolism, stress, disease, and factors affecting the expression of estrous behavior. There are relationships between regulation of the estrous cycle and energy balance, which can cause fertility problems in high producing dairy cows in negative energy balance (for reviews see [START_REF] Diskin | Effects of nutrition and metabolic status on circulating hormones and ovarian follicle development in cattle[END_REF][START_REF] Roche | The effect of nutritional management of the dairy cow on reproductive efficiency[END_REF]). Changes in reproductive performance that are associated with high milk production may in part be explained by elevated P4 and E2 clearance rates, as described in the physiological model of [START_REF] Wiltbank | Changes in reproductive physiology of lactating dairy cows due to elevated steroid metabolism[END_REF]. In this physiological model, clearance rates of hormones by the liver of cows with high milk production are increased as a result of elevated feed intake, leading to an increased liver blood flow and metabolic activity. With a similar level of hormone production, circulating hormone levels would thus be lower. Lameness, an example of a stress inducing condition, was found to inhibit the LH surge and ovulation, whereas incidence of estrous behavior (although with less intensity) was not reduced. These observations suggest that stress, caused by lameness, reduces P4 exposure before estrus and/or E2 production by the dominant follicle [START_REF] Dobson | Why is it getting more difficult to successfully artificially inseminate dairy cows?[END_REF][START_REF] Walker | Chronic stress, hormone profiles and estrus intensity in dairy cattle[END_REF]. Further, a normal endocrinological cycle is prerequisite for appropriate expression of estrous behavior. The relationships found between P4, E2 and intensity of estrous behavior show that hormones involved in regulation of the estrous cycle also affect the expression of estrous behavior [START_REF] Lyimo | Relationship among estradiol, cortisol and intensity of estrous behavior in dairy cattle[END_REF][START_REF] Roelofs | Various behavioral signs of oestrus and their relationship with time of ovulation in dairy cattle[END_REF]. These and other findings and hypotheses about regulation of the bovine estrous cycle could be translated into mathematical equations or modified parameterization and incorporated in the current model.

Appendix A. List of Hill functions 474

The Hill functions listed below are the full notations of the Hill functions mentioned in Section 3.2 and represent the mechanisms shown in Figure 2.

H - 1 (P4&E2) := m P4&E2 • h -(P4(t); T GnRH ,1 P4 , 2) + h -(E2(t), T GnRH,1 E2 , 2) -h -(P4(t); T GnRH ,1 P4 , 2) • h -(E2(t), T GnRH ,1 E2 , 2) H - 2 (P4) := m GnRH ,2 P4 • h -(P4(t), T GnRH ,2

P4

, 2)

H + 3 (E2) := m GnRH ,2 E2 • h + (E2(t), T GnRH ,2 E2 , 5) H - 4 (Inh τ ) := m FSH Inh • h -(Inh(t -τ Inh ), T FSH Inh , 2) H + 5 (P4) := m FSH P4 • h + (P4(t); T FSH P4 , 2) H - 6 (E2) := m FSH E2 • h -(E2(t); T FSH E2 , 2) H + 7 (GnRH Pit ) := m FSH GnRH • h -(GnRH Pit (t); T FSH GnRH , 1) H + 8 (E2) := m LH E2 • h + (E2(t); T LH E2 , 2) H - 9 (P4) := m LH P4 • h -(P4(t); T LH P4 , 2) H + 10 (GnRH Pit ) := m LH GnRH • h + (GnRH Pit (t); T LH GnRH , 2) H + 11 (FSH ) := m Foll FSH • h + (FSH Blood (t); T Foll FSH (t), 2), T Foll FSH (t) := T Foll FSH • h -(Foll (t); T FSH Foll , 1) H + 12 (P4) := m Foll P4 • h + (P4(t); T Foll P4 , 2) H + 13 (LH ) := m Foll LH • h -(LH Blood (t); T Foll LH , 2) H + 14 (P4 τ ) := m PGF2 α,1 P4 • h + (P4(t -τ P4,1 ), T PGF2 α P4 , 2) H + 15 (P4 τ ) := m PGF2 α,2 P4 • h + (P4(t -τ P4,2 ), T PGF2 α P4 , 10) H + 16 (LH τ ) := m CL LH • h + (LH Blood (t -τ LH ); T CL LH , 2) H + 17 (CL) := m CL CL • h + (CL(t); T CL CL , 2) H + 18 (PGF2 α) := m CL PGF2 α • h + (PGF2 α(t); T CL PGF2 α , 1)

Figure 1 :

 1 Figure 1: Schematic representation of the compartments in the model of the bovine estrous cycle.

Figure 2 :

 2 Figure 2: Complete mechanisms of the bovine model. Boxes represent the 12 key components of the system. Differential equations are derived for these 12 components. Arrows denote functional dependencies. Stimulating and inhibiting effects are indicated by + and -respectively.

Figure 3 :

 3 Figure 3: Scaled positive Hill functions with different steepness.

Figure 4 :

 4 Figure 4: Simulated curves of the closed model together with the data points used for parameter estimation. Panels 4(a), 4(b) and 4(c) show data points based on qualitative behavior of hormones as described in literature ([68]). Panels 4(d), 4(e) and 4(f) show data points obtained from the input curves. Day zero corresponds to the day of LH peak.

Figure 5 :

 5 Figure 5: Output curves of follicular function (Foll) and CL function (CL) over time for one cycle (a) and in consecutive cycles (b).

Figure 6 :

 6 Figure 6: Output curves of serum concentrations of E2 and LH, and portal concentration of GnRH over time for one cycle (a) and in consecutive cycles (b).

Figure 7 :

 7 Figure 7: Output curves of serum concentrations of Inh and FSH over time for one cycle (a) and in consecutive cycles (b).

Figure 8 :

 8 Figure 8: Output curves of serum concentrations of P4 and PGF2α over time for one cycle (a) and in consecutive cycles (b).

  The sensitivity analysis shows that parameters36, 20 and 11 are the most sensitive parameters of the model, which means that a small change in the value of one of these parameters will have a large effect on the model solution. Parameter 36 (delay of P4 until stimulating PGF2α increase) is possibly a sensitive parameter because CL life span is critical for the duration of the cycle. Parameter 20 (FSH clearance rate constant) and Parameter 11 (delay of Inh in FSH synthesis) are possibly sensitive parameters because FSH and Inh serum levels have an important effect on the progress of follicle development.

  ]. The complete

			Hypothalamus
			GnRH synthesis and release
			f(peripheral E2 and P4)
		P4 E2	Portal blood system
			GnRH
			Pituitary
	Peripheral blood	P4 E2 Inh	Synthesis of LH and FSH f(P4,E2,Inh) Release of LH and FSH f(P4,E2,GnRH)			LH FSH Peripheral
	Clearance of					blood
	P4, E2, Inh			LH		Clearance of
				FSH		LH, FSH
			Ovaries
	E2 P4		Follicle growth (E2 synthesis) f(LH,FSH, follicular size)		P4
	Inh		Inhibin synthesis			Uterus
			f(follicular size)	PGF2	α	PGF2	α
			CL growth			synthesis
			α f(LH,CL size,PGF2 )			f(CL size)

Table B

 B Continued on next page... Continued on next page... • (h -(y 10 (t); p 5 , 2) + h -(y 11 (t); p 4 , 2) -h -(y 10 (t); p 5 , 2) • h -(y 11 (t); p 4 , 2)) + p 6 • h -(y 10 (t); p 7 , 2) • y 1 (t) d dt y 2 (t) = p 3 • h -(y 10 (t); p 5 , 2) + h -(y 11 (t), p 4 , 2) -h -(y 10 (t); p 5 , 2) • h -(y 11 (t); p 4 , 2) +p 6 • h -(y 10 (t); p 7 , 2) • y 1 (t) • p 8 • h + (y 11 (t); p 9 , 5)p 10 • y 2 (t) d dt y 3 (t) = p 12 • h -(y 12 (tp 11 ); p 13 , 2)p 14 • h + (y 10 (t); p 15 , 2) +p 16 • h -(y 11 (t); p 17 , 2) + p 18 • h + (y 2 (t); p 19 , 1) • y 3 (t) (t) = p 37 • h + (y 10 (tp 36 ); p 38 , 2)p 40 • h + (y 10 (tp 39 ); p 41 , 10) • y 8 (t) d dt y 9 (t) = p 43 • h + (y 6 (tp 42 ); p 44 , 2) + p 45 • h + (y 9 (t); p 46 , 1) p 47 • h + (y 8 (t); p 48 , 2)

		.1 -continued from previous page Table B.1 -continued from previous page Table B.1 -continued from previous page Table B.1 -continued from previous page
	No. Symbol m GnRH ,2 E2 T GnRH ,2 E2 c GnRH ,2 τ Inh m FSH Inh T FSH Inh m FSH P4 T F SH P4 m FSH E2 T FSH E2 No. Symbol T FSH GnRH c FSH m LH E2 T LH E2 m LH P4 T LH P4 m LH GnRH T LH GnRH b LH c LH m Foll T PGF2 α,2 P4 c Inh Foll FSH No. Symbol T FSH Foll m Foll P4 T Foll P4 m Foll LH T Foll LH τ P4,1 m PGF2 α,1 P4 T PGF2 α,1 P4 τ P4,2 m PGF2 α,2 P4 No. Symbol m CL LH T CL LH m CL CL T CL CL m CL PGF2 α T CL PGF2 α c P4 CL c P4 c E2 Foll c E2 d dt y 1 (t) = p 2 • 1 -Value Quantity 1.5 [GnRH Pit ] [GnRH Hypo ] 1.276 [E2] 1.299 1/[t] 1.5 [t] 1 [FSH]/[t] 0.06 [Inh] 2 1/[t] 0.0966 [P4] 0.3 1/[t] 2.846 [E2] Value Quantity 0.4 [GnRH] 0.8 1/[t] 1.5 [LH]/[t] 0.1 [E2] 4.5 [LH]/[t] 0.0322 [P4] 4 1/[t] 4 [GnRH] 0.05 1/[t] 11 1/[t] 0.0966 [P4] 4.8 1/[t] [Inh]/[Foll] 0.8 [Foll]/[t] Value Quantity 0.3 [Foll] 2.5 1/[t] 0.1127 [P4] 2.8 1/[t] 0.525 [LH] 12 [t] 0.3 [PGF2α]/[t] maximum growth rate of PGF2α Explanation maximum scaling of pituitary sensi-tivity for GnRH threshold of E2 to increase pituitary sensitivity for GnRH GnRH clearance rate constant in the pituitary delay of Inh in FSH synthesis maximum FSH synthesis rate in the pituitary in the absence of Inh threshold of Inh for inhibition of FSH synthesis maximum part of FSH release rate that is stimulated by P4 threshold of P4 to stimulate FSH re-lease maximum part of FSH release rate that is inhibited by E2 threshold of E2 to inhibit FSH re-lease Explanation threshold of GnRH to stimulate FSH release FSH clearance rate constant maximum part of LH synthesis that is stimulated by E2 threshold of E2 to stimulate LH syn-thesis maximum part of LH synthesis that is inhibited by P4 threshold of P4 to inhibit LH syn-thesis maximum part of LH release rate that is stimulated by GnRH threshold of GnRH to stimulate LH release basal LH release rate constant LH clearance rate constant tion stimulated by FSH threshold of P4 to stimulate PGF2α proportionality factor of delayed fol-maximum increase of follicular func-Explanation threshold of follicular function to downscale FSH threshold maximum part of follicular decay stimulated by P4 threshold of P4 to stimulate de-crease of follicular function maximum part of follicular decay stimulated by LH threshold of LH to stimulate de-crease of follicular function delay of P4 until stimulating PGF2α increase 0.1672 [P4] threshold of P4 to stimulate PGF2α increase 17 [t] delay of P4 until stimulating PGF2α decrease 11 [PGF2α]/[t] maximum decay rate of PGF2α Value Quantity Explanation 0.334 [CL]/[t] maximum increase of CL stimulated by LH 1.2 [LH] threshold of LH to stimulate CL in-crease 0.0334 [CL]/[t] maximum increase of CL stimulated by itself 0.0651 [CL] threshold of CL to stimulate self-growth 6.536 1/[t] maximum decrease of CL stimu-lated by PGF2α 2 [PGF2α] threshold of PGF2α to initiate de-crease of CL 3.856 [P4]/[CL] 1/[t] proportionality factor of CL in P4 increase 2.737 1/[t] P4 clearance rate constant 1.9 [E2]/[Foll] 1/[t] proportionality factor of follicular function in E2 increase 0.9 1/[t] E2 clearance rate constant of y 1 (t) p 1 -p 3 d dt y d dt y 8
	m FSH GnRH T Foll FSH	3 0.8	1/[t] [FSH]	maximum part of FSH release rate threshold of FSH to stimulate follic-decrease licular function in Inh increase
	τ LH c Inh	4.5 4	[t] 1/[t]	that is stimulated by GnRH ular function delay of LH in CL Inh clearance rate constant
				Continued on next page...

4 (t) = p 14 • h + (y 10 (t); p 15 , 2) + p 16 • h -(y 11 (t); p 17 , 2) +p 18 • h + (y 2 (t); p 19 , 1) • y 3 (t)p 20 • y 4 (t) d dt y 5 (t) = p 21 • h + (y 11 (t); p 22 , 2) + p 23 • h -(y 10 (t); p 24 , 2) p 27 + p 25 • h + (y 2 (t); p 26 , 2) • y 5 (t) d dt y 6 (t) = (p 27 + p 25 • h + (y 2 (t); p 26 , 2)) • y 5 (t)p 28 • y 6 (t) d dt y 7 (t) = p 29 • h + (y 4 (t); p 30 • h -(y 7 (t); p 31 , 1), 2) p 32 • h + (y 10 (t); p 33 , 2) + p 34 • h + (y 6 (t); p 35 , 2) • y 7 (t) d dt y 10 (t) = p 49 • y 9 (t)p 50 • y 10 (t) d dt y 11 (t) = p 51 • y 7 (t)p 52 • y 11 (t) d dt y 12 (t) = p 53 • y 7 (t)p 54 • y 12 (t)

Acknowledgements

The authors would like to thank Dr B. Beerda, Prof Dr B. Kemp and Prof Dr M. Smits for their helpful comments on the manuscript. These results are obtained through IP/OP: Systems Biology, financially supported by the Dutch Ministry of Agriculture, Nature and Food Quality (BAS no. 4434660700). S. Röblitz and C. Stötzel have been supported by the DFG Research Center Matheon "Mathematics for Key Technologies" in Berlin.

Appendix B. List of parameters

In our model, [•] stands for the unit of the substance, usually a concentration, and can be specified from measurements. Typical units are