Application of growth models to describe the lactation curves for test-day milk production in Holstein cows

ABSTRACT

In order to describe the lactation curves of milk yield traits, six standard growth models (Brody, logistic, Gompertz, Schumacher, Von Bertalanffy and Morgan) were used. Data were 911,144 test-day records for unadjusted milk yield , 4% fat-corrected milk yield and energy-corrected milk yield from the first lactation of Iranian Holstein cows, which were collected on 834 dairy herds in the period from 2000 to 2011. Each model was fitted to monthly production records of dairy cows using the NLIN and MODEL procedures in SAS. The models were tested for goodness of fit using adjusted coefficient of determination, root means square error, Durbin–Watson statistic, Akaike's information criterion (AIC) and Bayesian information criterion (BIC). The Morgan model provided the best fit of the lactation curves due to the lower values of AIC and BIC than other models, but the Brody model did not fit the production data as well as the other equations. The results showed that the Morgan equation was able to estimate the time to the peak and peak yield more accurately than the other equations. Overall, evaluation of the different growth equations indicated the potential of the nonlinear functions for fitting monthly productive records of Holstein cows.

KEYWORDS:

• Milk yield
• growth model
• lactation curve
• dairy cow

1. Introduction

The mathematical description of the temporary evolution of milk production in ruminants kept for milk production indicates one of the most important applications of mathematical equations in animal science (France & Thornley 1984; Pulina & Nudda 2001). Several reasons can be found for the requirement of a mathematical modelling of the lactation pattern. Mathematical models of the lactation curve and, in general, of the mammary gland represent a valuable tool for basic research studies aimed at enhancing the scientific knowledge of complex physiological mechanisms that underlie the process of milk secretion (Dimauro et al. 2007). Lactation curves may be used by physiologists, nutritionists and other research to mimic the lactation process and to evaluate the relationships existing between secretary cells, hormones, energy supply and environmental factors affecting the milk production process.

The term lactation curve refers to a graphical presentation of the ratios between milk production and lactation time starting at parturition (Papajcsik & Bodero 1988). Functions that describe milk production in time can be very applicable in genetic breeding programmes, herd nutritional management, decision-making on the culling cows and simulation systems of milk production over the lactation. The lactation curve is also important because its wide characterization of the animal production throughout lactation allows estimating the peak yield (PY), persistency of lactation and days in milk (Ferreira & Bearzoti 2003). There are different mathematical models describing lactation curves in dairy cows, from the more empirical equations that relate input to output statistically with little consideration of the biology of lactation (e.g. Wood 1967; Rook et al. 1993), to the more mechanistic ones which describe the lactation curve based on the biology of lactation (e.g. Dijkstra et al. 1997; Pollott 2000).

Growth functions which have been used to describe growth in various animal (Lopez et al. 2000) and microbial (Lopez et al. 2004) species can also be used as simple empirical equations to describe cumulative yield of milk over lactation by analogy with growth of body tissues and microbial populations because the lactation curve is similar to the plot of growth rate (daily weight gain) against time (Thornley & France 2007). The use of mathematical growth models provides a good way of condensing the information contained in a data series into a few parameters with biological meaning, in order to facilitate both the interpretation and the understanding of the phenomenon (Fitzhugh 1976).

Milk production per cow has been considerably improved within the last decades, but the increase in milk yield has been accompanied by a higher occurrence of health and fertility problems worldwide (Ghavi Hossein-Zadeh 2013). The negative energy balance in early lactation is considered to play a decisive role in this context. Most of the health problems in highly selected dairy cows occur at the beginning of lactation and have been related to the extent and duration of the energy-deficit post-partum (Collard et al. 2000; Ingvartsen & Andersen 2000). Cows require sufficient energy intake to support their physiological functions, namely to produce milk and grow while reproducing and staying healthy (Banos et al. 2006). However, selection for milk traits has resulted in a partial shift of available energy towards milk production at the cost of other biological pathways (Collard et al. 2000). There is a general agreement that energy balance should be included into future breeding programmes, but direct measurements are costly and difficult to be implemented in practice (Hüttmann et al. 2009; Ghavi Hossein-Zadeh 2012). Thus, there is a large interest in low-cost traits describing the energy status in an appropriate way. Energy-corrected 305-d milk yield (ECM) and 4% fat-corrected milk yield (FCM) account for the true protein content of milk and determines the amount of energy in the milk based upon milk, fat and protein. The ECM and FCM could be considered as criteria for the energy status of dairy cows (Ghavi Hossein-Zadeh 2012).

Six standard empirical mathematical models (Wood, Dhanoa, Sikka, Nelder, Hayashi and Dijkstra) were used to fit 305 d standard lactations of Iranian Holsteins in the previous study (Ghavi Hossein-Zadeh 2014). But, the aim of the current study was to compare the performance of six standard growth functions (Brody, logistic, Gompertz, Schumacher, Von Bertalanffy and Morgan) by fitting these equations to monthly records of unadjusted milk yield, 4% FCM and ECM for standard 305-d lactation in Iranian Holstein cows.

2. Materials and methods

2.1. Data

Data were 911,144 test-day records for unadjusted milk yield (MY), 4% FCM and ECM from the first lactation of Iranian Holstein cows. Data were recorded on 834 dairy herds in the period from 2000 to 2011 by the Animal Breeding Center of Iran. Outliers and out-of-range productive records were deleted from the analyses. Records from DIM <5 and >305 d were eliminated. Records were also eliminated if no registration number was present for a given cow. Ages at first calving were required to be between 20 and 40 months and individual daily milk production should be between 3 and 90 kg. Also, fat and protein percentages should be in a range from 1 to 9%. The following formula was used for the calculation of 4% FCM of each cow (Gaines 1928): FCM = [(0.4 × kg milk) + (0.15 × kg milk × fat%)]. Also, the yield of energy-corrected milk (ECM) was calculated by the following formula (Sjaunja et al. 1990):

2.2. Lactation curve models

The non-linear growth equations used to describe the lactation curves for milk traits are presented in Table 1. The Brody, logistic, Gompertz, Schumacher, Von Bertalanffy and Morgan equations are growth functions (Thornley & France 2007). A typical lactation curve rises to a peak before falling away, which is the same trajectory mapped by the slope of a sigmoidal growth function. Therefore, growth functions expressed as a function of time have potential application as lactation equations and are simple to use (Fathi Nasri et al. 2008). For all models, PY was assumed as the maximum test-day MY, FCM and ECM, and peak time (PT) was accepted as the test time, at which daily yield was maximum. Predicted 305-d MY were obtained for every model using the following equation (Vargas et al. 2000) with substitution of by the corresponding model equation (Table 1):where 305TY denotes predicted 305-d MY and y(t)  represents MY at day t (5, … , 305) estimated by corresponding lactation models.

Table 1. Growth equations used to describe the lactation curve of Holstein cows.

Equation	Functional form
Brody
Logistic
Gompertz
Schumacher
Von Bertalanffy
Morgan

Notes: y represents milk yield traits; PY, the maximum value for milk yield traits; PT , the PT; a, b, c and k are parameters that define the scale and shape of the lactation curve; t , the time from parturition.

The persistency measures used in this study were ratios between different parts of the lactation (P_2:1, P_3:1 and P_Weller). The P_2:1 and P_3:1 were proposed by Johansson and Hansson (1940). P_2:1 and P_3:1 are the ratios between the MYs of the second and third 100 days of lactation, respectively, and that of the first 100 days. Also, Weller et al. (2006) defined milk persistency as estimated milk production at 180 day after peak divided by estimated peak production in percent as follows:where PROD(270) and PROD(90) are milk production at 270 and 90 days in milk, respectively.

2.3. Statistical analyses

Each model was fitted to monthly productive records of dairy cows using the NLIN and MODEL procedures in SAS (SAS Institute 2002) and the parameters were estimated. When nonlinear functions were fitted, the Gauss–Newton method was used as the iteration method. The models were tested for goodness of fit (quality of prediction) using adjusted coefficient of determination (), residual standard deviation or root means square error (RMSE), Durbin–Watson statistic (DW), Akaike's information criterion (AIC) and Bayesian information criterion (BIC).

was calculated using the following formula:where is the coefficient of determinationTSS is total sum of squares, RSS is residual sum of squares, n is the number of observations (data points) and p is the number of parameters in the equation. The coefficient of determination lies always between 0 and 1, and the fit of a model is satisfactory if is close to unity.

RMSE is a kind of generalized standard deviation and was calculated as follows:where RSS is the square root of residual sum of squares, n is the number of observations (data points) and p is the number of parameters in the equation. The best model is the one with the lowest RMSE.

DW was used to detect the presence of autocorrelation in the residuals from the regression analysis. The DW ranges in value from 0 to 4. A value near two indicates non-autocorrelation; a value towards 0 indicates positive autocorrelation; a value toward 4 indicates negative autocorrelation. DW was calculated using the following formula:where is the residual at time t, and is residual at time t − 1.

AIC was calculated using the following equation (Burnham & Anderson 2002):A smaller numerical value of AIC indicates a better fit when comparing models. Also, BIC was calculated as follows:A smaller numerical value of BIC indicates a better fit when comparing models.

3. Results and discussion

Estimated parameters of non-linear growth equations for the first parity dairy cows are shown in Table 2. Also, goodness-of-fit statistics for the six growth functions fitted to average standard curves of MY, FCM and ECM are presented in Table 3. Morgan equation provided the lowest values of RMSE, but Brody and Von Bertalanffy models had the greatest values. For MY, Morgan equation provided the lowest value of DW, but Brody and Von Bertalanffy models had the greatest values. For FCM, Morgan, logistic, Gompertz and Schummacher models provided the lowest value of DW, but Brody and Von Bertalanffy models had the greatest values. For ECM, Morgan, Gompertz and Schummacher models provided the lowest value of DW, but Brody and Von Bertalanffy models had the greatest values (Table 3). values were similar among models for all traits. AIC and BIC values were the lowest for Morgan model and were the greatest for the Brody equation. DW values obtained from fitting growth models in this study were moderate and this indicated that there was positive autocorrelation between the residuals for all non-linear models. Therefore, it appears that residual autocorrelation provided no problem, or only a very slight one. Positive autocorrelation is a serial correlation in which a positive error for one observation increases the chances of a positive error for another observation (Ghavi Hossein-Zadeh 2014). Comparing the current results with the pervious study of the same author, Ghavi Hossein-Zadeh (2014) indicated that non-linear growth functions used in the current study provided better fit of the lactation curve for first lactation MY than the standard non-linear models which were routinely used for describing the lactation curve. The reason for this conclusion was the lower AIC values of growth models compared with the values obtained for standard lactation models (Ghavi Hossein-Zadeh 2014).

Table 2. Parameter estimates for the different lactation equations of the first parity dairy cows (standard errors are in parentheses).

Trait	Parameter	Model
		Brody	Logistic	Gompertz	Schumacher	Von Bertalanffy	Morgan
MY	a	31.48 (0.01)	0.00500 (0.00002)	4301.00 (20.66)	1647.50 (9.41)	31.480 (0.008)	38,057.5 (120.1)
	b	0.50 (0.01)	8327.40 (36.22)	0.00700 (0.00003)	188.00 (0.63)	−1.68 (0.02)	793.60 (3.22)
	c	–	24,264.90 (61.88)	0.00400 (0.00001)	0.01600 (0.00008)	–	1.1600 (0.0007)
	k	0.12 (0.002)	–	–	–	0.130 (0.002)	–
FCM	a	28.06 (0.01)	0.00400 (0.00002)	5962.80 (45.33)	2727.10 (23.02)	28.060 (0.008)	49,667.2 (280.4)
	b	0.29 (0.01)	10,595.60 (74.97)	0.00500 (0.00003)	279.20 (1.56)	−2.28 (0.05)	1290.90 (9.13)
	c	–	27,394.6 (127.1)	0.00300 (0.00001)	0.01000 (0.00007)	–	1.1000 (0.0008)
	k	0.150 (0.005)	–	–	–	0.150 (0.005)	–
ECM	a	27.78 (0.01)	0.00400 (0.00002)	5718.10 (41.13)	2710.80 (22.09)	27.78 (0.01)	53005.8 (315.2)
	b	0.23 (0.01)	10085.10 (66.77)	0.00500 (0.00003)	286.50 (1.56)	−2.50 (0.04)	1404.10 (10.31)
	c	–	26,830.0 (116.9)	0.00300 (0.00001)	0.01000 (0.00007)	–	1.1000 (0.0008)
	k	0.120 (0.004)	–	–	–	0.120 (0.004)	–

Notes: MY, Unadjusted milk yield; FCM, 4% Fat-corrected milk yield; ECM, energy-corrected milk yield; a, b, c and k are parameters that define the scale and shape of the lactation curve.

Table 3. Comparing goodness of fit for average standard curves of milk yield traits for Brody, Logistic, Gompertz, Schumacher, Von Bertalanffy and Morgan models.

Trait	Statistics	Model
		Brody	Logistic	Gompertz	Schumacher	Von Bertalanffy	Morgan
MY		0.95	0.95	0.95	0.95	0.95	0.95
	DW	0.95	0.91	0.90	0.89	0.95	0.88
	RMSE	7.34	7.19	7.17	7.14	7.34	7.11
	AIC	16,136,045	16,097,425	16,093,776	16,086,287	1,613,5969	16,079,399
	BIC	3,632,947	3,594,327	3,590,678	3,583,188	3,632,870	3,576,300
FCM		0.94	0.94	0.94	0.94	0.94	0.94
	DW	1.08	1.05	1.05	1.05	1.08	1.05
	RMSE	7.15	7.07	7.07	7.06	7.15	7.05
	AIC	16,088,019	16,067,598	16,066,849	16,065,307	16,088,013	16,063,194
	BIC	3,584,921	3,564,499	3,563,750	3,562,209	3,584,915	3,560,096
ECM		0.94	0.94	0.94	0.94	0.94	0.94
	DW	1.06	1.05	1.04	1.04	1.06	1.04
	RMSE	6.80	6.73	6.73	6.73	6.80	6.72
	AIC	15,995,572	15,978,511	15,977,761	15,976,391	15,995,565	15,975,020
	BIC	3,492,474	3,475,413	3,474,662	3,473,293	3,492,466	3,471,921

Notes: MY, unadjusted milk yield; FCM, 4% fat-corrected milk yield; ECM, energy-corrected milk yield; , adjusted coefficient of determination; RMSE, root means square error; DW, Durbin–Watson; AIC, Akaike information criteria; BIC, Bayesian Information Criteria

Empirical and more mechanistic equations describing the lactation curve with different features and complexity have been reported in several studies. Six non-linear growth models representing various complexities were evaluated using production data obtained from Iranian dairy herds. Comparing the predictive ability of these models allows identification of a growth function capable of characterizing and providing a better perspective on the lactation curve shape of the Holstein cows in Iran. When an appropriate model of the lactation curve is selected, selection emphasis can then be placed solely on the special parts of the lactation curve. The provision of an optimal strategy to attain a desired lactation shape through modifying the parameters of lactation model is very important. If possible, shape of the lactation curve is determined, dairy cows can be classified based on their expected lactation milk productions. Then, more appropriate feeding programmes can be assigned to each group of animals according to their requirements by taking into consideration the variations among the groups.

Observed MY, FCM and ECM increased from Day 5 of lactation to a peak a few weeks later, decreasing thereafter until Day 305. The average MY, FCM and ECM increased from 23.8, 24.3 and 24.4 kg/day at day 5 of lactation to the PY of 33.8, 29.7 and 29.2 kg/day on day 78 and subsequently declined to 27.3, 25.4 and 25.5 kg/day on day 305, respectively. Time to the peak (PT) and production at peak (PY) for average standard lactations of MY, FCM and ECM predicted by the six non-linear growth functions are shown in Table 4. Also, predicted MY, FCM and ECM lactation curves by Brody, logistic, Gompertz, Schumacher, Von Bertalanffy and Morgan equations are depicted in Figure 1. Evaluation of the lactation features showed that the Morgan equation was able to estimate PY and PT more accurately than the other functions. Time to the peak for MY, FCM and ECM was over-predicted by the logistic, Gompertz and Schumacher models and under-predicted by the Brody and Von Bertalanffy models. For MY, the Gompertz model over-predicted PY, but Brody, logistic, Schumacher and Von Bertalanffy models under-predicted PY. The Gompertz and Schumacher models were over-predicted PY for FCM and ECM, but Brody, logistic and Von Bertalanffy models were under-predicted PY.

Figure 1. Predicted MY (a), FCM (b) and ECM (c) lactation curves by Brody, logistic, Gompertz, Schumacher, Von Bertalanffy and Morgan equations.

Table 4. PT and PY for average standard lactations of MY, FCM and ECM predicted by Brody, Logistic, Gompertz, Schumacher, Von Bertalanffy and Morgan models.

Trait	Statistics	Model
		Brody	Logistic	Gompertz	Schumacher	Von Bertalanffy	Morgan
MY	PT	68	130	140	95	63	84
	PY	31.48	30.33	36.42	31.93	31.48	33.68
FCM	PT	50	115	170	111	50	81
	PY	28.06	27.39	34.84	30.90	28.06	29.24
ECM	PT	60	127	170	124	61	88
	PY	27.78	26.83	33.41	31.37	27.78	28.69

Note: MY, unadjusted milk yield; FCM, 4% fat-corrected milk yield; ECM, energy-corrected milk yield; PY, peak yield; PT, peak time.

Several authors have shown variations in the general shape of the lactation curve (e.g. Landete-Castillejos & Gallego 2000; Fathi Nasri et al. 2008; Ghavi Hossein-Zadeh 2014), the most common shape being a rapid increase after parturition to a peak a few weeks later, followed by a slow decrease until the cow is dried off. The other shape is a gradual decline from calving. Fathi Nasri et al. (2008) reported PY was under-predicted by the logistic, Gompertz, Morgan and Schumacher models and consistent with the current study, they reported PT was over-predicted by the logistic, Gompertz and Schumacher equations, but Morgan equation could predict it accurately.

Different measures of persistency and 305-d yields for average standard lactations of dairy cows predicted by various growth functions are presented in Table 5. The Gompertz model predicted the greatest 305-d MY, FCM and ECM, but the Logistic model provided the lowest 305-d MY, FCM and ECM. According to the persistency measures obtained from different models (Table 5), the Schumacher and Morgan equations provided the most persistent lactation curves of MY based on P_2:1, but the Gompertz model provided lactations with better persistency than other equations according to P_3:1. Also, the Morgan model provided the most persistent lactation curves for FCM and ECM based on P_2:1 and Brody and Von Bertalanffy models provided lactations with better persistency for FCM than other equations according to P_3:1. More persistent lactations for ECM were obtained by the Brody, Von Bertalanffy and Schumacher models based on P_3:1. In addition, the Brody and Von Bertalanffy functions provided lactations for MY, FCM and ECM with better persistency than other equations according to P_Weller.

Table 5. Different measures of persistency and 305-d yields for average standard lactations of dairy cows predicted by Brody, Logistic, Gompertz, Schumacher, Von Bertalanffy and Morgan models.

Trait	Variable	Model
		Brody	Logistic	Gompertz	Schumacher	Von Bertalanffy	Morgan
MY	305TY (kg)	9399	8692	10,342	8838	9404	9358
	P2:1	1.02	1.03	1.07	1.00	1.02	1.00
	P3:1	1.04	0.95	0.99	0.89	1.03	0.91
	PWeller (%)	100	90	91	82	100	84
FCM	305TY (kg)	8418	7959	10,113	8878	8418	8348
	P2:1	1.01	1.01	1.07	1.02	1.01	1.00
	P3:1	1.02	0.95	1.05	0.96	1.02	0.94
	PWeller (%)	100	91	99	90	100	89
ECM	305TY (kg)	8331	7818	9698	9057	8329	8248
	P2:1	1.01	1.02	1.07	1.04	1.01	1.00
	P3:1	1.02	0.97	1.05	0.98	1.02	0.96
	PWeller (%)	100	93	99	93	100	91

Notes: MY, unadjusted milk yield; FCM, 4% fat-corrected milk yield; ECM, energy-corrected milk yield; 305TY, predicted 305-d yield for MY, FCM and ECM; P_2:1, Ratio between the milk yields of the second 100 days of lactation and those of the first 100 days; P_3:1,  ratio between the milk yields of the third 100 days of lactation and those of the first 100 days; P_Weller , milk yield persistency measure proposed by Weller et al. (2006).

In general, there was a positive relationship between different persistency measures and 305TY, predicted by different growth models, in this study. Persistency of lactation is associated with yields, especially total yields, but this relationship depends on the type of measures used. The ratio persistency measures (such as the measures in the current study) indicate a positive one, whereas the variation persistency measures show a negative relationship. The reason for this could be that the first are greatly influenced by the level of production and the second are affected by variation in production, with this variation more important for high producing dairy cows (Gengler 1996). In addition, positive relationship between PY and 305TY suggested that cows with high PY can produce more total MY over the lactation than lower producing cow at peak day. Therefore, it was suggested that dairy cows could be selected according to their PYs (Ghavi Hossein-Zadeh 2014).

Decreasing the cost of the dairy production system is possible through improving the lactation persistency, because it is associated with revenue from 305-d milk production, reproductive performance, disease resistance and health and nutrition costs (Dekkers et al. 1996, 1998). More persistent lactations are positively associated with flat lactation curves in dairy cows and the occurrence of reproductive and metabolic problems that originate from the physiological stress of high milk production would be lower in cows with flatter lactation curves (Tekerli et al. 2000; Ghavi Hossein-Zadeh 2014). In this situation the ratio of roughage in the diet could be increased, thus decreasing the costs of production (Tekerli et al. 2000). Therefore, a genetic modification towards a persistent lactation curve could be applied as a means to decline the disease susceptibility in dairy cows. Because of the favourable association between persistency and reduced negative energy balance, it is likely that persistent cows might lose less body weights.

Although different standard models were compared for describing lactation curve of MY in previous studies, the novel items in this study included the ability to use non-conventional models, as growth models, to model lactation curves using a very large data set of Holstein cows and the interest of modelling lactation curve for different measures of MY (FCM and ECM) with various growth functions. Also, different measures of persistency, PT, PY and 305-d yields of MY, FCM and ECM were predicted by different growth models in the current study.

4. Conclusion

Based on criteria to measure goodness of fit, the results of this study showed that selected empirical growth functions were able to fit monthly MY, FCM and ECM records. Of the six growth models investigated, the Morgan model provided the best fit of the lactation curve for MY, FCM and ECM, due to the lower values of AIC than other models, but the Brody model fitted the production data less well than the other equations. In general, there was a positive relationship between different persistency measures and predicted 305-d yield and between PY and 305-d yield for MY, FCM and ECM.

Acknowledgement

Author wishes to acknowledge the Animal Breeding Center of Iran for providing the data used in this study.

ORCID

Navid Ghavi Hossein-Zadeh http://orcid.org/0000-0001-9458-5860