|
|
||||||||


* Virginia Polytechnic Institute and State University, Blacksburg, 24061
Dexcel Ltd., Hamilton, New Zealand
1 Corresponding author: mhanigan{at}vt.edu
| ABSTRACT |
|---|
|
|
|---|
Key Words: model lactation dairy cow milk composition
| INTRODUCTION |
|---|
|
|
|---|
When simulating diets with varying nutrient content, milk yield responses are underpredicted, and BW responses are overpredicted as compared with the observed responses (McNamara and Baldwin, 2000; Hanigan et al., 2006, 2007). That is, the model underpredicts milk yield and overpredicts BW gain when dietary energy content is high and the reverse when low-energy diets are simulated (demonstrated in Figure 1
). These errors suggest that milk synthesis capacity is biologically regulated in response to nutritional state, causing partitioning of more energy to milk as nutritional state improves and protecting body stores at the expense of milk production as nutritional state declines. Insulin and the somatotropin axis have been observed to play a role in regulating milk synthesis (Asimov and Krouze, 1937; McGuire et al., 1995b), allowing the animal to alter mammary metabolite use to match its nutrient supply. Because these endocrines are responsive to nutritional state (McGuire et al., 1992, 1995a), failure to consider their effects on mammary synthetic capacity in Molly likely explains at least part of the observed prediction errors with respect to observed nutritional responses.
|
The objective of this work was to evaluate whether an alternative representation of mammary cell numbers, mammary enzyme activity, and the somatotropic axis would provide better predictions of milk and milk component yields and BW changes by Molly.
| MATERIALS AND METHODS |
|---|
|
|
|---|
All simulations and parameter estimations were conducted using ACSL Optimize (AEgis Technologies Group, Austin, TX) using a variable-step second-order Runge-Kutta-Fehlberg numerical integrator. A maximum integration interval was set at 0.05 d. Parameters were estimated using a Nelder-Mead simplex optimization algorithm to maximize the log-likelihood function. Residual errors were assumed to be homogeneous.
Revisions undertaken in Molly2006 were as follows.
Total Mammary Cells
The aggregated representation of mammary synthetic activity originally described by Neal and Thornley (1983) and utilized in Molly95 was replaced with a deaggregated representation of active and quiescent mammary cells based on the work of Vetharaniam et al. (2003a, b) with modifications of the representations of cells, cell division, and cell death.
In the model of Vetharaniam et al. (2003a, b), cell division was assumed to occur in binary fashion (i.e., secretory cells were derived solely from a fixed number of stem or progenitor cells), and cell division rates decayed to zero after parturition. However, the work of Dijkstra et al. (1997) demonstrated that cell division was exponential for many different species including goats (i.e., daughter cells also undergo division, generating additional secretory cells). Thus, it seems reasonable to conclude that cell division is exponential for cattle until evidence suggests otherwise.
In the model of Vetharaniam et al. (2003a, b), secretory cells were assumed to exist in active and quiescent states, with the latter subject to apoptosis. Cycling of cells to the quiescent pool was driven by reduced milking frequency, and the rate of cell senescence was affected by residence time in the inactive pool. In the Vetharaniam et al. (2003a, b) model, it was assumed that energy status of the animal did not affect cycling between the active and quiescent states or the rate of senescence; however, when the model was fitted to data derived from differing nutritional states, the rate of cell senescence was found to be significantly reduced on high-energy diets. Thus, both milking infrequently and dietary energy restriction were predicted to result in greater rates of cell death as compared with frequently milked or energy-sufficient states. These effects are manifested as a reduction in lactational persistency.
Aston et al. (1995) subjected several groups of cows to varying dietary energy inputs, which resulted in significant changes in milk yield. According to the observations of Vetharaniam et al. (2003a, b), the observed reductions in milk yield associated with low-energy diets should have resulted from greater cell death and lost productive capacity. But, cows that were energy-restricted in the first period of the study produced just as much milk in the second period as cows that were not energy-restricted. If cell apoptotic rates had been stimulated in period 1 as observed by Vetharaniam et al. (2003a, b), carryover effects should have been observed in the second period of the study. Because such effects were not observed, predicted changes in cell apoptotic rates are apparently not universally supported. As Vetharaniam et al. (2003b) observed such effects for studies conducted in New Zealand, it is unclear whether such a mechanism should be included. Given that inclusion of such a mechanism would clearly result in biased predictions for at least the observations of Aston et al. (1995), it was omitted in the current work.
Barnes et al. (1990) observed similar rates of decline in lactation yield (kg/mo) as lactation progressed for cows milked 3 times per day as for cows milked 2 times per day. These results suggest that the effects of milking frequency on cell senescence are not restricted to cells in the quiescent state as suggested by Vetharaniam et al. (2003a), because this would result in divergent lactational persistencies, which were not observed. These observations suggest that both active and quiescent cells are subject to apoptosis, with no apparent differences in the rates between the 2 pools.
Based on these observations, total mammary cells were represented as the balance of exponential cell division decaying as the animal approaches parturition and mass-action apoptosis. This balance was represented using a modified version of the model of Dijkstra et al. (1997), in which apoptosis was considered to occur in both the pre- and postparturient periods. Additionally, the equation was generalized for the entire lactation:
![]() | [1] |
where QCells = the total number of mammary cells at time t, which was expressed as DIM with preparturient DIM assuming negative values; T0 = the time of parturition (DIM = 0); Sign assumed a value of –1 for DIM < 0 and 1 for DIM
0; µDivision(T0) = the cell division rate at T0; KApoptosis = the rate of cell death at any point in time; and KDecay = the rate of decay in µ with respect to t.
In the representation of Dijkstra et al. (1997), KDecay could assume differing values pre- and postpartum. That representation was maintained herein, and the value of KDecay was set by DIM using a conditional statement.
Because data for prepartum mammary growth are not currently available for cattle, the prepartum KDecay was set to a value of 0.009 as observed by Dijkstra et al. (1997) for goats. Attempts were made to derive the postpartum KDecay and µDivision(T0), but the data were not adequate to uniquely derive both parameters. Because µDivision(T0) represents the rate of cell growth, it has little influence on the shape of the postpartum cell growth curve. Therefore, it was fixed to 0.03 as observed by Dijkstra et al. (1997), and the postpartum value for KDecay and KApoptosis was derived from observed data herein.
The default value for QCells(T0) was arbitrarily set to 792, because this was the value derived for that parameter when Molly95 was fitted to observations from animals of the New Zealand genotype. However, this value should be set to represent the genetic potential of the group of animals being simulated, and thus, accommodation was made for derivation of a separate value for the North American genotype.
Active and Quiescent Mammary Cells
The proportion of total mammary cells that were in the active state (PActive) was defined on a fractional basis. The differential describing the change in PActive with respect to t was as follows:
![]() | [2] |
where FQuiescent,Active and FActive,Quiescent = the fractional flux of cells from the quiescent pool to the active pool and the fractional flux of active cells to the quiescent pool, respectively. These fluxes were defined as follows:
![]() | [3] |
![]() | [4] |
where KActive,Quiescent and KQuiescent,Active = the rate parameters for cycling between the active and quiescent states. KActive,Quiescent was set to 0.11 based on the observations of Vetharaniam et al. (2003a), and KQuiescent,Active was set to 0.3, which was less than that derived by Vetharaniam et al. (2003a). This reduction was adopted to reflect the reference state of twice-daily milking and allow for significant increases in activity in response to more frequent milking. This change is arbitrary at this point and should be derived from observational data. Additionally, KFill = the scalar for inactivity conversion associated with udder fill as originally defined by Neal and Thornley (1983) and applied in Molly95. KFill assumes a value of 1 with continuous milking and a value of less than 1 for intermittent milking.
Furthermore, PActive at any point in time was determined by numerical integration of equation 2 from an initial starting proportion (iPActive):
![]() | [5] |
where iPActive was set to 0.75, because this closely approximated the average value assumed during a lactation simulation. This value is less than that of Vetharaniam et al. (2003a), reflecting the altered setting for KQuiescent,Active. The proportion of quiescent cells (PQuiescent) was then derived by difference:
![]() | [6] |
where PActive and PQuiescent were used to calculate the number of active (QActive) and quiescent (QQuiescent) cells:
![]() | [7] |
![]() | [8] |
Mammary Enzyme Activity
Equations 1 to 8 describe the mass of active mammary cells present at any point in time. However, to maintain continuity with the original representation of mammary enzyme in Molly95, a factor relating mammary enzyme and active cell numbers was required. Because changes in cell numbers with respect to stage of lactation account for the general increase and decline in mammary capacity in the revised model, total mammary enzyme activity was related to active cells via a scalar (PEnz,Cell) that represented the enzyme activity per active cell with modifications in activity influenced by lactation hormone (QLHor, defined below):
![]() | [9] |
where
was used to adjust sensitivity to QLHor. Additionally, PEnz,Cell was initially set to 12 to yield enzyme quantities roughly equivalent to the representation in Molly95. This parameter was subsequently derived by fitting to experimental data.
It was assumed that the observed changes in milk protein and fat content with respect to stage of lactation are driven by alterations in either the osmotic balance relative to lactose concentrations or that lactose yield per mammary cell is not constant as assumed for the other milk components. To address the problem, the maximal velocity (Vm) for milk lactose synthesis (VmGl,Lm) was described as a function of DIM using the equation of Dijkstra et al. (1997):
![]() | [10] |
where VmGl,Lm(T=0) = the Vm at parturition and was initially set equivalent to the value of the original VmGl,Lm (0.0025) as described by Baldwin et al. (1987b).
Additionally, kVm,Syn, kVm,Decay, and kVm,Deg were initially set to 0.005, 0.03, and 0.0005 and subsequently derived.
Lactation Hormone
Lactation hormone was altered from its original representation to reflect aspects of the somatotropin axis (primarily IGF-I), including a representation of the effects of photoperiod. The differential describing QLHor with respect to t was:
![]() | [11] |
where QLHor was determined from equation [11] by numerical integration starting with an initial mass of lactation hormone (iQLHor):
![]() | [12] |
where iQLHor was set to 1.0, which was defined as the reference state for Molly95. Synthesis (FLHor,Syn) and catabolism (FLHor,Deg) of QLHor were defined as follows:
![]() | [13] |
and
![]() | [14] |
where
and
were included to allow sensitivity adjustments;
was set to a value of 2 based on an appraisal of responses of QLHor to nutrient supply, and
was set to 2.97 based on a preliminary fit to the observed data. To derive
would require knowledge of the independent responses of the somatotropic axis to AA and glucose, which is beyond the scope of this work. The concentrations of glucose (CGlc) and AA (CAA) and the mass of adipose tissue (QAdipose) were represented as positive affecters of lactation hormone synthesis. The first 2 terms reflect the effects of energy and AA balance on somatotropin and IGF-I secretion (Chew et al., 1984; McGuire et al., 1992; Hatfield et al., 1999; Kobayashi et al., 2002). The adipose mass term reflects the relationship among fat mass, leptin secretion, and the subsequent effect of leptin on the somatotropic axis (Chilliard et al., 2000; McMahon et al., 2001; Morrison et al., 2001). The latter element reflects the concept of a set point or homeostasis for adiposity. Such a concept is supported by observational data (Chilliard et al., 2000). Additionally when such a concept is omitted from the model, stability issues occur, including a propensity to gain or lose excessive amounts of fat mass prior to changes in milk yield (Figure 1
). Because kAdipose essentially represents the set point for adipose mass, it must be scaled to BW, which was accomplished using a derivation of the equation of Waltner et al. (1994):
![]() | [15] |
where BWT0 = the postpartum BW and BCSTarget = the target BCS. Setting the latter to a value of 3.0 was found to yield appropriate declines and recovery in fat mass as lactation progressed.
It was found that VmLHor,Syn represented the maximal rate of synthesis and was arbitrarily set to a value of 4 to yield a synthesis rate of 1 in the reference state, and KLHor,Deg and KLHor,PP in equation 14 represented the basal rate of degradation of QLHor and the effects of daylength on degradation, respectively. Inclusion of the latter is supported by the observations of elevated IGF-I concentrations during long days, which were apparently associated with a reduction in clearance of IGF-I (Dahl and Petitclerc, 2003; Kendall et al., 2003). Consistent with the settings for the VmLHor,Syn, KLHor,Deg was set to a value of 1.
Additionally, KLHor,PP was calculated from daylength as follows:
![]() | [16] |
where KDaylength = a scalar for adjusting the magnitude of the effect. Daylength was calculated as:
![]() | [17] |
where DayofYear ranged from 1 (January 1) to 365 (December 31); Jan1toSprEq = the number of days from January 1 to the spring equinox (79 in the northern hemisphere and –101 in the southern hemisphere); and Latitude = the degrees of latitude at the location of the trial. SinDays was set to 0.017214 (calculated as 2
/365) to achieve a complete sin wave during the calendar year. No accommodation was made for manual manipulation of daylength, but such effects could easily be encoded in equation [17].
Parameter Estimation and Model Evaluations
One data set was used for initial model evaluations and subsequent parameter estimations, and a second independent data set was used for evaluations of the revised model after parameter estimations were completed. The first data set consisted of observations that were collected as part of an extended lactation trial conducted in New Zealand (Kolver et al., 2006). Multiparous Holstein-Friesian cows of North American or New Zealand genotypes were fed 0, 3, or 6 kg of concentrate DM daily throughout a 600-d lactation. Feed composition, milk yield, milk composition, BW, and BCS were assessed periodically throughout the trial. Blood glucose, NEFA, and urea N were also assessed at various points during the lactation. Intake was calculated using the equation of Holmes et al. (2002). The second data set consisted of the observations of Aston et al. (1995). Observed data included DM intake, diet composition, milk yield, milk composition, and BW.
For fitting purposes, observed dietary composition and intakes were averaged by week and treatment group and used as inputs to the model. Full nutrient inputs required by the model were derived from proximate analyses as described by Hanigan et al. (2005). Where needed, observed nutrient composition of ingredients was supplemented with NRC (2001) tabular values to provide a full input data set.
Observed BW and BCS at calving were used to define initial parameters for the model directly, or in the case of BCS, using the equations of Waltner et al. (1994). The model was evaluated against weekly mean observations by treatment.
Residuals were calculated as observed minus predicted, and root mean square prediction errors (RMSPE) and a decomposition of those errors were calculated from the residuals as previously described (Roseler et al., 1997). Slope bias was determined as a function of DIM.
| RESULTS AND DISCUSSION |
|---|
|
|
|---|
|
|
|
|
|
The estimated postpartum KDecay was found to be 0.44 ± 0.19. The relatively great standard deviation of the estimate is likely due to the rapidity of apparent cessation of mammary cell growth after parturition. A value of 0.44 results in cessation of cell growth by 6 DIM with a 5% increase in cell numbers after parturition. Because the difference in the growth curve is minimal for parameter estimates greater than 0.3, there is marginal ability to define the parameter (see Figure 4
). Dijkstra et al. (1997) estimated a value of 0.141 for KDecay when fitting to mammary DNA data derived from goats. Such a value is indicative of 18 d of postpartum cell growth and a 17% increase in total mammary DNA over that present at parturition. It is not clear whether the apparent differences in cow and goat data result from the method of derivation or reflect species differences.
|
It was found that KApoptosis was 2.25 x 10–3 for cows of the New Zealand genotype and 1.47 x 10–3 for cows of the North American genotype. Both estimates were less than the value of 4.0 x 10–3 observed by both Val-Arreola et al. (2004) for dairy cows and Dijkstra et al. (1997) for goats.
These 4 parameters define the shape of the lactation curve and are thus the most likely to be affected by parity, particularly primiparous vs. multiparous. However, because the data used for parameterization only included multiparous cows, the resulting parameter estimates should be evaluated with primiparous data to determine whether an alternate primiparous parameter set is warranted.
The changes undertaken in the representation of milk lactose resulted in milk composition that more accurately matched the patterns typically observed for dairy cattle. Milk fat and protein content were found to be elevated in early and late lactation and less at peak lactation, resulting in greater energy output in early and late lactation and reduced energy output at peak lactation when expressed per unit of milk volume. This contributed to greater accuracy of BW and BCS predictions.
Mammary enzyme activity was related to the predicted somatotropin axis as represented by the effects of
in equation [9], allowing greater changes in predicted milk component yields in response to nutrition than predicted by Molly95. Although the modified representation of endocrine effects improved model accuracy with respect to the observed data (see below) and thus appear to be consistent with the hypothesized mechanism of action and observed production responses to endocrine infusions, the effects of nutritional state and endocrine profiles on enzyme activities are not generally supported by more invasive measurements of mammary enzyme activity (Sorensen and Knight, 2002; Norgaard et al., 2005).
Although the existence of a set point or homeostatic point for adiposity has not been conclusively demonstrated, the hypothesis is generally consistent with observational data, including the positive relationship between adiposity and leptin secretion and the negative effect of leptin on intake (Chilliard et al., 2000). Because leptin is positively correlated with somatotropin secretion and somatotropin influences milk yield (Chilliard et al., 2000; McMahon et al., 2001; Morrison et al., 2001), the representation seems appropriate. In support of inclusion of such a concept, removal of the affects by altering the setting for
to a value of 0 results in a 33% reduction in the log-likelihood function derived from predictions of BW, BCS, milk yield and composition, and blood metabolite concentrations. The major changes were associated with increases in RMSPE of more than 50% for BW, BCS, and milk yield and an increase of 20% for predictions of blood glucose. Perhaps more importantly, slope bias for BW, BCS, milk yield, and blood glucose increased by 17, 32, 20, and 12 percentage units, respectively, indicating severe systematic prediction bias. And this degradation in accuracy was evident at
settings of 2 and 1 with log-likelihood reductions of 12 and 3%, respectively, and corresponding increases in RMSPE and the percentage of prediction errors associated with slope bias for body mass, BCS, and milk yield. Additional work is required to refine the estimate for
, which will be addressed in future work; however, assuming the current structure of Molly is appropriate, the necessity of inclusion of a set-point concept appears to be well-supported.
Photoperiod effects were also observed, as evidenced by a positive estimate for KDaylength, which is consistent with previous observations (Dahl and Petitclerc, 2003; Kendall et al., 2003). The variable milking interval portion of the model still requires parameterization, which is planned for future work. Additionally, the mammary cell growth curves for primiparous animals could be different than for the multiparous animals used for parameterization. If so, additional parameterization work would be required to derive the appropriate settings for the younger animals, because they were not represented in the current work.
Model Prediction Accuracy
Root mean square prediction errors associated with Molly2006 were generally halved for predictions of milk yield and composition, blood metabolite concentrations, and BW and BCS relative to those from Molly95. Exceptions were milk lactose content, blood NEFA concentrations, and BUN concentrations (see Table 3
). Because both models assume constant milk lactose content, it is not surprising that those predictions did not improve. And given that both models were fitted to the data, improvements in RMSPE reflect the changes in model coding and thus suggest that the changes were beneficial. More specifically, use of equation [10] provided the desired relationship among the milk components, resulting in high milk fat and protein in early and late lactation and low contents at peak lactation (Figure 3
), which significantly reduced prediction errors.
The proportions of RMSPE associated with mean bias were substantially reduced for the content of lactose, protein, and fat, and slope bias was reduced for protein content and milk protein and fat yields, reflecting the inherent problems with the Molly95 representation of milk component percentages.
Residual errors plotted against DIM for predictions of milk and milk lactose yields are presented in Figure 5
, and errors for milk protein and fat yield predictions are presented in Figure 6
. There were no apparent patterns to milk yield or milk component yield residuals by dietary treatment, supporting the hypothesis stated previously that rates of secretory cell apoptosis were not affected by dietary status (i.e., reductions in milk yield associated with nutritional insufficiency apparently result from changes in enzyme activity per cell and not the number of mammary cells). Such a finding is consistent with the observations of Sorensen and Knight (2002) but inconsistent with the observations of Vetharaniam et al. (2003a, b). The difference in the current findings relative to those of Vetharaniam et al. (2003a, b) likely reflects the more extensive representation of mammary enzyme activity in the current model as compared with that in the model of Vetharaniam et al. (2003a, b).
|
|
However, gestational nutrient use would only explain a portion of systematic errors in milk fat production. Rates of milk fat synthesis are underpredicted from 0 to 50 DIM, overpredicted from 50 to 200 DIM, and overpredicted again from about 550 DIM until the end of lactation. Overpredictions after 550 DIM may be due to gestational requirements, but obviously that is not the cause of errors during the first season when animals were not pregnant.
Underpredictions of milk fat yield and content in early lactation are likely due to underpredictions of blood NEFA concentrations and thus fatty acid removal and use for milk fat synthesis (Figure 7
). In the model, NEFA and triacylglycerol (TAG) are considered as a common pool. Because TAG concentrations were not measured in the extended lactation study, it cannot be fully determined whether changes in blood fatty acids were appropriate with respect to observed values. If the representation of fatty acids was deaggregated to represent NEFA and TAG independently, the rise in NEFA in early lactation could be used to drive the elevations in milk fat. Because NEFA concentrations fall, they would represent a lesser proportion of total fatty acids due to the more stable contribution of TAG to the total pool. Thus, it may be possible to drive milk fat output up in early lactation without significantly reducing output at peak lactation.
|
As might be expected, improved predictions of milk component yields and removal of systematic bias improved the accuracy of predictions of energy partitioning, and this is reflected in reductions in RMSPE for BW and BCS (Table 3
). In particular, reductions in errors of milk fat would be expected to have the greatest contribution to improvements in predictions of body energy stores. Such errors are consistent with the observation of mean and slope bias for Molly95, which was reduced for Molly2006. Although predictions were improved over those of Molly95, systematic bias in predictions is still apparent (Figure 8
) for both BW and BCS. In particular, the rate of BW and BCS loss in early lactation is underpredicted, resulting in overpredictions of score at peak lactation. Although the model appears to underpredict weight recovery as lactation progresses, the prediction errors for BCS do not change, suggesting that after peak lactation, the model predicts fat metabolism with some accuracy but exhibits bias in predictions of lean mass. The latter could certainly be associated with fetal growth. Some systematic bias in both BW and BCS predictions was associated with dietary supplementation rate. Simulations of the unsupplemented animals resulted in underpredictions of BW and BCS for most of the lactation.
|
|
|
In summary, predictions of milk and milk components, BW, and BCS were significantly improved by adopting a derivation of the representation of mammary cells and mammary enzymes described by Dijkstra et al. (1997) and Vetharaniam et al. (2003a, b) and by altering the representation of the somatotropin axis. Representation of mammary secretory cell loss as a function of time and unaffected by nutritional state was consistent with the observed data. The model appears to now respond appropriately to dietary energy inputs with respect to at least milk, milk lactose, and milk protein yields, and prediction errors for BW and BCS are significantly improved. Additional work on the representation of blood fatty acids and the effects of this pool on milk fat and body fat stores is required.
| ACKNOWLEDGEMENTS |
|---|
|
|
|---|
Received for publication January 12, 2007. Accepted for publication April 13, 2007.
| REFERENCES |
|---|
|
|
|---|
This article has been cited by other articles:
![]() |
P. C. Beukes, C. C. Palliser, K. A. Macdonald, J. A. S. Lancaster, G. Levy, B. S. Thorrold, and M. E. Wastney Evaluation of a Whole-Farm Model for Pasture-Based Dairy Systems J Dairy Sci, June 1, 2008; 91(6): 2353 - 2360. [Abstract] [Full Text] [PDF] |
||||
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |