Domestic buffalo have an important role in agricultural production in Asia and some Mediterranean, Balkan and African countries (Pasha and Hayat, Reference Pasha and Hayat2012; Zicarelli, Reference Zicarelli2020). 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., Reference Borriello, Capparelli, Bianco, Fenizia, Alfano, Capuano, Ercolini, Parisi, Roperto and Iannelli2006). According to Zicarelli (Reference Zicarelli2020), 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, Reference Pasha and Hayat2012). 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, Reference Zicarelli2020). 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, Reference Abd El-Salam and El-Shibiny2011; Zicarelli, Reference Zicarelli2020). 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., Reference Grossi, Addi, Borriello and Musco2013; Di Stasio and Brugiapaglia, Reference di Stasio and Brugiapaglia2021).
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., Reference Sarubbi, Auriemma, Baculo and Palomba2012) 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., Reference Quinn, Killen and Buckley2005).
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., Reference Şahin, Ulutaş, Arda, Yüksel and Serdar2015). 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., Reference Chaudhry, Khan, Mohiuddin and Mustafa2000; Barbosa et al., Reference Barbosa, Pereira, Santoro, Batista and Neto2007; Sarubbi et al., Reference Sarubbi, Auriemma, Baculo and Palomba2012). 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., Reference da Cunha, Pereira, e Silva, de Campos, Braga and Martuscello2010). 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., Reference Macciotta, Vicario and Cappio-Borlino2005; Şahin et al., Reference Şahin, Ulutaş, Arda, Yüksel and Serdar2015). These handicaps can be surmounted by using cumulative lactation curves (López et al., Reference López, France, Odongo, McBride, Kebreab, AlZahal, McBride and Dijkstra2015). 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., Reference Darmani Kuhi, Shabanpour, Mohit, Falahi and France2018, Reference Darmani Kuhi, France, López and Ghavi Hossein-Zadeh2019a, Reference Darmani Kuhi, Ghavi Hossein-Zadeh, López, Falahi and France2019b, Reference Darmani Kuhi, López, France, Mohit, Shabanpour, Ghavi Hossein-Zadeh and Falahi2019c). 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., Reference Sahoo, Sahoo, Singh, Shivahre, Singh, Dash and Dash2014). 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., Reference Sharma, Gajbhiye, Ramani, Ahlawat and Dongre2017). 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. (Reference Sahoo, Sahoo, Singh, Shivahre, Singh, Dash and Dash2014) (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 criteriaa, 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. (Reference Sahoo, Sahoo, Singh, Shivahre, Singh, Dash and Dash2014) are illustrated in Fig. 1, whereas Supplementary Figure S1 shows the results for the dataset reported by Sharma et al. (Reference Sharma, Gajbhiye, Ramani, Ahlawat and Dongre2017) during the years 2006–2010.
Table 2. Time (t*), cumulative production (y*) and daily milk yield (mp) at peak and milk production at 305 days lactation (y 305) 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., Reference López, France, Odongo, McBride, Kebreab, AlZahal, McBride and Dijkstra2015). 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, Reference Ghavi Hossein-Zadeh2014, Reference Ghavi Hossein-Zadeh2017; López et al., Reference López, France, Odongo, McBride, Kebreab, AlZahal, McBride and Dijkstra2015) and buffalo (Şahin et al., Reference Şahin, Ulutaş, Arda, Yüksel and Serdar2015; Ghavi Hossein-Zadeh, Reference Ghavi Hossein-Zadeh2016), 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, Reference Rekik and Gara2004; Macciotta et al., Reference Macciotta, Vicario and Cappio-Borlino2005; Silvestre et al., Reference Silvestre, Petim-Batista and Colaço2005; López et al., Reference López, France, Odongo, McBride, Kebreab, AlZahal, McBride and Dijkstra2015; Şahin et al., Reference Şahin, Ulutaş, Arda, Yüksel and Serdar2015). 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., Reference López, France, Odongo, McBride, Kebreab, AlZahal, McBride and Dijkstra2015). 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. (Reference López, France, Odongo, McBride, Kebreab, AlZahal, McBride and Dijkstra2015) 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. (Reference Ghavi Hossein-Zadeh, Darmani Kuhi, France and López2020) 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., Reference Ghavi Hossein-Zadeh, Darmani Kuhi, France and López2020).
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., Reference Abdel-Salam, Mekkawy, Hafez, Zaki and Abou-Bakr2011; Borghese et al., Reference Borghese, Boselli and Rosati2013). 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., Reference López, France, Odongo, McBride, Kebreab, AlZahal, McBride and Dijkstra2015). 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 0) 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 0 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 0 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 my +) 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., Reference Metry, Mourad, Wilk and McDaniel1994; Şahin et al., Reference Şahin, Ulutaş, Arda, Yüksel and Serdar2015; de Carvalho et al., Reference de Carvalho, de Moura, Vieira, Hurtado-Lugo, Aspilcueta-Borquis, Gomide, Tonhati, de Oliveira and de Araujo-Neto2020) 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., Reference Adediran, Malau-Aduli, Roche and Donaghy2007; Wasike et al., Reference Wasike, Kahi and Peters2011), 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., Reference López, France, Odongo, McBride, Kebreab, AlZahal, McBride and Dijkstra2015). 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., Reference Codrea, Højsgaard and Friggens2011). 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., Reference López, France, Odongo, McBride, Kebreab, AlZahal, McBride and Dijkstra2015).
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.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0022029924000062
