Comparison of nonlinear functions to describe lactation curves for cumulative milk production in buffalo

Abstract

The aim of this study was to examine the suitability of different growth functions (linear, sinusoidal, Gompertz, Schumacher and Richards) to fit cumulative milk production data from buffalo cows. Cumulative milk production at each day in milk was calculated from two published datasets reporting (i) fortnightly test-day milk yield records of the first lactation of Murrah buffalo that had calved during 1977–2012 and (ii) the first lactation records of Jaffarabadi buffalo collected from history-cum-pedigree registers for each quinquennium between 1991 and 2010. Each function was fitted to the lactation curves using nonlinear regression procedures. The Richards and sinusoidal equations provided the smallest root mean square error values, Akaike's and Bayesian information criteria and, therefore, the best fit for the cumulative lactation curves for milk yield. The Richards equation appeared to provide the most accurate estimate of the cumulative milk production at peak milk yield. Sinusoidal and flexible classical growth functions are appropriate to describe cumulative milk production curves and estimate lactation traits in buffalo.

Domestic buffalo have an important role in agricultural production in Asia and some Mediterranean, Balkan and African countries (Pasha and Hayat, 2012; Zicarelli, 2020). Buffalo are resistant to many diseases, being genetically resistant to some ailments, and are farmed in regions with a hot and humid climate where most cattle cannot be reared (Borriello et al., 2006). According to Zicarelli (2020), the world population of buffalo is approximately 200 million head (97% in Asia; 2% in Africa (mainly Egypt), 0.7% in South America and 0.2% in Australia and Europe, mainly Italy). The countries with the largest numbers of dairy buffalo are India, Pakistan, China, Egypt and Nepal (Pasha and Hayat, 2012). In global terms, in 2018 buffalo represented 5.2% of the total number of dairy ruminants and buffalo milk accounted for about 16% of total milk production (Zicarelli, 2020). Overall, buffalo and dairy cattle share several similarities but they have significant physiological differences that classify them as distinct species and therefore merit separate investigation. Buffalo milk has a higher fat, protein, lactose, vitamin and mineral content, and hence high nutritional value having twice as many calories as milk from dairy cattle (Abd El-Salam and El-Shibiny, 2011; Zicarelli, 2020). Buffalo carcasses, because of their larger size, have harder bones and their meat has a lower fat content and lower muscle pH than beef (Grossi et al., 2013; Di Stasio and Brugiapaglia, 2021).

Algebraic modelling of lactation curves offers a summary of milk yield patterns over lactation time, from which relevant lactation traits (e.g., time and yield at peak) and genetic parameters (Sarubbi et al., 2012) can be estimated or by which total milk yield may be predicted from incomplete data. Appropriate models provide useful information for breeding and management decisions at both the industry and farm levels and for comparison of alternative production strategies in bio-economic modelling. To ensure accurate decisions pertinent to individual animals or herds, it is important that milk yield is predicted at different stages of lactation with minimum error and from relatively few test dates, reducing the cost and inconvenience of milk recording. From a bio-economic viewpoint, a model for the lactation curve must accurately depict what is expected at the farm level so that a cost can be associated with each particular cow (Quinn et al., 2005).

In a typical buffalo lactation curve, milk production increases from calving to a peak then follows a gradual decline (whose slope is related to persistency) until the cows dry off (Şahin et al., 2015). The description of the lactation curve assists the breeder in predicting total milk yield from part lactations and in taking early decisions regarding culling or selection and the introduction of new genetic material (Chaudhry et al., 2000; Barbosa et al., 2007; Sarubbi et al., 2012). Thus, equations which describe milk production over time can be very useful in breeding programs, herd feeding management, decision-taking on the culling of animals and simulation of production systems. The lactation curve is also important to characterise the performance of the buffalo cow throughout lactation so that parameters such as peak milk yield or total lactation milk yield can be estimated (da Cunha et al., 2010). Unlike cattle and sheep, information on the shape of lactation curves in dairy buffalo is very limited. Modelling lactation data can become intricate when curves show an atypical trajectory or disruptions caused by feeding or health disorders (Macciotta et al., 2005; Şahin et al., 2015). These handicaps can be surmounted by using cumulative lactation curves (López et al., 2015). Currently, the suitability of sinusoidal functions, as newly applied equations, is confirmed for describing growth curves in broiler chicks, quail, dairy heifers and turkeys (Darmani Kuhi et al., 2018, 2019a, 2019b, 2019c). Therefore, the objective of this study was to investigate the suitability of sinusoidal and conventional growth functions (viz. linear, Gompertz, Schumacher and Richards) to describe the cumulative lactation curve in dairy buffalo.

Materials and methods

Data sources

Two published datasets of test-day records of milk yield of Murrah and Jaffarabadi buffalo were used to obtain cumulative lactation curves. The data used were those openly reported in each publication. One dataset was derived from 18871 fortnightly test-day milk yield records monitored during the first lactation of 961 Murrah buffalo that had calved between 1977 and 2012 (Sahoo et al., 2014). The data reported in the publication were the average milk yields (means of all the values recorded at a given time) recorded from 5 to 305 d in milk (DIM) with intervals of 15 d, resulting in lactation curves with a total of 21 datapoints. The other dataset was developed from first lactation records of 213 Jaffarabadi buffalo from history-cum-pedigree catalogues gathering data over a period of 20 years (from 1991 to 2010: Sharma et al., 2017). In this publication, four lactation curves were reported, corresponding to the average curves for each of the four quinquennia (1991–1995, 1996–2000, 2001–2005, 2006–2010) encompassed in the study. Each curve comprised 11 datapoints corresponding to the mean milk yields observed at day 5 after calving, and then every 30 d up to the end of lactation at 305 DIM. A detailed description of the data can be found in the original publications; some descriptive lactation curve traits of the datasets used in the study are shown in the supplementary material (Table S1). Average daily milk yield and cumulative production at each DIM were calculated following the guidelines of the International Committee for Animal Recording for buffalo milk recording (International Committee for Animal Recording – The Global Standard for Livestock Data, 2017).

Growth models

The equations used to describe the lactation curves are presented in Table S2 of the supplementary material. The linear, Gompertz, Schumacher, Richards and sinusoidal functions (the latter four represent sigmoidal curves with a fixed or flexible point of inflexion) were fitted to the data to model the relationship between cumulative milk production and DIM.

Statistical procedures

Each function was fitted using the SigmaPlot software version 12.0 (Systat Software Inc., 2012) to estimate model parameters. To fit the nonlinear functions, the Gauss–Newton algorithm was used as the iteration method. Once the fitting converged, the goodness-of-fit (quality of prediction) was assessed statistically by using the adjusted coefficient of determination ($R_{{\rm adj}}^2 $), residual standard deviation or root mean square error (RMSE), Akaike's information criterion (AIC) and the Bayesian information criterion (BIC), calculated as illustrated in the supplementary material (Table S3). The fit of a given model is more satisfactory if $R_{{\rm adj}}^2 $ is closer to unity and a better goodness-of-fit corresponds to smaller numerical values of RMSE, AIC or BIC.

Results

Estimated parameters of the equations used to fit the cumulative lactation curves of buffalo are shown in Supplementary Table S4. The Gompertz equation over-estimated initial milk yield in all datasets (estimated values ranging from 58 to 112 kg whereas observed values at 5 DIM were 2.5–4 kg/d, see Fig. 1 as an example). The linear function also over-estimated yield at calving (102 kg) in one dataset, and provided inconsistent estimates with other datasets (negative initial yields). In contrast, the Richards, Schumacher and sinusoidal values were close to the expected initial average milk yield (Fig. 1 and Fig. S1 of the supplementary material). There were significant differences among the different functions for final (asymptotic) cumulative milk production. The values estimated with the Schumacher and Richards were considerably higher than those obtained with the Gompertz and sinusoidal equations. Estimates of asymptotic cumulative milk production appeared to over-estimate total lactation yield, in particular with the Schumacher and Richards equations.

Figure 1. Plots of cumulative milk production (kg) against days in milk showing the fit of different functions to the data of Sahoo et al. (2014) (dots are the observed values and the solid line is the fitted curve for each function).

Statistical goodness of fit for the five functions fitted to the cumulative milk production curves is presented in Table 1. As can be seen, the $R_{{\rm adj}}^2 $ values (proportion of variance accounted for by the model) were always greater than 0.98, indicating suitable fits to the observed data, and showed little differences among the functions. Comparing the functions using RMSE, AIC and BIC (Table 1) showed differences in the goodness-of-fit. The Richards and sinusoidal equations gave the lowest values of RMSE, AIC and BIC, and therefore the best fit to the curves for cumulative milk production. The worst fit was for the linear function, showing the largest values of RMSE, AIC and BIC.

Table 1. Comparison between the general goodness-of-fit of the linear and nonlinear equations for different data profiles based on statistical criteria^{a, b}

^a $R_{{\rm adj}}^2 $, Adjusted coefficient of determination; RMSE, Root Mean Square Error; AIC, Akaike Information Criterion; BIC, Bayesian Information Criterion.

^b Best model is indicated in bold font.

Both cumulative production and DIM at peak milk yield for the different nonlinear functions fitted are shown in Table 2. Based on the comparison between the observed datapoints and the values predicted by each function, the Richards equation appeared to provide more accurate estimates of peak milk yield than the other functions. Limitations on space prevent presentation of the full results for all datasets. As examples of model behaviour in predicting the cumulative lactation curves for the five functions (supplementary material), observed datapoints and fitted curves with the data of Sahoo et al. (2014) are illustrated in Fig. 1, whereas Supplementary Figure S1 shows the results for the dataset reported by Sharma et al. (2017) during the years 2006–2010.

Table 2. Time (t*), cumulative production (y*) and daily milk yield (m_p) at peak and milk production at 305 days lactation (y ₃₀₅) for the models fitted to cumulative milk production curves in buffalo

Discussion

Modelling seeks to find parametric descriptors of the shape of the lactation curve to predict features such as peak milk yield, time to peak yield, and total milk production over lactation. Knowing when peak milk yield will occur can assist dairy farmers and managers in planning feeding strategies to maintain peak yield for as long as possible (López et al., 2015). Knowledge of the lactation curve shape is also a valuable tool in the management context, especially for time-dependent decisions. Although several studies have compared nonlinear models to fit lactation curves in dairy cows (Ghavi Hossein-Zadeh, 2014, 2017; López et al., 2015) and buffalo (Şahin et al., 2015; Ghavi Hossein-Zadeh, 2016), there are no reports on modelling the cumulative lactation curve in dairy buffalo. The shape and trajectory of lactation curves are similar in dairy cattle and buffalo, but there are differences in magnitude and rates of change that make it opportune to investigate the modelling of buffalo lactation curves. The current study is the first to report on fitting cumulative milk production data from buffalo, particularly assessing a sinusoidal equation as an alternative to conventional functions. Developing an appropriate strategy to reach the desired lactation shape through modifying the parameters of such nonlinear models can be of much value. These parameters may act as surrogates for various genetic and environmental factors, different combinations of which can account for differences in nonlinear curves.

Typically, the conventional lactation curve follows a skewed bell shape, with a sharp rise after parturition to a peak a few weeks later followed by a slow decline until the end of lactation. However, several researchers have reported non-standard lactation patterns such as continuously increasing or decreasing or reverse shape of the most typical curve. These atypical curves are not fitted easily using conventional lactation models (Rekik and Gara, 2004; Macciotta et al., 2005; Silvestre et al., 2005; López et al., 2015; Şahin et al., 2015). Under field conditions, unexpected incidents may cause disturbances in the course of lactation resulting in irregular curves that are more difficult to be fitted by the conventional functions. The monotonically increasing pattern of cumulative curves attenuates the possible impact of any disturbance on the curve and does not present major difficulties in fitting sigmoidal functions to cumulative milk production (López et al., 2015). The current results indicate that representative nonlinear growth functions, especially the Richards and the sinusoidal equations, provide suitable fits to cumulative lactation curves in dairy buffalo. Consistent with the results of this study, López et al. (2015) fitted six classical growth functions (monomolecular, Schumacher, Gompertz, logistic, Richards, and Morgan) to cumulative milk production curves of Canadian Holstein dairy cows. They concluded that the fitted classical growth functions could be an alternative to conventional models when analysing lactation data. Ghavi Hossein-Zadeh et al. (2020) investigated the appropriateness of a sinusoidal function by applying it to model the cumulative lactation curves for milk production and composition in primiparous Holstein cows in Iran and comparing it with three conventional growth models (linear, Richards and Morgan). Classical growth functions could be fitted accurately to cumulative lactation curves for production traits, but the new sinusoidal equation provided a better goodness of fit for some of the curves evaluated (Ghavi Hossein-Zadeh et al., 2020).

A cumulative lactation curve shows a monotonically increasing trend presenting a typical sigmoidal shape, with an accelerating phase related to early lactation when daily production rises, and a declining phase in which the slope decreases continuously (Abdel-Salam et al., 2011; Borghese et al., 2013). Yield over the course of lactation is simply obtained from this cumulative milk production curve by calculating its value at any specific DIM for each nonlinear function (López et al., 2015). This simple method of obtaining total production with cumulative curves differs from the method of calculating by integrating standard lactation curves to obtain total milk yield. The current study results show that monthly test-day production records provide a suitable input to the cumulative lactation curve. This will be a benefit if a daily system of recording production data in dairy farms is not available.

Initial milk yield at calving (y ₀) was over-estimated by the Gompertz and linear equations. Figure 1 and Fig. S1 show that the linear function systematically over- or under-estimates the observed values at certain sections of the curve (based on the small number of runs of sign) and thus it seems inappropriate for fitting cumulative lactation curves. The Gompertz equation is less flexible than the other sigmoidal functions, with a fixed inflexion point relative to the upper asymptote (see Table S2). To approach this fixed inflexion point, the fitted curve starts at a relative high value at DIM = 0 (thus over-estimating the observed value). The slow rate of increase (slope) at the initial stage of the curve counterweighs the initial high y ₀ estimate decreasing its impact on other parameters. As a result, the estimate of peak milk yield is similar to that obtained with other flexible sigmoidal equations (Richards or sinusoidal) providing lower y ₀ estimates but showing faster rates of increase at the initial DIM.

Comparison with the observed values for peak milk yield (calculated from the gradient between two consecutive datapoints and shown in Supplementary Table S1 as m_y ⁺) confirmed that all the sigmoidal functions provided accurate approximations to this key parameter, closely related to total milk production for the whole lactation. The Pearson correlation coefficients between observed and estimated values ranged between 0.988 for the Gompertz and 0.993 for the Richards (for the linear function it was 0.969). There were, however, substantial differences among functions in the estimates of DIM at peak milk yield, which was always at a later time with the Gompertz and at an earlier time with the Richards. Time at peak milk yield will determine total milk production before and after peak, and an accurate estimate is required if parameters such as persistency are to be assessed. However, to evaluate which model provides a better approximation of DIM at peak milk yield, lactation curves with daily records would be necessary to assess the time at which peak yield occurs. Estimates of DIM and milk yield at peak and of total milk production with the Richards, Schumacher and sinusoidal functions are comparable to those reported elsewhere (Metry et al., 1994; Şahin et al., 2015; de Carvalho et al., 2020) using conventional lactation models (e.g., Wood's equation) for different breeds of buffalo. Compared with conventional lactation models, which are sensitive to missing test-day records, especially at early DIM (Adediran et al., 2007; Wasike et al., 2011), fitting nonlinear growth functions to cumulative production records would be influenced to a much lesser degree by a few missing records (López et al., 2015). Therefore, the smoother trend of cumulative curves for milk yield is likely less susceptible to abnormal records or the occurrence of outliers. Abnormal records can be a result of several problems such as measurement faults and biological disruptions in animal performance occurring when cows are threatened by some severe conditions that limit the expression of their genetic potential, e.g., nutritional deficit, metabolic or infectious disorders (Codrea et al., 2011). The identification of these bias factors and correcting them is critical and this may be achieved by analysing daily records of production traits. However, as mentioned above, these disruptions in the lactation curve would be very problematic when working with conventional lactation models, whereas fitting growth functions to cumulative milk production and composition demonstrates the greater flexibility of such models to varied data attributes (López et al., 2015).

In conclusion, this study shows that sigmoidal growth functions can be fitted accurately to the cumulative lactation curves of dairy buffalo, resulting in acceptable statistical goodness-of-fit and accurate determination of production traits (DIM and daily milk yield at peak, total milk production). This information is useful for designing genetic programs or for management or feeding strategy planning. Although all the functions tested seem to be suitable to represent buffalo lactation curves, the results show enough evidence to assume that the Richards resulted in better accuracy of estimation of lactation traits. The sinusoidal equation had sufficient flexibility (variable inflexion point) and provided an accurate fit to the data, so it may be considered an alternative to classical growth functions.