Lactation curves of different cattle breeds under cooperative dairying conditions in Bangladesh

Abstract

This study was undertaken to compare and evaluate the parameters of 10 different mathematical models for their predictive ability in describing the lactation curves for daily milk yield of different genotypes under cooperative dairying in Bangladesh. A database consisting of 7340 herd-test records from 738 cows over the period 1999–2001 was assembled. Breed combinations included Pabna cattle, Australian-Friesian-Sahiwal×Pabna, Holstein×Pabna, Jersey×Pabna and Sahiwal×Pabna. The estimated parameters of the mathematical models, and the predicted lactation milk yields differed significantly between genotypes. The models were evaluated by four fit statistics: Akaike information criteria, coefficient of determination (R ²), root mean square prediction error (RMSPE) and concordance correlation coefficient (CCC). The AIC values indicated that the Ali model provided a good fit for all genotypes. The R ² value suggested that the Legendre polynomial and Ali models were the best fit for all genotypes. CCC and RMSPE values indicated that the best models for all genotypes were Nelder and Wood. Since the CCC value was considered the most informative of the four fit statistics, the Nelder model was the best model to predict the full lactation profile based on test-day records for all genotypes. Using the Nelder model, the predicted 270 day milk yield of the Australian-Friesian-Sahiwal×Pabna genotype was higher (1823 kg) than the other genotypes (1509, 1650, 1531 and 1627 kg for Pabna, Australian-Friesian–Sahiwhal×Pabna, Jersey×Pabna and Sahiwal×Pabna, respectively).

Keywords:

• lactation curves
• models
• genotypes
• predicted milk yield
• model fitness

Introduction

The members of the primary dairy cooperatives of the Bangladesh Milk Producer's Cooperative Union Limited (BMPCUL) milk mostly Pabna cattle, which were developed through the natural mating of local cows with Hariana, Red Sindhi and Sahiwal bulls. The Pabna cattle are regularly crossed with Sahiwal, Holstein-Friesian and Jersey breeds.

In Bangladesh, there is no well-developed recording system for dairy cattle genetic evaluation and improvement. However, non-governmental organisations and private farms maintain milk and health records. The Bangladesh Agricultural University, Bangladesh Livestock Research Institute (BLRI) and BMPCUL also maintain milk, pedigree and health records of their herds. But in the BMPCUL milk-shed area cows’ management and rearing conditions vary throughout the year. During the wet season, all dairy cattle are kept in unplanned semi-intensive housing, where they are fed mainly paddy straw. In the dry season, animals are grazed on pasture land, popularly known as Bathan. Throughout the year, lactating cows are fed 2–3 kg/day/cow of concentrates. For various nutritional experiments, genetic evaluations and other research projects, BMPCUL and BLRI sometimes keep individual cow productive and reproductive performance records of the different genotypes.

Mathematical models are required to predict the total amount of milk yield throughout the entire milking period using regularly recorded daily yields. It is also normal to standardise the number of days in the milking period for all cows. Mathematical models of lactation yields were developed in the 1960s when Wood (Wood 1967) proposed the incomplete gamma function. Before Wood, other workers (Gaines 1927; Nelder 1966), predicted lactation yield with simple mathematical models, which had limited accuracy. Since Wood, various models describing the lactation curve in dairy cows have been reported in the literature, but their applicability under tropical conditions has received limited investigation (Val-Arreola et al. 2004). In Bangladesh, no such studies have been undertaken.

The objectives of this study were to: (1) examine the lactation curves for Pabna and crossbred cows, (2) compare a number of mathematical models for simulating lactation curves and (3) determine the most suitable mathematical model(s), for predicting milk yield under cooperative dairying in Bangladesh.

Materials and methods

Data sources

Measurements of milk production traits of Pabna cattle, Australian-Friesian-Sahiwal×Pabna, Holstein×Pabna, Jersey×Pabna and Sahiwal×Pabna were obtained using records from the Animal Breeding Section of BMPCUL. The data covered the calving years from 1999 to 2001, with each year being divided into two calving seasons: dry (November to May) and wet (June to October). Milk yield data were recorded within 15 days of calving, and collected continuously until the end of lactation at an interval of 30 days for each individual cow. Occasionally, due to unavoidable circumstances, the milk recorder failed to record milk data on the desired day resulting in a small number of cows having records not taken exactly on the 30-day time schedule. Cows with less than four herd test records were deleted from the data, leaving 7340 records for analysis. In this study, nine test-day milk yield records on 738 cows were used to predict lactation milk yields (PLMY) over 270 days. In addition, lactation milk yields were calculated by BMPCUL staff as 15 days times the 15 day milk yield plus 30 days times 45 day milk yield etc. up to the end of lactation. This calculation takes no account of the shape of the lactation curve, nevertheless it will give a modestly accurate estimate of the total yield due to the regular recording of the milk yield. The averages of lactation milk yield and lactation length of all the genotypes were estimated and differences between means were performed using the least significant difference test (Steel et al. 1997).

Lactation equations

Ten mathematical models were used to model the lactation curve of daily milk yields recorded at 30 day intervals (Table 1). For the equations (Table 1), parameter a is always the y-intercept; it controls the vertical position of the curves when plotting daily milk yield against days in lactation. Generally, parameters b and c control the height of the maximum yield and parameters d and k shift the curve up or down the y-axis. The parameters a, b, c, d and k, together with days in milk define the shape and position of each curve.

Table 1. Mathematical models used to describe the lactation curve of dairy cows.

Model	Equation^a	Description
Gaines (Gaines 1927)	Y = ae ^−bt	A two-parameter model that takes no account of a rise to peak yield.
Wood (Wood 1967)	Y = at ^b e ^−ct	The model takes into account the rise to peak lactation, and the decay to drying off.
Polynomial	Y = a + bt + ct ^2	Polynomial can provide a very good fit to data. It is unbound; that is, as x (days in milk) increases indefinitely, the function, y increases or decreases without limit.
Sikka (Sikka 1950)		The model considers the initial rise in the milk yields.
Ali (Ali and Schaeffer 1987)	Y=a+b(t/270) + c(t/270)^2+dlog(270/t)+ klog(270/t)^2	Both random regressions and covariance functions were selected to compare the differences in model flexibility due to the number of parameters and the degree of correlation between them.
Legendre Polynomial (Brotherstone et al. 2000)	Y=aP 0 + bP 1 + cP 2	Functions of P 's of Legendre polynomials were calculated as P 0(t) = 1, P 1(t) = x and P 2(t) =1/2(3x ^2−1) where x=(2(t– (270 + 1)))/(270–1), and t is the ith day in milk.
Wilmink (Wilmink 1987)	Y=a+bt+ce ^−kt	This model combines exponential and linear functions. The parameters of this model can be reduced to a three-parameter linear model by assuming the k exponent is a suitable fixed value.
Nelder^b (Nelder 1966)	Y=t/[a+bt+ct ^2+dt ^3]	This model is a derivation of the Sikka model and an inverse exponential parabolic function.
Rook (Rook et al. 1993)	Y=a{1/[1 + b/(c+t)]}e ^−dt	The model represents the lactation curve as a multiplicative mixture of cell growth and death processes.
Dijkstra (Dijkstra et al. 1997)	Y=ae[b(1 – e ^−ct)/c – dt]	The model is based on a set of differential equations representing cell proliferation and cell death in mammary gland.
^a Y is milk yield (kg/d), t is days in milk, and a, b, c, d and k are parameters that define the shape of the lactation curve. It should be noted that the same symbol in different equations may represent a different component of lactation curves.
^bIn the original paper of Nelder (Nelder 1966), the model was presented without a cubic term (parameters a, b and c); however in this study with the addition of parameter d, a modified Nelder model was used.

Statistical analysis

The Ganies and Wood model functions were fitted using PROC NLMIXED of SAS (SAS 2000). The fifth-order polynomial, Legendre polynomial and Ali polynomial regression, were fitted using PROC MIXED of SAS to estimate their parameters. All other models were fitted using PROC NLIN of SAS.

The goodness-of-fit of predicted values to actual test-day records was tested using four model statistics: the Akaike information criteria (AIC); coefficient of determination (R ²); concordance correlation coefficient (CCC) and root mean square prediction error (RMSPE). The AIC combines the maximum likelihood (data-fitting) and the number of parameters of model, by penalising the (log) maximum likelihood with a term related to model complexity (Val-Arreola et al. 2004). Smaller values of AIC indicate a better fit when comparing the models.

The CCC is a comprehensive statistic used to test simultaneously and jointly for accuracy and precision of the models (Runze and Chow 2005) and assesses the systematic error of actual and predicted values (Pierre 2005), with high values indicating a superior fit. The CCC (Lin 1989) was calculated as:

A _i is the actual milk yield at day i and P _i is the predicted milk yield at day i.

The following statistical model was used to obtain the least square means (±standard error) for each parameter (a, b, c, d and k), estimated from the 10 mathematical models, the predicted 270 day lactation milk yields and for each of the fitness measures.

where Y _ijkl represents the parameter (predicted milk yield, a, b, c, d, k and fit statistics); µ is the overall mean; β_i is the effect of ith genotype; γ_j is the effect of jth year of calving; S _k is the effect of kth calving season; e _ijkl is the random residual effect distributed as N(0, σ²).

The model was fitted considering only the effects of genotype, year and season without interactions as the interaction of the year and calving season was found to be non-significant. To select the best model(s) for each breed(s), the fit statistics values were used.

Results and discussion

The number of cows, total lactation milk yield estimated by BMPCUL staff and lactation length of different genotypes varied (Table 2). There were only small numbers of cows for Australasian-Friesian-Sahiwal×Pabna and Jersey×Pabna genotypes. The calculation of lactation yields by BMPCUL technicians showed that the Australasian-Friesian-Sahiwal×Pabna and Holstein×Pabna genotypes produced higher lactation yields and had longer lactation lengths than other genotypes. Similar findings were reported by several authors (e.g. Bhuiyan et al. 1998; Khan and Khatun 1998; Hossain 2006). However, a study by Talukder et al. (2001) reported lower milk yields for Sahiwal×Local and Jersey×Local genotypes than in the current study, which might be due to the use of Deshi rather than Pabna cattle in their study. Furthermore, many findings (e.g. Cunningham and Syrstad 1987; Syrstad 1989; Madalena et al. 1990) from other tropical countries have also shown that the first cross of temperate breeds with tropical breeds produce more milk than pure breeds in a tropical environment.

Table 2. The number of cows, lactation yield calculated by BMPCUL technicians and lactation length of different breeds, estimated from test-day records.

Breed group	No. of cows	Lactation milk yield, kg (Mean±SE)	Lactation length, days (Mean±SE)
P	185	1503^y±31.15	268.63^y±3.10
A×P	53	1632^x±62.09	282.80^x±5.28
H×P	242	1619^x±28.84	280.35^xy±2.61
J×P	31	1517^y±65.39	267.48^y±2.61
S×P	207	1548^xy±24.89	271.30^y± 6.84
Note: Means with different superscripts within the same column are significantly different at p<0.05.
P, Pabna cattle; A, Australian-Friesian-Sahiwal; H, Holstein; J, Jersey; S, Sahiwal.

Model performance and model parameters

The 10 mathematical models fitting lactation curves to the test-day data were evaluated based on the four different fit statistics (Table 3). The Rook and Diskstra models were not considered further due to their unsuitable model parameters and also poor fit statistics value. The Gaines model was also excluded from further consideration, as it did not allow for the shape of the lactation curve. Four fit statistics (AIC, R ², CCC and RMSPE) values were used to evaluate the performance of the remaining seven models. Smaller AIC and RMSPE values indicated a better fit (Motulsky and Christopoulos 2003; Val-Arreola et al. 2004), but for R ² and CCC higher values indicated superior models (Motulsky and Christopoulos 2003; Pierre 2005). Within each model, the fit statistics were fairly consistent between breeds. Among the four fit statistics, the CCC value was considered the most informative due to its accuracy, precision and the possibility of assessing the systematic error for actual and predicted values (Runze and Chow 2005; Pierre 2005).

Table 3. Comparison of model performance using different fit statistics (±standard error).

		Equations
Breed group	Item	Wood	Legendre polynomial	Polynomial	Sikka	Ali	Nelder	Wilmink
P	AIC	21.47±0.71	43.04^xy±0.63	22.27±0.69	22.43±1.40	11.50±0.49	19.64±0.79	25.94±0.71
	R ^2	0.78^xy±0.01	0.95±0.004	0.82^x±0.01	0.71±0.02	0.89^y±0.011	0.82^x±0.01	0.78^x±0.017
	RMSPE	0.24^xy±0.08	1.30±0.06	0.63^x±0.09	1.2^y±0.52	0.55±0.03	0.20^y±0.05	0.70^y±0.22
	CCC	0.87^x±0.01	0.70^x±0.02	0.87^y±0.01	0.82±0.02	0.87^xy±0.005	0.88^xy±0.03	0.86^y±0.14
A×P	AIC	20.22±1.30	44.96^x±1.17	23.96±1.31	21.69±2.38	10.31±0.85	18.59±1.48	23.92±1.20
	R ^2	0.80^xy±0.03	0.95±0.007	0.82^x±0.02	0.67±0.03	0.93^x±0.018	0.84^x±0.02	0.73^y±0.03
	RMSPE	0.23±0.01	1.45±0.11	0.58^x±0.15	0.79^x±0.68	0.54±0.05	0.17±0.06	0.62±0.30
	CCC	0.89^x±0.02	0.65^y±0.03	0.90^x±0.02	0.77±0.03	0.89^y±0.01	0.90^xy±0.03	0.84^y±0.12
H×P	AIC	20.80±0.58	45.58^x±0.52	23.60±0.58	23.28±1.18	11.54±0.42	19.40±0.66	25.31±0.59
	R ^2	0.84^x±0.01	0.96±0.003	0.85^x±0.01	0.73±0.017	0.92^x±0.009	0.86^x±0.01	0.81^x±0.01
	RMSPE	0.21^x±0.07	1.31±0.05	0.66^y±0.08	0.80^x±0.48	0.54±0.022	0.17^x±0.04	0.62^y±0.20
	CCC	0.88^x±0.01	0.67^xy±0.01	0.86^y±0.02	0.80±0.015	0.88^y±0.003	0.92^x±0.01	0.84^y±0.12
J×P	AIC	19.65±1.74	40.67^y±1.56	21.99±1.76	20.34±3.09	11.27±1.10	20.12±1.95	24.17±1.56
	R ^2	0.73^y±0.03	0.94±0.009	0.78^y±0.03	0.74±0.04	0.85^y±0.023	0.75^y±0.03	0.78^x±0.04
	RMSPE	0.28^y±0.02	1.51±0.18	0.59^x±0.19	0.82^x±0.77	0.55±0.06	0.20^y±0.07	0.54^x±0.33
	CCC	0.85^y±0.009	0.66^y±0.03	0.86^y±0.02	0.78±0.04	0.84^x±0.02	0.87^y±0.03	0.87^x±0.31
S×P	AIC	21.08±0.54	44.62^x±0.48	23.16±0.54	21.56±1.18	10.68±0.42	18.76±0.61	25.50±0.59
	R ^2	0.80^x±0.01	0.96±0.003	0.85^x±0.01	0.72±0.017	0.91^x±0.009	0.84±0.01	0.79^x±0.014
	RMSPE	0.21^x±0.06	1.18±0.05	0.66^y±0.08	0.87^x±0.69	0.53±0.022	0.18^xy±0.04	0.65^y±0.20
	CCC	0.89^x±0.009	0.68^x±0.01	0.86^y±0.01	0.79±0.015	0.88^y±0.004	0.89^xy±0.01	0.85^y±0.12
Note: Means with different superscripts are different at p<0.05 between breeds within model.
A, Australian-Friesian-Sahiwal; H, Holstein; J, Jersey; P, Pabna cattle; S, Sahiwal; AIC, Akaike information criteria, R ^2, RSquare; RMSPE, root mean square prediction error; CCC, concordance correlation coefficient.

The poorest models were: AIC values – Legendre polynomials and Wilmink; R ² values-Sikka and Wilmink; RMSPE values – Legendre polynomials and Sikka and for CCC – Legendre Polynomial and Sikka models. The best models were: AIC values – Ali and Nelder; R ² – Legendre polynomials and Ali; RMSPE – Nelder and Wood and for CCC – Nelder and Wood models (Table 3). Therefore, according to the four fit statistics the best four (Nelder, Ali, Wood and Polynomials) models for five different genotypes were presented for the shape of the lactation curve, model performance and model parameters.

The variation of fit of the models between breeds may have arisen from the differences in breeds, mathematical functions of the models, differences of test-day yield and the amount of data. The differences in the lactation curves are a combination of environmental and genetic factors (Pérochon et al. 1996). The numbers of test-day records and the intervals between tests have effects on the estimation of fit statistics, and also the model fitting ability (Tekerli et al. 2000; Wiggans et al. 2002).

Results of this study indicated that the Nelder model was the most suitable function to model test-day milk yields for all genotypes based on higher CCC and RMSPE and they also performed well enough for the other two criteria. The Ali and Wood models were the next two best models.

Shape of the lactation curves

The raw test day milk yields of the five genotypes and the best fit lactation curve generated by the Nelder model were plotted (Figure 1). The figures showed that the time of peak milk yield calculated from the raw data was about 75 days in milk (DIM) for all genotypes, followed by a gradual decline. The Nelder model gave a peak at 45 DIM for the Pabna, Holstein×Pabna and Sahiwal×Pabna genotypes, 15 DIM for the Australian-Friesian-Sahiwal×Pabna genotype and 73 DIM for the Jersey×Pabna genotype.

Figure 1.  Lactation curves for all the cows of different genotypes obtained from the test-day data and by fitting the Nelder equation.

The gradual decline in the lactation curve has been reported in previous studies (Tekerli et al. 2000; Val-Arreola et al. 2004). The most common shape was a rapid increase, followed by a gradual decline until the cow is dried off. The other shape is a gradual decline from parturition. Several studies (e.g. Macciotta et al. 2005; Fathi Nasri et al. 2008), have shown that the differences in the general shape and the differences of model parameters are associated with yield at the beginning of lactation, the degree of incline and decline before and after peak yield, days in milk at peak and peak and lactation yields of the lactation curves.

The total lactation milk yields calculated from the raw data (Table 2) were similar to the total lactation yields predicted by the Nelder model (Table 4), except for the Holstein×Pabna genotype. The reasons for the similar milk yield between calculated and predicted were due to the regular milk testing by the BMPCUL staff. For the Holstein×Pabna genotype the Nelder model predicted 204 kg of milk more than that from the raw data. The predicted milk yield of Holstein×Pabna was significantly higher due to the breed effects and the predicted ability of the model parameters. Similar factors were reported for the differences in predicted milk yield from lactation curve studies from elsewhere (Pérochon et al. 1996; Vargas et al. 2000; Brown et al. 2001).

Table 4. Parameter estimates (±standard error) and other measures when Nelder model was fitted to the different breed groups in Bangladesh.

Breed group	Parameter	Estimated value
P	a	0.66^y±0.15
	b	0.12^x±0.04
	c	0.00019^y±0.00015
	d	1.465E-6^y±5.003E-6
	Plmy	1509.07^y±58.34
A×P	a	0.82^x±0.15
	b	0.06^y±0.017
	c	0.00061^x±0.00025
	d	−3.1E-7^x±8.6E-6
	Plmy	1649.91^x±98.46
H×P	a	0.55^y±0.11
	b	0.11^x±0.007
	c	0.00004^y±0.0001
	d	1.597E-6^y±3.481E-7
	Plmy	1823.43^x±48.88
J×P	a	0.73^y±0.51
	b	0.12^x±0.03
	c	0.00016^y±0.00047
	d	1.233E-6^y±1.583E-6
	Plmy	1531.01^y±126.84
S×P	a	0.54^y±0.08
	b	0.11^x±0.006
	c	0.00011^y±0.00008
	d	1.685E-6^y±2.896E-7
	Plmy	1626.90^xy±49.01
Note: Means with different superscripts are different at p<0.05 between breeds within model.
H, Holstein; J, Jersey; P, Pabna Cattle; S, Sahiwal, a, b, c and d are parameters that define the scale and shape of the lactation curve. Plmy, predicted lactation milk yield.

Effects of genotype

Differences in the shape of the lactation curve between different genotypes were observed, with the differences in day of peak yield, the steepness and also the flatness of the curves and these are the results of the differences for predicted model parameters of different models. These differences for predicted model parameters were also reported by other investigators (Pérochon et al. 1996; Vargas et al. 2000). In addition, the differences of curve shape could be attributed to few test-day records prior to the peak, and also to a low number of available records. Other studies (Rekik and Gara 2004; Silvestre et al. 2006) identified similar results for the estimation of lactation curve shape.

As estimated by the Nelder model (Table 4), the intercept (a) for Australian-Friesian-Sahiwal×Pabna, Holstein×Pabna and Sahiwal×Pabna, the curve shape (b) for Holstein×Pabna and Jersey×Pabna Sahiwal×Pabna and Pabna genotypes and the placement of the parabola (c and d) for the Australian-Friesian-Sahiwal×Pabna genotypes were significantly higher than other genotypes.

Conclusion

Results of this study indicate that the Nelder model was the most suitable to transform test-day milk yields into a 270-day predicted milk yield for all genotypes, based on higher CCC and lower RMSPE and AIC values. The Wood and Ali models were the next best models. The predicted lactation milk yield differed between genotypes, mathematical functions of the model and also the model parameters. The Australian-Friesian-Sahiwal×Pabna and Holstein×Pabna genotypes fitted better than the other genotypes for all models. There were some limitations in this study, such as the low number of records per genotype representing only three years. Confidence that the correct mathematical model has been chosen would be increased if additional data were to be collected and added to the current analyses.

Acknowledgements

The first author is funded by NZAID: Postgraduate scholarship for PhD studies sponsored by The New Zealand Ministry of Foreign Affairs and Trade. The authors are grateful to the Ministry for their financial support and they would like to thank the Institute of Veterinary, Animal and Biomedical Sciences, Massey University, Palmerston North, New Zealand for their active cooperation and support to pursue this study. The authors are also very much grateful to Mr Shakhawat Hussen, Manager, Animal Breeding Section, Bangladesh Milk producer's Cooperative Union Limited for providing the data.