Quantifying Shape of Lactation Curves , and Benchmark Curves for Common Dairy Breeds and Parities

The MilkBot® (DairySight LLC, Argyle, NY; http:// milkbot.com) lactation model provides a means of quantifying both shape and magnitude of lactation curves as a set of parameter values, each of which is associated with a single aspect of lactation curve shape. Lactation data may be fitted to the model to summarize a lactation as a set of parameter values which summarize the lactation as a whole. The scale parameter controls magnitude without changing the shape of the curve, the ramp parameter controls steepness of the post-parturient rise in milk production, the decay parameter controls the rate of late lactation decline, and the offset parameter defines a theoretical offset between the start of milk production and calving. The decay parameter is easily re-expressed mathematically as persistence to quantify the rate of decline in production after peak milk. Time and quantity of peak milk, or production for any day or period in the lactation may be calculated directly from parameter values. Aggregate normal lactation curves for mean and median milk production of Holstein, Jersey, Crossbred, Guernsey, Ayrshire, and Brown Swiss dairy cattle, stratified by parity, are calculated from a DHIA data set collected from 2005 through mid-2008 and covering over six million lactations and 51 million milk weights, mainly from farms in the eastern United States. These constitute benchmark curves that may be used as standards for normal milk production, or to quantify changes in normal productivity over time or with respect to other variables, or in econometric studies.


Introduction
Milk production in dairy cows is a sensitive measure of cow health, as well as being the most important determinant of a dairy farm's income.Almost everything controlled in managing a dairy herd, including health, feeding, environment, breeding, and welfare, will influence daily milk production.This suggests that by developing more precise and accurate tools for measuring changes in milk production, we can better measure the influence of many other variables.This includes measuring the herd response to an intervention, the key metric of production medicine.
A lactation model can correct for the confounding effect of the normal rise then fall of production after calving, i.e. the normal lactation curve.The mathematical model, which can be customized using a priori information about the animals studied, can be used in many ways, particularly in computer software or spreadsheet applications to project future production or compare actual to projected production.Various curve-fitting algorithms may also be used to fit a set of observed data points to the model, yielding parameter values, which constitute a condensed quantitative summary of the data set.
Many lactation models have been proposed and compared over the past half-century, 4,5,7,8,10 but none have achieved widespread acceptance outside a few specialized applications which are directed primarily at improving estimates of actual production from incomplete data sets.The MilkBot ®a model has close mathematical similarities to some of these earlier models, such as the Mitscherlich-Exponential model proposed by Rook 8 in 1993, but differs in that the model is derived from a theoretical-mechanistic hypothesis, 2 leading to parameters which can be interpreted both in terms of the effect that they have on the curve and in terms of the mechanistic hypothesis.This is of critical importance to the interpretability of fitted parameter values, which become metrics in their own right of the distribution of milk production within a lactation, or lactation curve shape.
A few metrics have been proposed to measure the distribution of production within a lactation that do not rely on an explicit lactation model.The term "persistency", for example, has been used to quantify the rate of decline in milk production in the later part of lactation.Unfortunately definitions of persistency vary widely, 1,6,9 and there is no standard definition in use.Similarly "peak milk", meaning the highest daily production, while strictly a measure of magnitude, can be used to quantify shape when compared to cumulative production.It can also be observed that some lactations rise more steeply after calving than others.One attempt to quantify this aspect of curve shape is "time to peak milk", which has the heavy disadvantage of being dependent on frequency of data collection (which is often monthly).Difficulties in calculation and accuracy of these measures, along with lack of clear definitions in some cases, have kept them from wide application. 3Currently there is little understanding of normal variation in lactation curve shape between individual animals and groups aside from the very basic observation that persistency is typically higher in primiparous animals.

Materials and Methods
The MilkBot ® Model (equation 1, below) predicts daily milk yield, Y(t) as a function of time (t).Four parameters, a (scale), b (ramp), c (offset), and d (decay), control the shape of the curve.The constant e is Euler's number (i.e. the root of natural logarithms, approximately 2.718).

Interpretation of MilkBot ® Parameters
Each of the parameters (a,b,c,d) in equation 1 dominates a particular aspect of lactation curve shape, and has a descriptive name related to that effect.
Parameter a is the scale parameter.It is a simple multiplier, which determines the overall magnitude of milk production.It can be expressed as pounds/day, kilograms/day, or similarly.This is the theoretical maximum daily milk yield, which is approached by actual peak production as ramp, and offset values approach zero (i.e. a lactation which peaks on the day of calving), or as decay approaches zero (infinite persistence).It must be a positive number.Changing the model to a different unit of measure for milk output only requires applying the appropriate conversion to the scale parameter, and all other parameter values remain unchanged.
Parameter b is the ramp parameter, controlling the rate of rise in milk production in early lactation.Ramp values are time, normally in days.Smaller ramp values imply faster creation of productive capacity and a steeper rise in early lactation.Ramp must be a positive number.
Parameter c is the offset parameter, and has relatively minor influence on the model.It represents the offset in time between calving and maximal growth rate of productive capacity.Its effect is slight, except in the first few days of a lactation.Offset values are time (days), and indicate the time of maximal creation of productive capacity.They may be positive, negative, or zero.
Parameter d is the decay parameter, controlling the loss of productive capacity, and analogous to the first-order decay constant common in pharmacokinetics.Decay is inverse-time (days -1 ).It should be constrained to positive values under normal circumstances, though it can be argued that there might be situations in which negative decay could be biologically feasible.
A first-order decay parameter, such as the decay parameter, may easily be converted to a half-life, which allows us to address the existing concept of lactation persistency.By this re-expression (equation 2, below), the fourth parameter becomes persistence, correspond- ing to the time it would take for production to drop by half, if we were to ignore the growth side of the model.Persistence is close to the actual half-life for productive capacity in late lactation.

Peak and Cumulative Production
Mathematical manipulation of equation 1 allows calculation of some useful results.By setting the derivative equal to zero we can calculate the day (tPeak) when milk production peaks, and then yPeak = Y (tPeak) gives peak milk production by substitution in equation 1.

3)
Cumulative production between two dates t 1 and t 2 can be calculated by integrating Y(t), substituting values for t 1 and t 2 in the integral, and finding the difference.M305 is the cumulative milk yield between calving and day 305 of the lactation, calculated in this way.

Milk Production Data
A large set of DHIA milk production records b from January 2005 through June 2008 was imported into a computer database.c This included 51,180,569 monthly milk weights from 6,459,942 lactations in 17,013 herds.Median and mean milk weights were calculated for each days-in-milk (DIM) value between six and 400 days, for parities 1, 2, and 3 of the six most common breeds.Herds contributing records are primarily in the eastern portion of the US.Data was not otherwise edited except that records from unspecified or minor breeds was discarded.

Curve Fitting
Data for each group was fitted to the MilkBot ® model using an implementation of the Levenberg-Marquardt algorithm d to minimize mean square error (MSE).In all cases, the fitting process converged readily and was insensitive to minor changes in input data.Fitted parameter values and MSE for each aggregate curve are reported in Table 1.An internet supplement e to this paper supplies high-resolution graphs showing each aggregate curve plotted with the data points to which it was fitted.Final MSE values varied between 0.32 lb (0.15 kg) and 9.9 lb (4.5 kg), depending mainly on the number of lactations in the group.

Fitted Parameter Values
Table 1 shows fitted parameter values for mean and median daily milk production for five dairy breeds, and parities 1, 2, and 3. Mean and median curves are similar in all cases, as are second and third parity lactations within each breed.There are larger differences in curve shape and magnitude between breeds and between first and later parities.
Scale rises with parity in all breeds, and varies considerably between breeds with Holsteins having greatest scale values followed by crossbred cows, then other breeds.It may be noted that the crossbred category is likely less homogeneous than the true breed categories.Scale values are a general measure of the productivity of the breed, but not tied to a particular lactation length.It remains to be seen whether scale is better or worse correlated with genetics or management than other measures of productivity, such as M305.
Ramp values seem to be lower for second parity lactations (that is a steeper rise in post-parturient production) than the other parity groups.Holstein lactations generally rise more slowly (greater ramp) than other breeds.A probable but unproven hypothesis is that ramp values will be influenced heavily by transition cow management.
Offset values cover a very narrow range, but Jerseys seem to have higher offset than other breeds.Since offset represents the difference between the recorded birth of the calf and the physiological start of lactation, one possible explanation would be the excellent calving ease expected from the Jersey breed.This is because a difficult birth could delay the calving date relative to the physiological onset of lactation, which would decrease the offset value for the lactation.This is only a speculative rationalization of a small observed difference, but it seems possible that a very large database could show differences in calving ease through the offset parameter.
Decay values show the expected pattern of being considerably lower (i.e. higher persistence) in lower parities.Jersey cattle seem to have the lowest decay (highest persistence), and Ayrshire cattle the highest decay or lowest persistence.Use of BST might have a confounding effect on decay that is difficult to predict.
Since data are aggregated before fitting, it is not possible to tell whether differences result from different distributions between breed populations, or consistent differences between individuals of different breeds.It is also possible that differences between breeds found here may be confounded by patterns of popularity of breeds in different climatic zones, or other factors, and so reflect differences in environment or management in addition to true breed characteristics.

Persistence
Calculation of persistence as a half-life to replace current measures of persistency has advantages and disadvantages.The largest disadvantage is that it is a subtle re-expression of an existing concept, so will introduce some confusion of terms.There are already multiple conflicting definitions of persistency in use, however, some are poorly defined or applicable to only a specific section of the lactation curve.By defining a metric that is an attribute of the lactation as a whole rather than a loosely defined post-peak zone, some problems of sampling bias are avoided, and the metric benefits from a more biologically-oriented derivation.

Graphed Curves
Normal median production for Holstein, Jersey, and crossbred cattle is graphed in Figure 1 (for first parity lactations) and Figure 2 (second parity).The parameter values used in plotting the curves are taken from Table 1.  Figure 3 plots data points for Holstein mean daily milk along with the fitted curve, showing an excellent general fit, and also some interesting anomalies.Similar patterns are visible in graphs for other groups (not shown, but available online).There appears to be a small inflection point in the observed data just before 300 DIM, which cannot be matched by the MilkBot ® mathematical model.In other words, persistency seems to improve slightly near 300 DIM.This effect is small, and could be an artifact of population dynamics rather than individual lactations.For example a subpopulation of cows with extended lactations of higher persistence would have a larger relative effect as cows in the main population become pregnant and are dried off.These hypothetical high-persistence lactations might be from infertile or do-not-breed cows, and possibly the use of BST.
Also there is an anomalously low data point at 30 DIM, which may be due to data errors in reporting of calving dates.If calving dates were occasionally entered as the previous test date rather than the actual calving date, it would have the effect of biasing points near 30 DIM downwards.That is, some of the cases nominally at 30 DIM may actually be earlier in their lactations, so producing less than they would when they do reach 30 DIM.This conjecture is supported by the observation that there are about 10% more data points at 30 DIM than surrounding DIM values, suggesting about a 10% error rate for calving date.

Benchmark Milk Yield
Parameter values reported in Table 1 can be used with the MilkBot ® model to predict milk for a particular breed and parity at any value of DIM, though values for greater than 400 DIM are speculative since the model was not fitted past 400 days.For example, using parameter values for median production of second parity Holsteins [a=109 lb or 49.4 kg (scale), b=27.2 days (ramp), c= -3.8 days (offset), and d=0.00215 days -1 (decay)] in equation 1, we can calculate median second parity Holstein production at 100 days at 86.9 lb (39.4 kg).Similarly we can calculate that in the week between 100 DIM and 107 DIM, production is expected to drop by 1.1 lb (0.50 kg).Formulas for these calculations, formatted for use in spreadsheet software, are shown in Appendix 1. Table 1 shows that this corresponds to a 305-day total of 23,165 lb (10,507 kg).Persistence of mean milk for that group is 322 days, with a peak of 91.3 lb (41.4 kg) at DIM=56 days.This is based on 12,675,546 data points from 1,583,397 second parity Holstein lactations.

An Example Simulation
It is sometimes asserted that increasing peak milk production in a herd will increase total lactation production by some multiplier, typically around 250 lb (113 kg) of milk per incremental pound at peak.The exact increment will depend on the shape of the lactation curve, and the questionable assumption that whatever management changes altered peak milk will not alter lactation curve shape.As an example simulation, the standard curves in Table 1 were used to plot the relationship between a change in peak milk and annualized milk per cow, for several breeds and parities.Simulations of this type could easily be included in a spreadsheet or other computer software.
Total milk over a lactation of any length can be calculated similarly to the way M305 is calculated, by substituting the actual lactation length into equation 4 wherever the number 305 occurs.If that total is divided by the sum (lactation length plus dry period length), we get average milk per day during the full period of the lactation, then multiplying by 365 gives annualized milk.If this is then divided by peak milk, it can be interpreted as the expected change in annualized milk for a small incremental change in peak milk, assuming the lactation curve shape does not change.This is plotted in Figure 4, for selected breeds and parities.
The simulation indicates a higher effect for an incremental change in peak milk for parity 1 (about 285-fold) than parity 2 (about 245-fold).This could have been predicted from the flatter lactation curve of heifers.There is little difference between breeds, and the incremental effect is lower for longer lactations in cows, but minimally so for heifers.Results also depend on dry period length, which was fixed at 50 days for this simulation.This exercise is meant as an illustration of the sort of calculations that are possible from a lactation model, but not a full exploration of the scenario.

Internet Supplement
An internet supplement to this paper, available through the members area of the AABP website e includes individual graphs for each parity and breed (like Figure 3), and also graphs similar to Figures 1 and 2 for Guernsey, Brown Swiss, and Ayrshire breeds.

Conclusion
The following equation is proposed: to describe normal daily milk production as a function of DIM.Specific effects that the parameters a (scale), b (ramp), c (offset), and d (decay) have on lactation curve shape are described.The decay parameter is easily transformed into a measure of persistence.The mathematical model may be used with fitting algorithms to generate quantitative measures of shape and magnitude of lactation curves.This MilkBot ® model was fitted to a large DHIA data set to generate parameter values that describe aggregate curves for major breeds and parities.These constitute benchmark curves, which are easily incorporated into computer applications or spreadsheets as standards of normal production.
MilkBot ® parameters constitute a precise professional language suitable for describing lactation curve shape.With practice, it is not difficult to visualize a set of parameter values as the 2-dimensional shape of a lactation curve, or to estimate parameter values by eye from a lactation graph.

Figure 3 .
Figure 3. Holstein second parity mean daily milk with fitted lactation curve.

EndnotesaFigure 4 .
Figure 4. Effect of an incremental change in peak milk on annual milk assuming a 50-day dry period.

Table 1 .
Fitted MilkBot ® parameter values for major breeds and parities.