Model-based exploration of the impact of glucose metabolism on the estrous cycle dynamics in dairy cows

Background Nutrition plays a crucial role in regulating reproductive hormones and follicular development in cattle. This is visible particularly during the time of negative energy balance at the onset of milk production after calving. Here, elongated periods of anovulation have been observed, resulting from alterations in luteinizing hormone concentrations, likely caused by lower glucose and insulin concentrations in the blood. The mechanisms that result in a reduced fertility are not completely understood, although a close relationship to the glucose-insulin metabolism is widely supported. Results Following this idea, we developed a mathematical model of the hormonal network combining reproductive hormones and hormones that are coupled to the glucose compartments within the body of the cow. The model is built on ordinary differential equations and relies on previously introduced models on the bovine estrous cycle and the glucose-insulin dynamics. Necessary modifications and coupling mechanisms are thoroughly discussed. Depending on the composition and the amount of feed, in particular the glucose content in the dry matter, the model quantifies reproductive hormones and follicular development over time. Simulation results for different nutritional regimes in lactating and non-lactating dairy cows are examined and compared with experimental studies. The simulations describe realistically the effects of nutritional glucose supply on the ovulatory cycle of dairy cattle. Conclusions The mathematical model enables the user to explore the relationship between nutrition and reproduction by running simulations and performing parameter studies. Regarding its applicability, this work is an early attempt towards developing in silico feeding strategies and may eventually help to refine and reduce animal experiments. Reviewers This article was reviewed by John McNamara and Tin Pang (nominated by Martin Lercher).

For ruminants, the energy content in feed cannot be increased without limits due to the fermentative character of the digestive system [7]. High energetic feed with little fiber leads to an imbalance of microbes, rumen acidosis, and may even cause severe illness and death. Nevertheless, targeted feeding strategies are able to extenuate the NEB and to ensure animal health and welfare [5,8,9].
A number of experimental and clinical studies were performed to examine the relationship between the metabolic status and the fertility of cows, both, in qualitative and in quantitative manners, e.g., [10,11]. Reduced nutrition intake was observed to delay the onset of puberty in beef heifers [12][13][14], to change the growth pattern of the dominant follicles (maximal diameter, persistence, number of follicular waves) [15], and to increase the period to conception postpartum [16][17][18]. Studies in the postpartum period of dairy cows showed that the NEB is strongly correlated with low concentrations of glucose, insulin and IGF-1 in the blood [19][20][21]. Changes in the secretion of gonadotropins, caused by low glucose levels, lead to low FSH and LH concentrations [10,22], whereby missing LH peaks cause anovulation [4]. Non-regular estrous cycles are often associated with low average concentrations of insulin in the blood [23]. On the other hand, it was reported that good feed management, e.g., nutritional manipulation that causes increased insulin, reduces the incidence of non-regular estrous cycles [24]. This paper focuses on glucose, as part of the feed and as one of the main energy sources of the body. The aim is to develop a mathematical model that represents metabolic processes as well as reproductive regulation, thus allowing to analyze the impact of glucose originating from the feed on the reproductive hormones and the follicular development.
Previous modeling efforts mainly focused on either the bovine estrous cycle [25][26][27][28][29] or the nutritional strategies [30][31][32], yet there are a few approaches that combine the two topics. The most recent model, named "Jenny", was developed by McNamara and Shields [33]. It connects the reproductive cycle (given by differential equations from [25,26]) with nutrition (implemented by a rather sophisticated model called Molly [31]) via the ATP to ADP reduction reaction. Martin et al. [34] introduced an empirical model that includes nutritional effects on the reproduction. Pring et al. [27] modeled different nutritional scenarios by varying parameters in an estrous cycle model. A more conceptual model was suggested by Scaramuzzi et al. [35], where the coupling between nutrition and reproduction is realized by IGF-1, the glucose-insulin system, and leptin.
None of these models, except [33], captures the dynamics between nutrition, hormonal regulation, and milk yield, mechanisms that are of particular interest in cows. The model of McNamara and Shields [33] contains these elements, as it is based not only on Molly but on the BovCycle model [25,26]. The effort by McNamara and Shields and the effort here are closely related and complementary. However, McNamara and Shields [33] do not include the full reproductive process. The model introduced here aims at understanding the involved interactions and time evolution on a more detailed level. It includes compartments for the nutrient intake, the glucose-insulin system [36], the milk production, and the reproductive hormones [26]. Based on that model, it is analyzed how changes in dietary intake, which usually happen on the time-scale of days, affect the behavior of the estrous cycle on the scale of weeks and months.
The paper is organized as follows. The glucose-insulin model and its coupling to the estrous cycle model are presented in the "Methods" section. The "Results and discussion" section deals with the simulation for nonlactating and lactating cows and compares the outcome with data from literature. Finally, the results are summarized again and limitations of the model are presented in the Conclusion. The model was implemented in MAT-LAB (release 2014b). The code is available in Additional file 1.

Methods
The model that is developed in this section and, later, used for simulations in the "Results and discussion" section is built on two major pillars. The one is the glucose-insulin dynamics in dairy cows, which was modeled in [36] utilizing the Systems Biology Markup Language [37] and CellDesigner [38]. The other is the bovine estrous cycle, modeled by a system of differential equations (BovCycle) that quantifies reproductive hormones and other relevant compartments, representing follicles and corpora lutea [25,26].
The model here consists entirely of ordinary differential equations (ODEs), which are solved for problemspecific initial conditions and parameter values. One half of the model ( Fig. 1 and r.h.s. of Fig. 2) implements the mechanisms explained in [36], which allows for simulating the time-evolution of glucose and insulin for different dietary inputs in lactating as well as nonlactating cows. The other half (l.h.s. of Fig. 2) implements the biological feedback mechanisms between hypothalamus, pituitary gland and ovaries, which produces periodic estrous cycles of constant duration, similar to [26]. However, modifications needed to be implemented as the mechanisms suggested in [26] are not tailored to cows during pregnancy, calving and lactation. In these stages the interaction between hormones is somewhat different. To simulate the onset of lactation, oxytocin is included in the model; this hormone peaks during delivery [39], and it is required for milk ejection [40][41][42].

Metabolic model
The metabolic model to be developed in this section is based on an improved version of the glucoseinsulin model in [36]. It involves six components (Glu blood , Glu liver , Glu store , Fat, Ins, Gluca; see Table 2) and, as formulated here in terms of ODEs, their explicit interaction over time. Initial conditions are chosen based on the following calculation. For a cow of weight 600 kg and body condition score 3.5, the total body fat can be estimated by 25% of the total body weight [7,43]. That is, Fig. 2 Schematic representation of the coupled metabolic-reproductive model. The coupled model links the metabolic model (right hand side) to the bovine estrous cycle model [26] (left hand side). Red arrows depict the sites where both models are coupled. Insulin acts on the site of anterior pituitary influencing LH and FSH release to the blood circulation. Insulin stimulates IGF-1 levels in the blood. Progesterone inhibits IGF-1 secretion which in turn decreases the responsiveness of follicular cells to LH Gluca 50-120 ng/L [81,100] 150 kg is taken as initial value for Fat. Typical physiological ranges for Glu blood , Ins and Gluca are listed in Table 1.
As long as the initial values are within these ranges, they do not affect the performance of the model. The model only involves the most basic mechanisms that regulate the flow of glucose through the body. It starts with the feed, continues with the digestive system and the blood, and ends up with glucose usage. Glucose and glucogenic substances are ingested with the dry matter intake (DMI). In the liver, the glucogenic substances are converted to glucose via gluconeogenesis. Glucose is used for maintanance and milk production, it is stored as glycogen or, after conversion, as fat. The compartments of the model and their interactions are illustrated in Fig. 1. Flows and regulatory mechanisms are summarized in Table 3 and explained in detail in the following subsections.

Feed intake
The first step involves the quantification of the amount of substances in the DMI that are either available for gluconeogenesis in the liver or directly absorbable as glucose into the blood. There exist empirical formulas that estimate the DMI needed to meet the energy requirements; these formulas are based on the cow's body weight (BW) and the net energy (NE) of the diet; see, e.g., [7]. Throughout the paper, a standard cow with body weight 600 kg is considered, and the value for DMI of 11700 gram per day (g/d) is adopted from [36]. This value also results from a formula in [7], assuming a diet's net energy of 1.32 Mcal/kg. 1 Ruminants digestion involves fermentation, which makes consumption of a high-fiber diet possible and necessary [44,45]. In the default setting, the fraction of glucose and glucogenic substances in the DMI, glu pool , is assumed to be 8% of the total DMI, where c 0 is a mass-fraction parameter (with default value c 0 = 0.08) that allows for varying the total amount of  The initial values are used to solve the differential equations glucose and glucogenic substances that can be extracted from DMI. This fraction combines glucose precursor substances such as short chain fatty acids, which are converted to glucose in the liver by gluconeogenesis, as well as glucose that can directly be absorbed from the digestive tract into the blood [45][46][47]. In the cow, only very little glucose is available for direct absorption from the digestive tract [48]. From the total amount of glucose and glucogenic substances in the DMI (glu pool ), the portion of glucose was estimated to be less than 10% [49][50][51], whereas up to 90% of glu pool are glucogenic substances. The flow of absorbable glucose that goes directly to the systemic circulation is incorporated into the model via the rate The flow of glucose precursor substances that are converted to glucose by gluconeogenesis in the liver is incorporated into the model via the rate The default parameter value is c 1 = 0.08 (cf. Table 4). It is assumed here that there is no loss from the glucose pool (the flows sum up to 1·glu pool ), i.e., the processes take place with 100% efficiency. If some loss was included here, the simulation results presented further below would be t he same but correspond to higher values of c 0 (the amount of glucose and glucose precurser substances in the feed).

Insulin and glucagon
The blood glucose concentration is maintained at normal levels primarily through the action of two hormones, namely insulin and glucagon. Any elevation in the blood glucose concentration leads to the production of insulin in the pancreatic beta cells. Insulin promotes glucose uptake in target cells, e.g., those in the liver, muscles and fat tissue, and it promotes the conversion of glucose to glycogen (glycogenesis) in the liver [52]. When the glucose blood concentration is low, the pancreatic alpha-cells produce glucagon. Glucagon increases the plasma glucose concentration by stimulating the generation of glucose from non-carbohydrate substrates (gluconeogenesis) and the breakdown of glycogen to glucose (glycogenolysis) in the liver [52]. In the model here, the dynamics of the blood insulin and glucagon concentrations are determined by their secretion rates (ins sec , gluca sec ) and their degradation rates (ins deg , gluca deg ), with linear degradation rates It is assumed that the insulin secretion rate increases when the glucose concentration in the blood is above a certain threshold value (T 1 = 0.5 g/L = 2,77 mmol/L), whereas the glucagon secretion rate decreases whenever the glucose concentration in the blood is above that threshold value (T 2 = 0.5 g/L = 2,77 mmol/L), ins sec = c 2 · H + (Glu blood , T 1 , 10), The symbols H + and H − denote a positive and a negative Hill function, which are used to model threshold-dependent stimulatory or inhibitory effects. Here, S ≥ 0 denotes the substance, T ≥ 0 the threshold, and n ≥ 1 the Hill coefficient. A Hill function is a sigmoidal function between zero and one that switches at the threshold S = T from one level to the other with a slope specified by n and T. Threshold kinetics were selected to account for rapid adaptivity, which is an important mechanism to keep the plasma glucose concentration within the physiological range. There are no reference values for the individual rate constants c 2,3,4,5 , but their values were chosen such that a constant glucose blood concentration of Glu blood = T 1 = T 2 = 0.5 g/L (resulting in a Hill function value of 0.5) would give rise to constant insulin and glucagon concentrations that are within the physiological range, namely 0.5 · c 2 /c 3 = 20 mU/L and 0.5 · c 4 /c 5 = 100 ng/L, respectively, compare Table 1.

Glucose production and storage in the liver
When the glucose blood level rises above a certain threshold (T 3 = 0.45 g/L = 2,77 mmol/L), insulin promotes the absorption of glucose from the blood into liver cells (rate glu bl−lv ), Insulin also stimulates the conversion of glucose available in the liver (Glu liver ) to glycogen (glycogenesis rate glu lv−st ). It is assumed here that this rate decreases when the cow produces more than a certain amount of milk (threshold T 4 = 10 L) per day in order to make more glucose available for milk production. In addition, the rate glu lv−st is switched off when the glycogen store, Glu store , reaches the maximal carrying capacity K = 1000g 2 . The equation that describes this process is given by In addition, insulin promotes the absorption of glucose into fat cells and its conversion into triglycerides via lipogenesis. It is assumed here that this process is enhanced once the glycogen storage Glu store is full (threshold T 6 = 1000g). Again, similar to the glycogenesis rate glu lv−st , the rate is assumed to decrease when the cow produces more than a certain amount of milk (threshold T 5 = 10 L) per day, When nutritional supply with glucose is insufficient, the glucagon concentration increases and stimulates the breakdown of glycogen to glucose in the liver (glycogenolysis) to maintain blood glucose homeostasis [55]. This process is assumed to slow down when the glycogen store is below a certain threshold (T 7 = 10g), In this case, i.e., when the glycogen store falls below a threshold (T 8 = 10g), glucagon additionally stimulates the breakdown of lipids into glycerol and free fatty acids 2 Berg et al. [53] estimated that 2% of the weight of the muscle tissue is formed by glycogen, and 10% of the liver weight. Is was also estimated that for a cow with 600 kg body-weight the mass of muscle, liver and kidney is around 280 kg, whereof 9 kg is liver weight [54]. According to these numbers, the liver stores about 900 g glycogen.
(lipolysis) in adipose tissue and the conversion of glycerol into glucose via gluconeogenesis in the liver. This rate is assumed to slowly decrease whenever the total body fat becomes smaller than a certain threshold (T 9 = 150 kg), Finally, glucagon stimulates the release of glucose synthesized in the liver (Glu liver ) into the blood, In the equations above, threshold kinetics were selected for Glu store to differentiate between full end empty store, without modifying the rates in dependence on the actual amount of glycogen in the store.
There are no reference values for the rate constants c 6 to c 11 . They were fitted manually such that the simulation results qualitatively agree with the results reported in literature.

Glucose utilization
All organs and tissues of dairy cows use glucose, except adipose tissue which cannot directly convert glucose to fatty acids [45]. Glucose provides energy for maintenance and production. In the milk producing dairy cow, glucose utilization is dominated by the requirements of the mammary gland for milk synthesis [56]. These requirements increase rapidly right after parturition [57]. Glucose utilization is modeled here in terms of two different sink terms, one from Glu liver , and one from Glu blood , The sink term from Glu blood accounts for maintenance (1st term) and milk production (2nd term). Maintenance refers to glucose utilization by non-mammary tissues including brain and skeletal muscle, but excluding liver. For example, glucose that is stored in skeletal muscle as glycogen cannot be released back into the bloodstream due to the absence of glucose-6-phosphatase. It is assumed here that the glucose consumption for maintenance decreases when the glucose blood level drops below a certain threshold (T 10 = 0.5 g/L = 2,77 mmol/L). The second term accounts for glucose utilized for milk production, including substance and energy. The variable Milk quantifies the daily milk yield in kg/day, whereas the parameter c 13 = 72 g/kg [58] quantifies the amount of glucose (in gram) per kg of milk. Hence, the mammary glucose requirement in a cow with a daily milk yield of 40 kg would be about 3 kg per day. There is no reference value for the non-mammary glucose requirement, but according to the literature [56] this value should be significantly lower (here, c 12 = 1 kg/day was chosen).

The system of differential equations
The final set of ordinary differential equations modeling the dynamics of the glucose exchange reads where V = 22.8 L is the extracellular volume of blood [36]. The ordinary differential equations were solved using the software MATLAB. The parameters and the initial values are listed in Tables 2, 4, and 5, respectively.

A metabolic-reproductive model
Several studies have shown that the metabolic status has a large influence on growing cattle and on reproductive performance in dairy cows. During negative energy balance, which can be caused, e.g., by dietary restrictions or high milk yield, a remarkable change occurs in the levels of the metabolic components IGF-1, insulin, and glucose in the systemic circulation, which in turn influences the levels of reproductive hormones and follicular development [19][20][21]. The aim is to reproduce these observations by coupling the metabolic model and the reproductive model BovCycle introduced in [25,26]. The initial values for the species in the BovCycle model are listed in Table 6. The flowchart for the coupled model is presented in Fig. 2. Detailed explanations of the coupling mechanisms are given in the three sections below.

IGF-1 and insulin
Kawashima et al. [59] reported that IGF-1 is positively correlated with the level of feed intake. The authors argue that the plasma IGF-1 concentration increases transiently during the follicular phase and decreases during the luteal phase of the estrous cycle, i.e., IGF-1 levels decrease when progesterone (P4) increases. On the other hand, IGF-1 is lowest during early lactation when there is no P4 in circulation, and highest in late lactation [60]. In particular, a decrease in blood insulin and glucose concentrations in postpartum cattle is associated with the decrease in IGF-I [21]. In addition, acute dietary restrictions reduce both insulin and IGF-1 concentrations in the blood [4,61]. Even if these are only empirical observations and evidence for mechanistic relationships is missing, these observations are incorporated into the equation for IGF-1 as follows, where c 17 accounts for the basal IGF-1 synthesis rate. The rate constants c 17,18,19 were determined such that the simulated IGF-1 concentrations match with the experimental data from 13 Holstein cows [59], see Fig. 3b. Moreover, in order to fit the simulated progesterone concentrations to the data (Fig. 3a), the basal P4 production rate had to be increased from c P4 = 0 in the original model [26] to c P4 = 0.1. This is consistent with reports about baseline progesterone levels [62]. A change in plasma IGF-1 has an impact on follicular cell development and responsiveness to hormonal signals. In particular, experimental studies demonstrated that reduced IGF-1 reduces ovarian responsiveness to LH stimulation [21,63]. To include this mechanism in the model, the term in [26] that models the follicular cell responsiveness to LH, was improved as follows. The LH blood concentration that is required for an ovarian response (threshold T 13 ) is made dependent on IGF-1, Such a dependency was chosen because it allows for LH concentrations to increase in response to IGF-1 being below a certain threshold, T 14 . This mechanism is essential to ensure appropriate ovarian responses to IGF-1.
Insulin serves as a metabolic signal influencing the release of LH and FSH from the anterior pituitary into the blood [21,64]. This mechanisms is included in the model by a stimulatory effect of insulin on the synthesis rates of LH and FSH. The equations for LH and FSH in [26] are changed to where LH syn , FSH syn , LH rel , and FSH rel are the synthesis and release rates of LH and FSH, respectively, as described in [26].
Hence, if insulin levels drop below a certain threshold (T 15 = T 16 = 21 mU/L), the synthesis of LH and FSH halts.

Lactation
Pregnancy and calving are characterized by a complex interplay of hormones. One of these hormones is oxytocin. The release of this hormone and milk yield are positively correlated [41]. Overall as well as peak concentrations of oxytocin decrease over one ongoing lactation [65]; earlier studies reported similar dynamics [66][67][68][69]. According to measurements in those studies, peak concentrations of oxytocin during early lactation are more than twice the magnitude of those during late lactation.
The BovCycle model [26] does not capture changes in oxytocin concentrations during pregnancy and calving. Holstein dairy cows (red dots) were collected and kindly provided by Kawashima et al. [59] To this end, the model was extended by introducing an additional term Oxy lac into the equation of oxytocin, with This is the simplest form of a non-negative decreasing function, namely a Gaussian function, see Fig. 4. The parameter value c 25 = 0.0007 determines the width of the curve and was adopted to the approximate length of the early lactation period, whereas the parameter value c 24 = 1.5 was fitted so that Oxy(t) during early lactation is about twice as high as Oxy(t) during late lactation.

Reparametrization of the BovCycle model
The changes in the equations of the original BovCycle model [26] required changes of some of the original parameter values in order to be able to recover regular estrous cycles. In addition, the original BovCycle model [26] was challenged with the scenario of adding exogenous oxytocin early in the cycle. In a study by Donaldson et. al [70], it was shown that daily oxytocin injections to eight non-lactating cows starting on day two of the cycle reduced the estrous cycle length to nine days. The slow increase in plasma P4 concentration during the first five days of the cycle was not altered significantly, but plasma P4 concentrations decreased again to low values after day five. These results confirmed earlier studies [71,72]. However, the original BovCycle model [26] did not reproduce these results. Hence, changes were made on parameters that describe the interaction of oxytocin and enzymes with prostaglandin F2 α and the interovarian factor such that the recalibrated model correctly reflects the effects of oxytocin administration on the length of the estrus cycle and plasma P4 concentrations. Parameters that required changes are listed in Table 7.

Sensitivity analysis
Sensitivity analysis aims at determining the model input parameters which mostly contribute to a quantity of interest depending on the model output. Let us denote the model input parameter vector as p = (p 1 , . . . , p d The model here is an ordinary differential equation model of the form and a quantity of interest, y, can be any observable depending on the model output x, y = y(x(t, p)).
This quantity can be for instance the value of a specific output variable x j at a specific time point t, or the variance of x j over a specific time interval. These are examples for scalar outputs. For the sake of simplicity, the study here is restricted to a scalar output y. The sensitivity of y with respect to input parameter p i is given by To account for differences in physical units among variables and parameter, often relative sensitivities are used, If the exact derivative is difficult to compute, the sensitivity can be approximated by a finite difference scheme, where is the size of the perturbation and e i is a vector of the canonical base. Often, is a relative perturbation, i.e., = · p i for some small number (e.g. = 0.1) corresponds to a perturbation by 10%. In this case, the relative sensitivity is approximated bŷ  This is a local sensitivity in the sense that it describes the influence of a specific local perturbation of parameter p i on the model output. Sampling or sampling pairs of input and output variables would allow for a global sensitivity analysis, but this is computationally much more demanding and the results are often difficult to interpret. For details on global sensitivity analysis, the reader is referred to [73].
For the metabolic-reproductive model presented here, sensitivity analysis is performed to determine the model parameters that are most important for the onset of luteal activity after calving. Hence, the observable y is chosen as the earliest time point at which the (relative) P4-level is larger than a threshold T P4 = 1, y(x(t, p) The results of this analysis are presented in the following section.

Results and discussion
The aim of this study was to analyze the impact of supplied glucose, represented by the parameter c 0 , on the estrous cycle dynamics in both lactating and non-lactating cows. For this purpose, the model was simulated for different feeding scenarios, including short and long time dietary restrictions. For a cow of 600 kg BW, DMI at maintenance is set to its default value of 11.7 kg/d [36]. This is the reference value corresponding to 100% DMI throughout the following, and variations to this value are stated accordingly.

Non-lactating cows
To model these cows, the value of Milk in Eq. (14) is set to zero. The numerical experiments for acute and chronic dietary restrictions are designed according to three experimental feeding studies from Mackey et al. [74], Murphy et al. [15] and Richards et al. [75]. Since these are studies in beef heifers and anestrus beef cows, respectively, the results are expected to agree only qualitatively, not necessarily quantitatively.

Varying the glucose content in the DMI
The effect of varying glucose content in the DMI on the glucose-insulin dynamics is analyzed by changing the value of the parameter c 0 (glucose content in the DMI) between 4%, 8% and 16%. Simulation results are presented in Figs. 5 and 6.
At maintenance intake, i.e. c 0 = 0.08, the model calculates the non-mammary usage to be slightly less then 400 gram per day (Fig. 6d). This number is in qualitative agreement with Danfaer et al [76], who estimated the amount of glucose required for maintenance in a non-lactating cow with a slightly lower body weight of 500 kg to be 290-380 gram per day. The amount of glucose absorbed from the digestive tract directly into the blood is calculated to be 75 g/d (Fig. 6a). The calculated amount of glucose released from the liver into the blood is about 800 g/d (Fig. 6c). This means that the total amount of glucose available in the blood is around 875 g/d, whereas the glucose uptake into liver cells (Fig. 6f ) and the non-mammary usage (Fig. 6d) sum up to the same amount. This balance between input and consumption of glucose leads to stable glucose and  . 6 Simulated metabolic rates in non-lactating cows at maintenance. Glucose content in the DMI was fixed at 8%. The figure illustrates glucose input, storage, and usage in terms of the amount of glucose absorbed via the digestive tract (a), glucose generated from glucogenic substances in the feed (b), glucose released from the liver into the blood (c), glucose absorbed into liver cells (f), and glucose used for body maintenance (d) and for metabolic processes in the liver (e). At maintenance intake, the cow is able to cover the daily glucose requirement, which results in stable levels of glucose in the different compartments insulin levels in the blood (Fig. 5a, d). In addition, this leads to stable glycogen and fat levels in the respective storage components (Fig. 5b, e).
With increasing glucose content in the DMI (c 0 = 0.16), more glucogenic substances are available and lead to an increased gluconeogenesis [45]. This increases glucose and insulin concentrations in the blood, but they are still within their physiological range (Fig. 5a, d). Excess glucose in the system is stored as glycogen or fat reserves (Fig. 5b,  e). When the glucose content in the DMI is decreased to 4%, blood glucose and insulin levels decrease towards their lower physiological bounds within two days (Fig. 5a,  d), compare Table 1. As a result, the stored glycogen and the fat reserves (Fig. 5b, e) are reduced as well.

Acute nutritional restriction
To simulate the effect of acute nutritional restriction on the estrous cycle, a numerical experiment was designed according to the study of Mackey et al. [74], who reported about the effect of nutritional deprivation for a period of 13-15 days. Heifers with 406 ±5 kg body weight were allocated to a diet with a DMI of 1.2% of body weight for maintenance and then reduced to a diet with a DMI of 0.4% of body weight. In the model here, this reduction to 1/3 of the default diet corresponds to a reduction in the DMI from 11.7 kg/d to 3.84 kg/d. This acute nutritional restriction is applied immediately after ovulation. The simulation results show increased levels of P4 (Fig. 7d), indicating a failure of luteolysis. Anovulation can be attributed to the absence of LH pulses (Fig. 7a) and lower FSH levels (Fig. 7b), as a result of decreased insulin levels (Fig. 7f ). In addition, IGF-1 is decreased during the dietary restriction (Fig. 7e), which negatively influences the responsiveness of follicular cells to LH [20].

Chronic nutritional restriction
To simulate the effect of chronic nutritional restriction on the estrous cycle, numerical experiments were designed according to the studies of Murphy et al. [15] and Richards et al. [75]. Murphy et al. [15] examined the effect of chronic dietary restriction on the estrous cycle over 10 weeks. In this study, heifers with 375 ±5 kg body weight were allocated to a maintenance diet with an amount of DMI corresponding to 1.2% of the body weight and a reduced diet with 0.7% of the body weight. In the model here, this reduction to 58% of the maintenance diet corresponds to a reduction in the DMI from 11.7 kg/d to 6.79 kg/d. In the experiment by Richards et al. [75], multiparous non-lactating Hereford cows underwent a chronic nutritional restriction for 30 weeks. They were fed to lose 1% of their bodyweight weekly. After the restriction period, the diet was increased to 160% of the maintenance diet.
The simulation was adapted to these two scenarios as follows. The nutritional restriction starts after ovulation. From then on, the model was simulated with 58% of the maintenance DMI within a time interval of 30 weeks. Simulation results (Fig. 8) show that the cow exhibits normal estrous cycles over a period of 15 weeks. During the chronic restriction period, the glycogen store (Fig. 8g) and the insulin in blood (Fig. 8f) decrease. LH (Fig. 8a), FSH (Fig. 8b) and IGF-1 (Fig. 8e) pulses decrease in frequency and amplitude, resulting in cessation of cyclicity after 15 weeks of feed restriction. The fat compartment loses around 10%. After 15 weeks, P4 decreases to a low level for the remaining 15 weeks, indicating the onset of anestrus. FSH and E2 exhibit changes in their wave patterns, that is, the number of waves per cycle increases. A similar tendency was observed in [15].
Murphy et al. [15] examined ultrasound data and serum P4 between week 6 and 9. They found no alteration in CL growth, whereas P4 in restricted cows was numerically higher than in cows on maintenance diet. No anestrus was observed in the first 10 weeks of the restriction period, which is in agreement with the simulation results.
During the first weeks of restriction in the experiment by Richards et al. [75], P4 concentration increased as well. After losing 24.0 ± 0.9% of their initial body weight, cows had decreased luteal activity measured via P4, and cessation of the estrous cycle was observed in 54% of the cows after 26 weeks. The authors reported that estrous cycles were re-initiated by week 40 in 64% of the restricted cows, feeding 160% of maintenance diet. The model predicts reinitiation of cyclicity by week 32, feeding 160% of DMI at maintenance.

Lactating cows
To investigate the effect of lactational metabolism and NEB on fertility hormones, different scenarios were simulated with the metabolic-reproductive model. As model input, interpolated time series data of DMI and milk yield from a study by Friggens et al. [77] were used, see Fig. 9. Each kilogram of milk produced requires around 72 gram glucose (parameter c 13 in Eq. (14) [58]. Hence, the production of 41 kg milk per day requires about 3 kg of glucose per day. This was confirmed by Reynolds et al. [78], who predicted the glucose usage for milk to be between 2500 g/d and 3000 g/d. Milk production and the provided DMI in this study were 41 kg/d and 21 kg/d, respectively, averaged over 5 Holstein cows with an average body weight of 647 kg.
Energy balance is usually calculated as energy input minus output, requiring measurements of feed intake and energy output sources (milk, maintenance, activity, growth, and pregnancy) [79]. Alternatively, the energy balance can be calculated based on changes in the body reserves, using body weight and body condition score Fig. 7 Effect of acute dietary restriction on the bovine estrous cycle in non-lactating dairy cows. On day 43, DMI is reduced from 100% (11.7 kg/d) to 33% (3.84 kg/d) for 15 days (the time period bounded by the two red lines). During the restriction period, one can observe a decrease of glucose in the store (g), insulin in the blood (f) and IGF-1 (e), an absence of LH pulses (a), and a decrease of amplitude in the FSH waves (b), leading to anovulation and failure in luteolysis with increasing P4 (d). The cycle re-starts soon after the end of the restriction period [79,80]. Since the model presented here does not explicitly calculate the body weight, the change in body fat is considered as an indicator of the energy balance, This approach was also used in [81].

Varying the glucose content in the DMI
To explore the metabolic processes during lactation, simulations were performed for different values of glucose content in the DMI (parameter c 0 ). The results are compared qualitatively with the studies by Elliot [82] and Reynolds et al. [78]. The changes in the glucose-insulin dynamics, body fat reserves, and metabolic rates are illustrated in Figs. 10, 11, and 12, respectively. The simulation results clearly show a non-linear relationship between glucose content in the DMI and the values of glucose in blood and storage as well as insulin in blood at peak milk. Decreasing the glucose content in the DMI, starting from c 0 = 0.3, first leads to a slow decrease in glucose and insulin levels, followed by a rapid decrease if c 0 approaches the value 0.2.
For a high amount of glucose in the DMI (30%, c 0 = 0.3), glucose and insulin levels in the blood are maintained within their physiological range (Fig. 10a, d). After the peak milk phase, the cow is even able to store glucose and fat (Fig. 10b, e). Consequently, the overall energy balance is positive throughout the lactation period (Fig. 11d). The model calculates the amount of glucose available in the circulation by direct absorption from the digestive tract (rate glu feed−bl ) to be between 500 and 600 g/d (Fig. 12a). This is in agreement with Elliot [82], who estimated that for a cow with 600 kg BW and a milk yield of 40 kg/d, the amount of glucose absorbed from the digestive tract is around 600 g glucose per day.
For medium amounts of glucose in the DMI (22.5% or 25%), glucose and insulin levels are still kept within their physiological range (Fig. 10a, d), but the period of negative energy balance is prolonged (Fig. 11b, c).
If the amount of glucose in the DMI is decreased even further (20%, c 0 = 0.2), one can observe an extended phase of negative energy balance with glucose and insulin dropping towards their lower physiological limits around peak milk (Fig. 10a, d). High demand and low input trigger the mobilization of body reserves, represented in the model by glycogen and fat in the store (Fig. 10b, e).
When c 0 is varied between 0.2 and 0.3, the calculated amount of glucose released from the liver (glu prod ) Fig. 8 Effect of chronic dietary restriction on the bovine estrous cycles in non-lactating cows. DMI is reduced to 58% for 30 weeks (period between the red lines) and increased to 160% afterwards. During the restriction period, the glycogen store (g) and insulin in blood (f) decrease. LH (a), FSH (b) and IGF-1 (e) pulses decrease in frequency and amplitude, resulting in cessation of cyclicity after 15 weeks of feed restriction within the first 83 days post partum is 2500-4400 g/d (Fig. 12c). These numbers are in qualitative agreement with Reynolds et al. [78], who estimated the daily glucose production in the liver within the first 83 days post partum to be between 2700 and 3600 g/d. On can also observe that the mammary glucose usage gets prioritized compared to the non-mammary usage (Fig. 12f, d), and that this effect becomes more pronounced for low glucose diets.

The effect of changing glucose in the DMI on the estrous cycle
The glucose content in the DMI (parameter c 0 ) has an effect on the estrous cycle. In the previous subsection, it When c 0 = 0.2 (corresponding to 20% glucose content), glucose levels during peak milk drop towards the physiological limit (0.39 g/L) (a). In general, low amounts of glucose lead to a rapid depletion of the store (b), accompanied by a decrease in body fat (e), indicating a negative energy balance due to high milk production was shown that decreasing c 0 from 0.3 to 0.2 prolongs the phase of negative energy balance. A decrease in blood glucose and insulin concentrations is associated with a decrease in IGF-I [83][84][85]. As a consequence, elongated postpartum anestrus periods occur [86][87][88][89]. Similarly, Walsh et al. [5] resumed that NEB leads to low insulin concentrations in blood, which in turn prevents an increase in IGF-1 secretion, resulting in delayed resumption of ovarian cyclicity [90].
The simulation results (Fig. 13) agree with those observations. Increasing the relative amount of glucose in the DMI from c 0 = 0.2 to 0.3 increases the IGF-1 concentration. This stimulates the responsiveness of follicles to LH, thereby shortening the postpartum anestrus interval from about 150 to 40 days (Fig. 14). Accordingly, the oxytocin level becomes cyclic again at the end of the anestrus interval, after having significantly decreased over the postpartum period (Fig. 15).
The length of the postpartum anestrus in the simulations agrees with the literature. In studies based on postpartum progesterone profiles, it was demonstrated that 90 to 95% of post partum dairy cows have resumed ovarian cycles by day 50 after calving [91][92][93]. Hence, a postpartum dairy cow is considered 'normal' if it has resumed ovarian cyclicity by day 50 post partum and continues cycling at regular intervals of approximately 21 days [94]. The simulations also show that estradiol levels at the beginning of the lactation period are within their normal range. This was confirmed by several studies. The authors in [75] found that restricted nutrition leads to anovulation but does not alter estradiol blood concentrations. Although ovulation and luteal development do not occur in anestrus cows, follicular growth is not totally impaired by restricted nutrient intake. In a review, Diskin et al. [10] suggested that NEB in early lactation does not affect the follicle population but does affect the ovulatory fate of the first dominant follicle. The authors summarized that low IGF-I and insulin cumulatively reduce follicular Fig. 13 Simulated levels of P4, IGF-1, LH and estradiol during lactation for different values of glucose content in the DMI. Hormonal cycles were simulated over the lactation period for different fractions of glucose in DMI (parameter c 0 ). A lower glucose content results in negative energy balance (Fig. 11), thereby prolonging the anestrus period. A higher glucose content results in an improved energy balance, which leads to increased insulin and IGF-1 levels and an earlier re-start of the estrous cycle . Simulated data (red dots), which represents the estimated incidence of first ovulation, is determined by the time of first LH peak followed by an increase in progesterone production above baseline. The blue line represents the fitted curve f (x) = a · exp(−b · x) + c to the data with a = 45581, b = 0.30317, c = 35.644. A lower glucose content results in a late ovulation, whereas a higher glucose content results in an early ovulation responsiveness to LH and ultimately suppress follicular oestradiol production.
There is evidence that a good management of the diet can reduce the incidence of abnormal estrous cycles [23,24,27]. Improving postpartum nutrition increases the blood concentration of insulin and IGF-I, which ultimately enhance LH pulsatility [19,85]. Higher IGF-1 levels during the first two weeks postpartum lead to an earlier re-start of the estrous cycle [5]. It was demonstrated in a study that providing a diet high in starch promotes an increased insulin release with a subsequent rise from 55% to 90% in the number of cows that ovulated within 50 days postpartum [24], a time interval that is considered to be an indicator for good reproductive performance [91]. In sum, resumption of cyclicity during lactation is crucial for good fertility in dairy farming. It can be influenced by feed intake, but also depends on many other factors such as uterine health, metabolic status, milk yield and overall condition.

The effect of changing model parameters on the estrous cycle
A local sensitivity analysis as described in Eq. 27, was performed to assess the influence of all model parameters on the time of first ovulation after calving, characterized by the onset of luteal activity (increased P4 levels). Throughout the calculations, glucose content in the DMI was fixed at c 0 = 0.25, which resulted in an onset of luteal activity at day 50 post partum. The parameters' impact on the timepoint of ovulation is illustrated in Fig. 16. Figure 16a shows the change in the timepoint of first ovulation after perturbation of single parameters by +10%, whereas Fig. 16b shows the change in the timepoint of first ovulation after   27). a shows the change in the timepoint of first ovulation after perturbation of single parameters by +10%, whereas b shows the change in the timepoint of first ovulation after perturbation by -10%. Note that in the two subplots a and b only the numerator ofŜ i y is plotted since the denominator is independent of the parameter index i. The two most influential parameters are T 1 (parameter 91) and T Foll P4 (parameter 33). A change of the parameter T 1 by +10% results in a later occurrence of ovulation (Fig. 16c). On the other hand, a decrease in the value of T Foll P4 by -10% stimulats the decay of follicular function, which causes a prolongation of the anovulatory period to day ≈ 120 perturbation by -10%. Note that in the two subplots (A) and (B) only the numerator ofŜ i y is plotted, since the denominator is independent of the parameter index i. The two most influential parameters are T 1 (parameter number 91) and T Foll P4 (parameter 33, described in [26]). The first one describes the threshold of glucose in the blood to stimulate insulin secretion, while the second one is the threshold of P4 to stimulate decrease of follicular function. A change of the parameters 91 and 33 by +10% and -10%, respectively, results in a later occurrence of ovulation (Fig. 16c). Indeed, an increase in the value of T 1 by 10% limits the secretion rate of insulin. As insulin influences the release of LH, LH pulses are suppressed, which delays the ovulation to day ≈ 90. On the other hand, a decrease in the value of T Foll P4 by -10% stimulates the decay of follicular function, which causes a prolongation of the anovulatory period to day ≈ 120.

Conclusion
In the previous sections, the relationship between fertility and metabolism was explored based on two validated models [26,36]. These models were slightly modified and coupled to simulate the interplay of follicular development and its hormonal regulation with the glucose-insulin system. Information about the mechanistic interactions between fertility and metabolism, if taken straight from the literature, is sometimes contradictory and/or redundant. Therefore, only a small number of mechanisms was included, sufficient to realize the coupling of the two models.
With the coupled model, acute and chronic dietary restriction scenarios were simulated, intending to reproduce clinical study findings for non-lactating cows [74,75,95]. Furthermore, numerical experiments were run by varying the amount of DMI and the glucose content in the DMI for both lactating and non-lactating cows, and the effect of dietary restrictions on the estrous cycle was analyzed in lactating cows. The simulation results agree with the findings from the clinical studies, at least on a qualitative level.
The graphs presented here show the same qualitative behavior as the graphs supplied by McNamara and Shields [33]. This is not surprising since the model in [33] uses parts of the BovCycle model [25,26] and combines it with the more complicated Molly model. Despite using independent data sets (compared to the study by McNamara and Shields), the parameter values and trends came out very similarly.
The model here has also some limitations. Increasing (decreasing) the glucose content in the DMI, given by the parameter c 0 , results in the same simulation output as increasing (decreasing) total DMI, because only the product c 0 · DMI is contained in the model equations but not the individual factors. In reality, this is certainly not true. A way out would be to relate DMI directly or indirectly (e.g., via metabolic activity as in [36]) to one of the other variables. However, this would have complicated the model structure which, from the authors' point of view, is not necessary for the modeling purpose in this paper. Furthermore, the model presented here only describes processes in a single representative cow. In its current form, the model is not able to display inter-or intraindividual variability. However, since the implemented mechanisms are universal, variability could easily be included by adapting parameter values to individual measurements, once such measurements are available.
One could also criticize the model for its restriction to glucose as the only feed component. Hence, the protein content should be included in addition to glucose and fat to complete energetic composition of DMI. This would provide one with a more realistic nutrient supply, change of body composition and body weight as well as milk production and composition.
In addition, experimental research is gaining more and more insights into the effect of nutrition on follicular development. With an improved follicle model, similar to the one introduced in [96], further simulations can be conducted to explore the effect of nutrition on multiple follicles in more detail.
So far, it is fair to say that the model presented here is only a starting point. It will certainly be modified and improved in future. However, by conducting numerical simulations relying on it, it was confirmed that an appropriate nutritional intake is fundamental in mitigating the effects and the extension of NEB in order to reduce the incidence of metabolic disorders in high producing cows and to avoid subsequent fertility problems [1,5,8,97]. To understand the interaction between nutrition, metabolism and reproduction, a unified approach was followed, similar to [33,98], where these fields of interest are integrated in one mathematical framework. The model here, formulated in terms of differential equations, enables the user to explore the relationship between nutrition and reproduction by performing related parameter studies. The local sensitivity analysis with respect to the onset of luteal activity after calving is just one example for such an analysis, which can easily be extended to other quantities of interest.

Reviewer's comments
We thank the editor and the reviewers for their time and effort to handle our manuscript. We revised our manuscript according to their recommendations.

Reviewer summary
The work is very nice in an area sorely in need of a quantitative approach. It is a sound approach to a very complex problem. It is original in that they are the first to put this level of detail into the system, and it is valid and significant to the field of dairy cattle nutrition and reproduction.

Recommendations to authors
Abstract. The line 'simulations confirm that an appropriate nutritional... ' really should state 'simulations describe realistically the effects of nutritional glucose supply on the ovulatory cycle of lactating dairy cattle' (no model can confirm what happens in reality, it is the other way around).
That is right, we edited the text accordingly. Page 4, line 29/30 or 50/51: "none of these models....in dairy cows. Excepting the model of McNamara and Shields, which contains many of the same elements herein, as it was based not only on Molly but on our BovCycle model (ref ); however the McNamara and Shields model did not include the full reproductive process model." (I think the authors as well as the authors of the referred to article (McNamara and Shields) should have a more in depth reference, as the efforts are closely related and complementary).
We agree and added a more detailed discussion about the model by McNamara and Shields and its relation to our modeling effort in the introductory section.
Lines 50-51. It is not clear what the 'exact value' refers to here. I think they need to be specific in that 'as long as the initial values are within this range, they do not affect the performance of the model' (which is kind of neat, as that is how reality works too).
We followed this suggestion. Page 11, line 179. All organs use glucose, except adipose tissue which cannot directly convert glucose to fatty acids.
Thank you for the correction, we adopted it. Line 286. I don't think the sole purpose was to examine the effect of glucose in the feed, but to describe the role of supplied glucose in affecting the ovulatory cycle, so it should be stated somehow like that.
In fact, it was already stated like this in the manuscript, but the sentence structure was misleading. Hence, we reformulated the sentence.
It is interesting, but not surprising, that the graphs look a lot like the simpler graphs supplied by McNamara and