The geometry of the lactation curve based on Wood’s equation: a two-step prediction

Lactation records from cows of the southwestern Paraná state, Brazil, form the dataset of this study. We applied the information-theoretic approach to evaluate the ability of the nonlinear Wood, Brody, Dijkstra, and Gamma functions to fit to these data by employing a two-step technique based on nonlinear mixed-effects models and generalized linear mixed-effects models. Wood’s equation was fitted with the combination of a first-order autoregressive correlation structure and a variance function to account for heteroscedasticity. This version was the best choice to mimic lactation records. Some geometric attributes of Wood’s model were deduced, mainly the ascending specific rate from parturition to peak milk yield and the descending specific rate as a measure of the lactation persistence of the milk yield at peak production. Breed and parity order of the cows were assumed as fixed effects to obtain a reliable model fitting process. Regardless of breed, first-order parity cows had greater persistency than their older counterparts, and the greater the ascending rate of milk yield from the parturition to the peak, the sharper the decrease in milk yield post-peak; therefore, the rates (absolute values) of ascending and descending phases correlated positively. Nonetheless, the actual estimated values of the descending phase rates are negative. Wood’s equation was flexible enough to mimic either concaveand convex-shaped lactation profiles. The correlations between both peak milk yield and random estimates for β with total milk yield per lactation were positive. However, peak milk yield might not be the only variable used for ranking cows; the total milk yield integrates all information of the lactation profile through the estimated parameters of Wood’s equation.


Introduction
Lactation functions provide an elegant example to study the applicability of mathematical models to explain nonlinear phenomena in animal science. Scientists have proposed several mathematical functions to mimic milk production records over time (Brody et al., 1924;Wood, 1967;Grossman and Koops, 1988;Rook et al., 1993;Dijkstra et al., 1997;Pollott, 2000). Specialized dairy breeds typically The geometry of the lactation curve based on Wood's equation: a two-step prediction Oliveira et al. 4 and descending phases about peak milk yield are those with a constrained parametric space, i.e. 0 < β ≤ 1; whereas those records that show only a descending phase starting from the onset of lactation at parturition are those for which the parametric space is constrained to β ≤ 0 (concave profiles). We remind the reader that equation 1 resumes to the exponential model for β = 0. The cumulative or total milk production (TMμ t ) is the area under the curve used to describe the lactation records. Here, we provide an analytical solution that yields the definite integral as follows: TMμ t = ∫ l u αt β exp(-λt)dt = -αλ -(β+1) (Γ[β + 1, uλ] -Γ[β + 1, lλ]) Eq. (6) Numerically, the upper limit (u) of the integral was set at 305 d, and the lower limit (l) can be a small value greater than zero (10 -6 ) for numerical integration. The incomplete gamma function is The gamma function Γ[β + 1] = ∫ 0 ∞ t β exp(-t)dt is a numerical result from a definite integral, and the lower incomplete gamma function is given by γ[β + 1, λt] = ∫ 0 t v β exp(-λv)dv.
The geometry of the function depicted in Figure 1 can yield a fractional rate constant of decline with reference to peak milk yield. The dashed descending line tangent to the inflection point after the time for peak milk yield is the geometric place defined by equation 7a: in which (t', μ t' ) is the point where equation 7a crosses the tangent line to the point (t p , μ t p ), and (dμ t ⁄dt) t i is the value of the first derivative of equation (1) at t i (Figure 1a). The tangent to (t p , μ t p ) is parallel to the abscissa axis, i.e., μ t' = μ t p , and the abscissa coordinate t' can be isolated in equation 7a and solved for t = t i . Finally, if we divide both sides of equation 7a by μ t p and manipulate signs for κ d because (dμ t ⁄dt) t i < 0, for not changing the trend of equation 7a, we have the persistency trend over time (P t ) as follows: On panel a, we displayed the whole lactation (solid line) and indicated the inflection point past peak milk yield with abscissa t i as the day at the inflection, and the respective expected milk production rate at the ordinate (μ ti ). Panel b is the enlarged part from day 0 to day 72 of the same lactation profile. The coordinates of the points of initial milk production (μ 10-6 ) at day 0 (t = 10 -6 ≈ 0), peak milk yield (μ tp ) at the day of peak (t p ), and where the tangent to the inflection point (descending dashed line) intercepts the abscissa parallel (dotted lines) that passes the ordinate (μ t' = μ tp ) defines the arbitrary beginning of the descending exponential phase. These two points are important to compute the fractional descending rate (κ d ). The double dotted and dashed lines mark t p , t', and t i . The ascending phase is defined by the cord that passes the coordinates (10 -6 , μ 10-6 ) and (t p , μ tp ). Two dotted lines parallel the abscissa and passes through these two points also united by an ascending dashed line. These lines are important to define the angle δ. This angle is important to compute the fractional ascending rate (κ α ). The circles mark the coordinates on both panels. The expected function was built on parameters estimated for fifth parity order cows of the Holstein breed as part of our results. Days in milk The geometry of the lactation curve based on Wood's equation: a two-step prediction Oliveira et al.
Therefore, the specific rate constant of decline after peak milk yield arbitrarily becomes the positive value, as follows: The new parameter κ d can be d -1 or w -1 and dependent on t i (equation 5) as an intrinsic attribute of equation 1. We can relate the fractional rate described by equation 8 to the phenomenon of stromal apoptosis in the mammary gland (Stefanon et al., 2002;Dijkstra et al., 1997Dijkstra et al., , 2010. The relative rate of decline midway between peak milk yield and the arbitrary end point of lactation (t f = 305) is an additional rate computed from the parameters of Wood's equation, namely r d = 2βλ(β -λt f ) -1 -λ (Dijkstra et al., 2010).
Another critical factor of the lactation might be the rate of DNA proliferation in the mammary gland, which would be responsible for the rise in milk production after parturition (Dijkstra et al., 2010). We can geometrically generate a point estimator for an ascending rate of milk production by calculating the average milk production during the ascending phase until peak milk yield, as follows: . By computing the tangent of the δ angle (Figure 1b), as equation 10, we have: Eq. (10) and the following ratio yields a constant, fractional ascendant rate of milk production given by equation 11. κ α = tan δ⁄μ̅ 10 -6 ,t p Eq. (11) This rate can be d -1 or w -1 . The rate κ α can be associated with the process of cell proliferation until peak milk yield (Dijkstra et al., 1997(Dijkstra et al., , 2010. The constraint 0 < β ≤ 1 must hold for a valid κ α ; nonetheless, even though estimates for the fixed parameter β are negative, there are possible occurrences of positive values from the random cow effect over parameter β that may result in typical ascending-peakdescending (convex-shaped) lactation profiles.

Brody's lactation equation
Brody and colleagues proposed a model to mimic the standard course of milk flow of dairy cows, that is, the rise of milk flow after parturition and its steady decline after peak milk yield (Brody et al., 1924). Their model describes a biphasic phenomenon, in which the rise follows a monomolecular change and, simultaneously, an exponential decay follows the peak milk yield. Mathematically, the final equation has the form: The mean μ t has the same meaning as previously described, but θ 1 and θ 2 are scale parameters, κ 1 (d -1 or w -1 ) is the characteristic constant of the decline of milk production after peak milk yield ("after the second month" in those authors' words), and κ 2 (d -1 or w -1 ) is the characteristic diminishing constant of the rising course of lactation. It is possible that equation 12 mimics concave-shaped lactation profiles, but this will depend on final parametric estimates: if θ 1 > 0 and θ 2 ≤ 0, the equation resumes to an exponential decay. All derivations applicable to Wood's equation also apply to the Brody's equation.

Dijkstra's lactation equation
There is an elegant mechanistic equation based on the pool size of DNA in the cell population of the mammary gland during lactation. This model roots on cell proliferation and death and mimics the typical convex-shaped lactation profile as follows: Parameter μ 0 (kg • d -1 or kg • w -1 ) is the theoretical initial milk production at parturition, θ T (d -1 or w -1 ) represents the specific rate of cell proliferation at parturition, κ 3 (d -1 or w -1 ) is a decay parameter, and κ 4 (d -1 or w -1 ) is the specific rate of cell death. Dijkstra et al. (1997) showed several attributes of equation 13 as general nonlinear functions of its parameters.

The Gamma density function used as a lactation equation
The gamma probability density function (PDF) models time-to-event data (Mood et al., 1974;Stroup, 2013). Because Wood (1967) reparametrized the gamma density function as previously mentioned, we believe that it is worthy evaluating this model in its original form because the properties of the gamma PDF may be of use to describe lactation records. We can compute the mean and variance based on this lactation function. The model can have the following form: The mean μ t represents the lactation course, and α and r are scale and shape parameters, respectively. Parameter λ (d -1 or w -1 ) has the same meaning as previously described. The constraints α 0 > 0, r > 0, and λ > 0 must hold. The equivalence between equations 1 and 14 is trivial, because of β = r -1, α = α 0 λ r ⁄Γ[r], and λ is the parameter in common. One can use the same rationale for Wood's model to obtain the nonlinear parametric functions of interest.

Model fitting: first-step prediction
We used the nlme function of R to fit equations 1, 12, 13, and 14 to the lactation records. The gnls and nlme functions belong to the nlme package (Pinheiro et al., 2017). The stochastic version of the models was set as follows: The record y t ijk corresponds to the milk production rate (kg • d -1 ) for the k-th cow of the i-th breed, during the j-th parity order, and recorded at time t. The basic nonlinear function that fits y t ijk is g(Θ m , t ijk ) = μ t ijk , ∀μ t related to equations 1, 12, 13, and 14. Parameter Θ m represents the vector of parameters of the m-th model. However, prior to the fit of the model equations, we had to group lactation records by the random cow and cow × parity order intersecting effects by using the groupedData function of nlme. The resulting formula for the grouped data was y ~ Time|cow/cpo, i.e., y is daily milk yield, time is t, cow and cpo are the random cow and cow × parity order (cpo) intersection. Breed (B), parity order (PO), and their interaction (BPO) were the fixed effects associated with each parameter of μ t for each model evaluated.
We modeled the variance associated with the error term (e t ijk ) according to four conditions: σ t 2 = σ 2 (ω + |g(Θ m , t ijk )| ψ ) 2 , and Eq. (18) cov(e t ijk , e t' ijk │u m The geometry of the lactation curve based on Wood's equation: a two-step prediction Oliveira et al. 7 Parameters represent the residual variance (σ 2 ), the power (ψ) that scales the absolute mean |g(Θ m , t ijk )|, and the assumed positive correlation parameter (0 < ϕ < 1). The model fitted with the variance function in equation 17 challenges the fitted one based on the traditional homoscedastic assumption described by equation 16. Parameter ω is the intercept of the VarConstPower function (equation 18), which allowed modeling the variances when |g(Θ m , t ijk )| = 0 or equal to a baseline or background value, and when the variance scales to the mean for |g(Θ m , t ijk )| > 0 or greater than a baseline value (Pinheiro and Bates, 2000;Rohem Júnior et al., 2020;Vieira et al., 2020). The fit of the model with the first-order autoregressive correlation structure given by equation 19 challenges the traditional assumption of independent time records on the same subject (Vonesh, 2012); this was accomplished by fitting the corCAR1 function of nlme package (Pinheiro and Bates, 2000). We also checked for possible diagonal and symmetric covariance structures for random effects by using the pdDiag (diagonal, D) and pdSymm (unstructured, U) covariance classes available in the R software (Pinheiro and Bates, 2000).
The models were initially fitted with the gnls function with no random effects introduced. Thence, by using the nlme function, we gradually introduced random effects (cow and cow × parity order) in the different fixed parameters contained in g(Θ m , t ijk ) creating several model versions. An example is the introduction of the random effects over parameter λ in Wood's equation, in sequence over parameters λ and β, and so on for all equations studied, which formed the different model versions. In addition, we assessed the quality of fit of the model solutions by using the I-T approach (Buckland et al., 1997;Anderson, 2004, 2014;Burnham et al., 2011a,b;Sober, 2002). This methodology is based on the corrected Akaike Information Criterion as AICc m (Akaike, 1974;Sugiura, 1978;Hurvich and Tsai, 1989;Cavanaugh, 1997), and the derived measures ∆ m (m-th Akaike differences), p m (m-th model probabilities), and ER m (m-th evidence ratios) computed for each m-th feasible model version fitted. We discarded solutions that yielded heavy-tailed 0.95 confidence intervals (0.95CI) for variance components of the pdDiag or pdSymm structures, as well as for parameters of the correlation structure, the power-of-the-mean function, or both (Pinheiro and Bates, 2000).

Generalized linear mixed-effects models: second step prediction
Only one solution satisfied all criteria for model selection (Table 2). Therefore, we obtained the random effects as outputs of the nlme function and computed the nonlinear geometric functions of the parameters. In sequence, we fitted the values for each cow and cpo for all BPO, which generated new random variables to be analyzed by using the GLIMMIX procedure (SAS University Edition, SAS Systems Inc., Cary, NC, USA). We checked the patterns of the Pearson residuals to evaluate the quality of fit of the model and the assumed probability density function. The model definition was as follows: The linear predictor is η ijk and the Greek uppercase letters represent the i-th breed (B i ), the j-th parity order (Ρ j ), and their interaction (ΒΡ ij ). We assumed that variables could be y ijk |c k ∼ Gamma(μ ijk , Φ) or y ijk |c k ∼ Normal(μ ijk , σ 2 ), and the random effect of cow is identically and independently distributed as Greek capital Φ is a scale parameter, and the link functions used were as follows (Stroup, 2013): the identity link or η ijk = μ ijk was used for the Normal distribution, whereas the inverse and log links, namely η ijk = μ ijk -1 or η ijk = log (μ ijk ), were used for the Gamma PDF. Because the Normal distribution and the Gamma distribution produced equal residual patterns, the alternative Gamma PDF was chosen as the best solution. Therefore, the Gamma was used as a generalizing distribution for variables in the domain (0, ∞).

Results
We challenged the traditional assumption of independence and homoscedasticity among errors. Initially, we fitted the nonlinear models corresponding to the general mean with fixed effects of breed (B), parity order (PO), and their interaction (BPO), without the correlation (equation 18) among repeated measures, and with a simple variance function (equation 16). The gradual introduction of fixed The geometry of the lactation curve based on Wood's equation: a two-step prediction Oliveira et al. 9 effects resulted in an improvement of the model quality of fit (Table 2). The addition of a random effect was a necessary improvement, as well as the introduction of a first-order autoregressive correlation (equation 19 or corCAR1 function, parameter ϕ) and the power-of-the-mean variance function (equation 17, parameter ψ). Nonetheless, the introduction of a single random effect to parameter β of equation 1 was the best solution after an exhaustive and careful choice among several possible and feasible model versions, given the data (Tables 2 and 3). Therefore, the best model in the set of all models studied was formed by combining equation 1 fitted with 62 parameters related to the BPO interactions as the fixed effects, the variance function described by equation 17, the corCAR1 function, and a random effect associated with the shape parameter β. We also presented the standard errors of the variance components of the random effects and their respective 0.95CI related to parameter β, namely, the standard deviations related to the cows (σ cow,β ), the intersecting cow × parity order (σ cpo,β ), and the residual error (σ). There is one parameter estimated from equation 17 and another parameter from equation 19, namely ψ and ϕ, respectively (Table 3).
Equation 1 presented convergence in all cases. In contrast, the other models sometimes failed at convergence or yielded non-positive definite Hessians for variance-covariance matrices. Other problems observed were unreliable (heavy-tailed) confidence intervals for variance-covariance parameters, and we discarded those solutions as recommended by Pinheiro and Bates (2000). Equations 12 and 13 yielded poor performances if compared with the best solution for equation 1. Only model versions with fixed effects were feasible for those models.
The nonlinear geometric functions of the fixed parameters (α, β, and λ) of the best model (equation 1) and respective standard errors were estimated for each BPO (Table 3). The estimate obtained for parameter β of BPO 41 was negative. The sign of the estimate of parameter β determines the shape of the lactation function, so that the curve was concave for β ≤ 0 (Figure 2, panel e) or presented a typical convex profile Typical ascending-peak-descending (convex) lactation profiles are presented on panels a, b, c, d, and f. The breeds shown are Holstein (i = 2) and Jersey (i = 3). An example of a concave, non-peak lactation profile is presented on panel e for first parity order crossbred (i = 4) cows. Interestingly, the same breed (crossbred cows) exhibited the two patterns.   The geometry of the lactation curve based on Wood's equation: a two-step prediction Oliveira et al.

11
for 0 < β ≤ 1 (Figure 2, panels a, b, c, d, and f). A group of 385 cows presented a final negative random value for β; therefore, 6369 cows presented typical convex-shaped lactation profiles. The curves fitted with fixed parameter β ≤ 0 did not present peak milk yields as observed for convex-shaped profiles, and we found those results for the first and third parity orders of the crossbreds, namely BPO 41 and BPO 43 . Possibly, the low number of records and cows of the third and fourth parity orders limit the strength of the inferences taken.
Other attributes computed from the two-step prediction, such as μ t p , κ α , κ d , t p , etc., are presented with standard errors of the means and respective 0.95CI for all BPO combinations (Tables 4 and 5). Sixth parity order cows of the Holstein breed produced more milk per lactation than any other parity order for Holsteins or breeds, as one can check by comparing the confidence intervals reported. Therefore, Holsteins produced more milk than any other breed in the Southwest Paraná state. Fourth-, fifth-, and sixth-order Holstein cows presented the same peak milk yield, as demonstrated by overlapping confidence intervals (Table 4). The combination of the two-step statistical tools allowed the identification of these attributes for ranking animals within breeds and groups.
Older Holstein cows reached t p earlier than their first-and second-order counterparts did (Table 4). Brown Swiss cows followed the same Holstein patterns, but stayed in second place as milk producers. Jersey and crossbred cows did not follow the same patterns presented by Holsteins and Brown Swiss cows. Second-and third-order Jersey cows produced more milk than the ones of the first parity order. Girolando cows were low in rank, but the dataset contained only first and second parity order cows (Tables 1 and 4).
Holsteins presented an increasing pattern for κ α from the first to the fifth parity order, but the rate reduced at the sixth parity order ( Table 5). The same happened for Jerseys and crossbreds from first to the fourth lactation, whereas Brown Swiss cows showed an increase from first to the third lactation. The second-order Girolando cows presented an ascending rate of milk production faster than their firstorder counterparts did. Despite the β ≤ 0 estimates for BPO 41 and BPO 43 , the convex profiles generated from positive random effects over the fixed β allowed the estimation of κ α (Table 5).
Holsteins and Brown Swiss cows of the first parity order generally reached the time at the inflection point later than their older counterparts did. The exception was the sixth parity order Holsteins (Table 5). However, this pattern was not followed by the other breeds and the crossbred cows. In fact, Jersey cows increased t i as they got second and third parity orders, but fourth-order Jersey cows reached the inflection point earlier.
The negative fractional descending rate was higher for first parity order cows, irrespective of breed or crossbreed (Table 5). Older cows generally presented lower descending rates; therefore, they did not sustain the peak milk yield as first parity order cows did. There was a negative association between κ α and the negative κ d (ρ̂ = -0.896, P<0.001). Given the values of the SE reported, we can depict that κ α and κ d estimates for the fixed effects were precise. The relative variation for the ascending rate (100SE κ α /κ α ) ranged from 0.33 to 8.97%, whereas the relative variation for κ d ranged from 0.14 to 3.95%. On the other hand, the estimates for λ presented a greater relative variation for the fixed effects, i.e., from 2.18% to 221%. The introduction of a random effect on parameter λ was ineffective to reduce the AICc m substantially if compared with the random effect ascribed to β (Table 2). Therefore, only positive values occurred for λ as population estimates (Table 3). The ascending and descending rates are nonlinear functions of the model parameters, and the integrated information resulted in precise estimates, which favored the direct comparisons; isolated parameters were less precise (Tables 3 and 5).
For comparison purposes, we regressed the estimated least-squares means for r d (w -1 ) over the estimated least-squares means for κ d (w -1 , Eq. 8) in their original scales, i.e., both rates are characteristic of the descending phase after peak milk yield and are both negative. Because of r̂d and κ̂d domains (estimates < 0), we assumed r d normally distributed, and we estimated the simple linear regression by least squares (PROC REG, SAS University Edition): r̂d = 3.5 • 10 -3 + 1.5 • κ̂d. Therefore, for the proper estimation of r̂d one must use -κ̂d in the preceding regression (Table 5). The SE of the intercept was 4.11 • 10 -4 , the SE of the slope was 3.19 • 10 -2 , the standard error of the regression was 4.19 • 10 -4 , and  The geometry of the lactation curve based on Wood's equation: a two-step prediction Oliveira et al.

13
the R 2 = 0.9925. The estimated values for r d and κ d were precise and close to each other. Regardless of breed, cows of first parity order presented smaller arbitrarily positive κ d estimates in the first half of the down ranked least-squares means contained in Table 5.
The relationships among parameters α, β, and λ with, for instance, κ α , κ d , t p , μ t p , and TMμ t are mathematical. Therefore, the functional relationships established among those parameters explain, by and large, the correlation estimates among those same parameters. Nonetheless, the Pearson correlation for TMμ t × μ t p in our study amounted to 0.982 (P<0.001) regardless of breed and parity order. This was done for those cows (n c = 6369) with an overall random estimate for β ∈ (0, 1). Nonetheless, the random β estimates for all 6754 cows presented a positive correlation with TMμ t , namely 0.873 (P<0.001). Therefore, generally, even for cows with atypical concave profiles, the greater the random estimates of β, the greater the predictions for TMμ t . In our data, some examples of the Pearson correlation coefficients were: -0.519 (P<0.001) for TMμ t × κ α , 0.035 (P = 0.005) for TMμ t × κ d , -0.828 (P<0.001) for κ α × κ d , and -0.896 (P<0.001) for κ α × t p .

Discussion
The technique of the nonlinear mixed-effects models allows the ranking of dairy cows within each breed and each parity order (ranks not shown). Therefore, one can replace the unstructured G matrix and ascribe a phenotypic value by associating to it genetic merits by the two-step estimation technique with a genomic relationship matrix in the linear mixed model Soares et al., 2017). An extension to the generalized linear mixed-effects model is possible (Littell et al., 2006;Vonesh, 2012; Table 5 -Attributes 1 of the lactation records obtained from the two-step estimation based on Wood's equation: specific ascending rate of milk production until peak milk yield (κ α , w −1 ), time at the inflection point (t i , w), specific rate of milk production decline post peak milk yield (κ d , w −1 ), their respective standard errors (SE), and lower (L) and upper (U) 0.99 confidence limits The geometry of the lactation curve based on Wood's equation: a two-step prediction Oliveira et al.
14 Stroup, 2013). Nonetheless, we remind the reader that we fitted the random effects of cow (σ cow,β ) and cow × parity order (σ cpo,β ) successfully only to the shape parameter β, and this result may explain the fitting flexibility of Wood's equation over the other equations studied here. In addition, the combined powers of equations 17 and 19 accounted for variations attributable to either scale and correlations among repeated measures (Tables 2 and 3). If not discounted, those variations can be understood as pseudorandom variations that otherwise would inflate the residual variation (Vieira et al., 2018). The nlme function of R accounts for all these issues and its main advantage is to integrate all information in a single variance-covariance matrix. The overall benefits are an improved precision of parametric estimates and information loss minimization during the quantitative interpretation of nonlinear phenomena with mathematical models (Pinheiro and Bates, 2000).

Nonlinear models and their quality of fit
The Akaike criterion (Akaike, 1974) corrected for small samples (Sugiura, 1978;Hurvich and Tsai, 1989;Cavanaugh, 1997), namely AICc m , has been reported as one of the most important measures to evaluate the predictive power of mathematical models in a multiple hypotheses framework (Buckland et al., 1997;Sober, 2002;Burnham and Anderson, 2004). The AICc m and its derived measures constitute the I-T approach, which allows the comparison of the quality of fit of different models and, whenever necessary, establishes parameters for model averaging. If one takes the smallest AICc m value to choose a model over the others, does an incomplete procedure within the I-T framework (Burnham and Anderson, 2004;Vieira et al., 2018;. The smallest AICc m value indicates the lower information loss due to the fit of a given model among the feasible models evaluated. However, if one assumes that a given model is the best solution and ignores the others, depending on the p m for the chosen model, its uncertainty can be large. This is why we assumed that a low uncertainty for a given model is 1 -p m < 0.10 ( Burnham and Anderson, 2004;Vieira et al., 2018).
The frequentist statistical tools available nowadays allow the association of random factors to the fixed parameters of nonlinear models and the fit of the resultant models by maximum or restricted maximum likelihood. In addition, one can challenge the traditional assumptions of homoscedasticity and independence of the errors with these tools based on the theories of nonlinear and generalized linear mixed-effects models (Pinheiro and Bates, 2000;Vonesh, 2012;Stroup, 2013). In our study, the use of the available statistical tools (nlme from R and GLIMMIX of SAS) provided considerable improvements in model fitting and prediction. Based on the I-T approach, we observed that Wood's equation was superior in terms of quality of fit and the best choice to represent the diversity in shape and form of the lactation profiles, given the data. Nonetheless, many studies use the R 2 to compare models, but this measure is not sensitive to differences among nonlinear models, which make comparisons useless because of the difficulty to demonstrate differences. Therefore, models that actually differ based on I-T may be considered equal if one bases on R 2 or, alternatively, its adjusted version for model selection (Spiess and Neumeyer, 2010;Hossein-Zadeh, 2016).

Lactation curve models
In this study, we challenged the ability of some models to mimic actual lactation records over time. Although the abundant dataset (37810 time-records) might have favored both accuracy and precision of the estimates, considerable variation may have arisen due to different ethnic types, parity orders, low number of cows for some breed-parity order interactions, and incomplete or irregular lactation records to characterize an entire lactation (0-305 days). These sources of variation might have constrained some of the results and increased the difficulty of obtaining valid, convergent, and robust solutions for some models. Nonetheless, the dataset offered the opportunity for Wood's model to exhibit its flexible nature to fit both convex and concave lactation profiles. Concave profiles are rather frequent than one might expect; these concave profiles are characteristic of animals that did not present peak milk yield, and equation 1 mimicked those profiles when other models (equations 12-14) failed. Within the framework of nonlinear mixed-effects models, Wood's model even was able The geometry of the lactation curve based on Wood's equation: a two-step prediction Oliveira et al. 15 to fit lactation profiles with incomplete records, which are needed to characterize an entire lactation trend (Table 3 and Figure 2, panel e).
Several authors faced difficulties with fitting models to atypical lactation profiles, that is, those profiles that deviate from the standard lactation convex-shaped trend (Rekik and Ben Gara, 2004;Macciotta et al., 2005;López et al., 2015). Therefore, the greater flexibility of equation 1 provides the necessary accuracy and precision to the selection of cows based on ranking their performances as parametric attributes of the lactation profiles such as peak milk yield, mean milk production, ascending and descending rates, and more importantly, the total milk production per lactation.
There was a positive agreement between random β and total milk production estimates. Therefore, the use of the peak milk yield as a parameter for selection may result in the exclusion of cows that may produce high quantities of milk, even though the shape parameter β ≤ 0 for a given breed. Our findings favored the use of all lactation profiles from all cows. They avoided an equivocated exclusion of animals before the fit of the model, i.e., all animals were evaluated regardless of the occurrence of a peak milk yield or not (Rekik and Ben Gara, 2004). The poor performance of the other models (equations 12-14) may be explained by the no prior arbitrary exclusion of concave lactation records, as well as incomplete and irregular lactation profiles from our dataset. The advantage of the nonlinear mixed-effects models relies exactly upon the joint analysis of all lactation records employing a matrix of random effects, which encompasses the variation from the random sample of cows about the fixed parameters.
There are some reports that present a poor performance for equation 1 if compared with equations 12 and 13. However, they used the nonlinear least-squares method of estimation to fit the models, did not consider random factors over the fixed parameters, nor did challenged the independence of residuals and homogeneity of variances (Brody et al., 1924;Dijkstra et al., 1997;Hossein-Zadeh, 2016). Generally, convex-shaped lactation records are selected because they resemble the typical, standard lactation profile (Morant and Gnanasakthy, 1989;Dijkstra et al., 1997;Macciotta et al., 2005). Nonetheless, the final shape and form of the nonlinear model depend on its parametric estimates (Wood, 1967;Grossman and Koops, 1988;Morant and Gnanasakthy, 1989;Dijkstra et al., 1997;Macciotta et al., 2005). The sign and magnitude of the β estimate determine if the curve is flatter, presents a sharp peak, or decline monotonically at the onset of lactation at parturition (Congleton Jr. and Everett, 1980;Morant and Gnanasakthy, 1989). In addition, it is important to note that the inflection point in the descending phase is part of a typical convex-shaped lactation profile (Druet et al., 2003).
Congleton Jr. and Everett (1980) and Macciotta et al. (2006) believed that the lack of points at the ascending phase do not allow the proper characterization of the typical convex-shaped lactation profiles with a characteristic peak milk yield. Nonetheless, our dataset contains concave-shaped lactation records with several time points, as the better-characterized convex-shaped profiles do. This is why the accommodation of the random factor to parameter β resulted in the superiority of equation 1 to represent the data, because final random estimates can be either 0 < β ≤ 1 and β ≤ 0. Because of the factors mentioned previously, we believe that those variations occurred when cows presented the peak milk yield on the day of parturition or nearly after parturition. These types of lactation profiles occurred for groups BPO 41 and BPO 43 , which might be associated with the poor dairy temperament of Zebu breeds and its crossbreeds with dairy breeds (Bangar and Verma, 2017). Sometimes, even cows that exhibited a typical convex-shaped lactation profile in a previous lactation fail to present a peak milk yield in a further lactation, because of environmental factors such as metabolic disorders, mastitis, other diseases, and nutritional or other management errors (Wood, 1968(Wood, , 1970(Wood, , 1972(Wood, , 1976(Wood, , 1980Macciotta et al., 2005;Hossein-Zadeh, 2016;Ahmed et al., 2019). Nevertheless, the cows that presented typical convex-shaped lactation profiles were those that produced more milk in the entire course of lactation than their counterparts did.

Geometric attributes of Wood's model
The ascending phase of equation 1 is represented by κ α and can be associated with the proliferation of secretory cells of the epithelial parenchyma in the mammary gland of cows (Capuco et al., 2001), dairy does (Knight and Peaker, 1984), and mice and rats (Dijkstra et al., 1997). The DNA content of the mammary gland is richer in the ascending phase, which results in greater RNA transcription and the consequent translation into enzymes responsible for the synthesis of milk components such as galactosyl transferase, fatty acid synthase, and acetyl-CoA carboxylase (Knight and Peaker, 1984;Boutinaud et al., 2004). The continuous milk production increase may be associated with the net number of secretory cells that result from the processes of cell mitosis or differentiation (Capuco et al., 2001;Dijkstra et al., 1997Dijkstra et al., , 2010. Homeorhetic mechanisms orchestrate changes in the mammary gland and other metabolic organs of the animals to provide the machinery and nutrients for milk synthesis by the secretory tissue. The net results are the changes observed for milk yields in the course of the lactation, including changes in other variables that dairy animals go through during the dry and transition periods (Bauman and Currie, 1980;Bauman et al., 1999;Bauman, 2000).
The smaller estimates for κ α of primiparous cows within each breed observed in our study may be explained by the incomplete development of their mammary glands. First parity order cows are younger and still grow after the first calving; therefore, there is a competition for nutrients between the processes of growth, gestation, and the subsequent lactation, which delays the growth of the mammary gland and its full secretory potential observed in a mature cow (Bauman and Currie, 1980;Bauman et al., 1999;Bauman, 2000). Indeed, the mammary glands of cows (Capuco et al., 2001) and dairy does (Knight and Peaker, 1984) continue to grow until and after the first parturition. The volume of secretory tissue in the mammary gland of second parity order dairy does is greater than that of the first parity order does, but the greater differentiation in the quantity of secretory parenchyma occurs during pregnancy (Fowler et al., 1990;Knight and Wilde, 1993). Those biological facts may explain the greater κ α estimates for multiparous cows and, consequently, their greater peak milk yields estimated in our study.
The peak milk yield can be used as a parametric reference for phenotypic selection of cows because of the good positive agreement between total milk production and peak milk yield in dairy cattle (Grossman and Koops, 1988;Tekerli et al., 2000;Hossein-Zadeh, 2014;López et al., 2015) and buffaloes (Hossein-Zadeh, 2016).
The descending phase occurs soon after the peak milk yield (Figure 1) (Wood, 1967(Wood, , 1976Knight and Peaker, 1984;Dijkstra et al., 1997Dijkstra et al., , 2010López et al., 2015). We used the inflection point at the descending phase to obtain an estimator for κ d . This specific constant rate can be an index of persistency: the greater the absolute values, the lower the lactation persistence, i.e., the sharper the decline after peak milk yield. Nonetheless, those cows that presented lactation profiles with greater κ α estimates also presented greater arbitrarily positive κ d estimates (smaller persistency). We observed a negative Pearson correlation between these two specific rates, which clearly demonstrated the relationship. Wood (1969) observed that as cows grow older, they start with greater milk yield records at parturition, but their rate of decline in milk production reduces as lactation advances post-peak milk yield (persistency reduces = larger absolute κ d ). Our results confirm that primiparous cows present flatter lactational profiles, which means a greater persistency if compared with their older counterparts (Grossman and Koops, 1988;Bauman et al., 1999;Cobuci and Costa, 2012;Hossein-Zadeh, 2014, 2016López et al., 2015). Therefore, the geometrization of equation 1 yielded parameter estimates that agreed with the biological trends reported in the literature.
The milk production decline post-peak yield may be associated to the reduction in the numbers of secretory cells in the mammary glands of cows (Capuco et al., 1997;Cobuci and Costa, 2012), dairy does (Knight and Peaker, 1984;Quarrie et al., 1994;Dijkstra et al., 1997), and small rodents (Walker et al., 1989;Strange et al., 1992;Dijkstra et al., 1997). The number of epithelial secretory cells recorded at the peak milk yield can be reduced by 40% until the end of lactation in cows (Capuco et al., 2001) and by 30 (Knight and Peaker, 1984) to 40% in dairy does (Fowler et al., 1990). Capuco et al. (2001) reported that the number of secretory cells results from the balance between cell proliferation and death, i.e., the quantity of secretory tissue reduces as the cell death rate overcomes the cell proliferation rate. This resultant cell loss occurs by some sort of a programmed mechanism for cell death also called cell apoptosis (Walker et al., 1989;Strange et al., 1992;Quarrie et al., 1994;

Author Contributions
Conceptualization