A meta-analysis of the relationship between milk protein production and absorbed amino acids and digested energy in dairy cattle.

Milk protein production is the largest draw on AA supplies for lactating dairy cattle. Prior NRC predictions of milk protein production have been absorbed protein (MP)-based and utilized a first-limiting nutrient concept to integrate the effects of energy and protein, which yielded poor accuracy and precision (root mean squared error (RMSE) > 21%). Using a meta-data set gathered, various alternative equation forms considering MP, absorbed total essential AA (EAA), absorbed individual EAA, and digested energy (DE) supplies as additive drivers of production were evaluated, and all were found to be superior in statistical performance to the first limitation approach (RMSE = 14–15%). Inclusion of DE intake and a quadratic term for MP or absorbed EAA supplies were found to be necessary to achieve intercept estimates (non-productive protein use) that were similar to the factorial estimates of NASEM. The partial linear slope for MP was found to be 0.409, which is consistent with the observed slope bias of −0.34g/g when a slope of 0.67 was used for MP efficiency in a first-limiting nutrient system.


INTRODUCTION
Perhaps due to the historical pricing schemes, most dairy nutrition models have only provided predictions of milk volume (NRC, 1966(NRC, , 1971(NRC, , 1978(NRC, , 1989(NRC, , 2001).Thus, one must assume milk composition to assign an economic value to the predicted milk production, and changes in composition that may be induced by the new ration are not prospectively considered.This is perhaps a greater deficiency for milk fat production given relatively high variation in milk fat concentration than for milk protein concentration, but the latter does vary, and the variation is economically important.Additionally, the biological order of operation is reversed when predicting milk volume first because milk volume is a function of the osmotic pull of the milk components A meta-analysis of the relationship between milk protein production and absorbed amino acids and digested energy in dairy cattle.
with most of the pull being generated by lactose (Jenness, 1974).Because lactose is synthesized from glucose by the enzyme lactose synthase which is composed of galactosyl-transferase and α-lactalbumin (Ebner and Schanbacher, 1974) with the latter being rate limiting (Jenness, 1974), lactose concentration and milk volume are driven by milk protein production which explains the low relative variation in milk protein concentration.Milk fat does not have significant osmotic potential, and thus it is much more variable in concentration.Predicting milk components first and using those to predict milk volume follows the biological process and is thus a more appropriate prediction approach.Additionally, that approach provides predictions of the primary value components of milk for use in nutrition models so that diets can be optimized to maximize product value minus input costs.
Several empirical and mechanistic models of milk protein production have been developed (Baldwin et al., 1994, Maas et al., 1997, Hanigan et al., 2001, Hristov et al., 2005a, Vyas and Erdman, 2009, Lean et al., 2018, Lapierre et al., 2020, Lapierre et al., 2023).Early efforts followed the conceptual scheme of Mitchell and Block (1946) where a single AA was assumed to be the sole limiter of milk or milk protein production.This approach was introduced for absorbed Lys and Met in NRC (2001).However, such an approach was shown to be an inappropriate representation at the tissue level (Hanigan et al., 2000).Over the past 20 years, the role of regulation of ribosomal protein synthesis has been shown to be an important determinant of protein synthesis (Davis et al., 2002, Arriola Apelo et al., 2014b, Luo et al., 2018).One of the signaling pathways involved in regulation of protein synthesis, mTOR, responds independently and additively to several EAA, energy, and hormonal stimulation (Rius et al., 2010b, Appuhamy et al., 2012, Appuhamy et al., 2014b, Arriola Apelo et al., 2014c, Castro et al., 2016b).Protein synthesis is also regulated by phosphorylated eIF2α, which is mediated by general control nonderepressible2 kinase in response to EAA supplies (Anthony et al., 2004).Thus, the concept of a single limiting AA or nutrient is not supported by the observational data.Additional work has demonstrated that milk protein gene expression also changes in response to AA and energy supply (Toerien et al., 2009, Rius et al., 2010a, Dessauge et al., 2011).
Although it has been demonstrated that representation of the independent effects of each AA provides a better prediction of milk protein production at the tissue level (Hanigan et al., 2000), such an approach has not been widely adopted at the whole animal level (Lean et al., 2018).Certainly it represents a major conceptual shift from application of the first-limiting nutri-ent approach (Mitchell and Block, 1946) with AA supplies expressed as proportions of MP or total absorbed EAA.The latter may have contributed to the lack of progress as that representation mathematically creates non-linearity in the independent variables, interactions among AA, and variation in the responses to each at a given AA supply due to variation in the denominator.Adding or subtracting a single AA results in changes in the ratios of all other AA given the change in MP or EAA supply elicited by the addition or subtraction.Changes in the ratios of the non-manipulated AA represent false signals as their supply was not altered by the single AA addition.
Following from the cell and tissue work, the objectives of this work were: 1) to assess the existing framework for prediction of milk protein production; 2) to develop and assess alternative, static, deterministic predictions of milk protein production using predicted absorbed nutrient supplies and observed animal characteristics as independent variables.It was hypothesized that the prediction model would be significantly improved when the supplies of digested energy (DE) and absorbed EAA were considered, and that expressing absorbed AA in absolute, additive terms would be superior to expressing relative to total absorbed EAA, MP, or DE.

MATERIALS AND METHODS
The work consisted of a mixed-model, meta-regression analysis of literature data.All data handling and statistical analyses was completed using R (v. 4.0.3;R Core Team, 2020) within R Studio (v 1.3.1093).

Meta Data
Data used for the analyses was assembled from data sets collected by the NRC (2001) committee with enhancements and additions by Bateman et al. (2008), Roman-Garcia et al. (2016), and Feng et al. (2018), the mammary data described by Hanigan et al. (2002), the total-tract fatty acid digestibility data of Daley et al. (2020), and data collected from studies describing the effects of infused AA and ruminally protected AA identified by a comprehensive literature search.In the latter case, the key words AA, dairy, cattle, and lactation were used to search within the Agricola and PubMed databases in 2019.Duplicates of studies already contained in the prior databases were excluded.Data were entered into an Excel-based database prepopulated with a single data format.In some cases, more than one experiment was described within a publication.In this case, each experiment was assigned a unique publication ID (PubID) assuming the experiments were From the total observations, 933 treatments from 221 experiments reported milk protein production or it could be derived from data that were reported, and were thus selected.Three treatments examining the effect of infused insulin were dropped.Seven treatments testing the effects of DL-2-hydroxy-(4-methylthio)-butyrate and its isopropyl ester derivative were removed due to uncertainty regarding the bioavailability of those Met analogs.
Dietary nutrient inputs were derived from the reported dietary ingredients for each treatment using the feed library from NASEM (2021) with adjustment for reported ingredient and dietary nutrient composition as described by Hanigan et al. (2021) and Fleming et al. (2019b).Three studies containing 12 treatments were removed because dietary nutrient composition derived by summation of the ingredient nutrients were more than 3 SD from the reported dietary nutrient concentrations suggesting data reporting may have been inaccurate in the publication.After initial work, 1 experiment containing 4 treatments was removed, as the reported milk protein production was identified as an outlier based on variance in the residual errors generated from model predictions.Two additional treatments from another experiment were removed due to missing input values, yielding 905 treatment means from 216 experiments that were used for the analyses.A listing of the studies describing the experiments that were used is provided in supplemental Table S1.
746 of the 905 treatments were from purebred Holstein, 7 from British Friesian, 4 from Jersey, 74 from Ayrshire, 4 from Brown Swiss, 8 from Polled Red or Red and White, and 4 from Pie Noir breeds.The remainder were Holstein crosses with Jersey, Simmental, or Brown Swiss breeds.
Milk protein was reported on a crude basis for 79 experiments, and the basis was unspecified for 76 experiments.In the latter case, 55 of the experimental descriptions indicated the analyses were conducted by infrared, generally at commercial laboratories.The milk analyzers in the US DHIA laboratories were calibrated to milk CP before the year 2000, and converted to a true protein (TP) basis thereafter (D.Barbano, personal communication).Therefore publications before 2000 were assumed to have been reported on a CP basis if unspecified, and thereafter on a TP basis.After characterization of the basis of measurement, all CP values were converted to a TP basis assuming TP was 95.1% of CP (DePeters and Cant, 1992).
Five experiments containing 16 treatments did not report milk fat concentration or production.These were filled with the mean milk fat concentration from the meta-data which was 3.60% ± 0.51.
Parity was defined as a continuous value based on the reported animal composition of the treatments with nulliparous animals coded as 0 (unused for the work herein), primiparous animals coded as 1, and multiparous animals as 2. Twenty-nine experiments containing 110 treatment means did not report parity.For those treatments, 1 was entered for treatments with reported BW less than the mean BW for primiparous animals (525 kg), 2 for treatments with reported BW greater than the mean BW for multiparous animals (610 kg), and the mean group parity (1.842) for the remaining missing values.
Body weight was not reported for 43 experiments representing 157 treatment means.Missing BW for these treatments were predicted using a regression equation developed from observed BW and reported parity (before replacement of missing parity values): .

NRC 2001 Evaluation
The NRC (2001) NEL and MP allowable milk protein predictions and the first-limiting nutrient approach were assessed as described in the Supplemental Materials (https: / / doi .org/ 10 .3168/jds .20XX-XXXXX).

Regression Techniques
Regression analyses were conducted using the lmer function of the lme4 package (R Development Core Team, 2020).Mixed models were fitted to the data by maximum likelihood with experiment ID (PubID) denoted as a random effect.The SE reported in many of the older publications were incorrect due to the use of Proc GLM from SAS which did not model covariance (Littell et al., 1998).This primarily affected the SE estimates for repeated measures and random effects models, although even the early estimates from studies that used neither were incorrect, and in the absence of covariance information, cannot be corrected.Therefore, the data were weighted by the square root of the number of observations represented in each treatment mean, which generally captures differences in experimental power, but does not fully capture the power of repeated measure designs nor reflect differences in precision among laboratories (see discussion by Hanigan et al., 2021).
Model performance was summarized by Akaike Information Criterion (AIC).Because the number of observations was large and fixed in size, there was no need to correct AIC for sample size, and it was thus calculated as: where N parms represented the number of parameters in the model, and serves as a penalty for more complicated models.Differences among nested equations were assessed statistically using the anova function of the stats package in R as described by Chambers and Hastie (1992).
The concordance correlation coefficient (CCC), root mean squared error (RMSE), and its partition were calculated from the residuals unadjusted for random study effects (provided by the predict function of lme4) as described by Lin (1989) and Bibby and Toutenberg (1977), respectively.The unadjusted values provide an estimate of the expected performance of the models on novel data where random study effects cannot be estimated a priori, which is the case for all applications of the model outside of the current data.RMSE and CCC were also calculated for the adjusted data (generated by the predict function of lme4 using the full mixed model solutions) where an estimate of the contribution of random effects could be calculated by difference.
A subset of the candidate models that were biologically consistent, plausible, and ranked high based on AIC and the RMSE analyses were further assessed for parameter stability across solves, potential collinearity problems, and parameter stability across subsets of the data.The latter was achieved using exhaustive Monte-Carlo based, cross-evaluations, which consisted of 500 iterations of fitting the model to random selections of 85% of the data, and evaluation of the resulting model using the remaining 15% of the data.This was carried out using the caret (v.3.0-8) and mlbench (v.2.1-1) packages of R. Data selection for the split between training and evaluation was completely random with no consideration of study.Means and SD of model parameters and model fit statistics were calculated from the results of the 500 cross-evaluation runs.Variation in parameter estimates across solutions less than 50% of the mean were denoted as stable and reliable, between 50% and 70% as moderately stable and reliable, and variance greater than 70% as very unreliable and not well supported by the data.These categories were supported by probabilities of the estimates differing from 0 where parameters with CV less than 50% generally had mean probabilities less than 0.05; CV between 50 and 70%, had mean probabilities generally ranging from 0.05 to 0.20; and CV greater than 70% having mean probabilities generally greater than 0.20.

Global Model Forms
Preliminary work to identify the best global model (those containing all of the terms considered for inclusion) forms was conducted using additive linear and quadratic terms for CP intake (CPIn), MP intake (MPIn), total absorbed EAA intake (EAAIn), DEIn and their interactions, and this was subsequently extended to consider deaggregated protein and energy terms.Based on the preliminary evaluations and the study objectives, 3 global mixed effects model equation forms were formulated using different methods of Hanigan et al.: Predicting Milk Protein representing the effects of absorbed AA on milk protein synthesis.They are described in the following sections.
Although DIM was found to be generally significant in those equations, only 5 studies had variation in DIM among treatments within the study (Polan et al., 1991, Schwab et al., 1992, Choung and Chamberlain, 1993, Volden, 1999, Samuelson et al., 2001); but, none of the studies were designed to test for stage of lactation effects and thus all were confounded by time and diet.Inclusion of DIM in the current models resulted in highly non-biological intercepts (maintenance requirement) suggesting it was confounded with study effects that could not be addressed in the absence of any within study comparisons.Year of publication was also found to be significant with approximately a 6 g/year gain in milk protein output, however, inclusion of the factor did not affect other parameters, and because no data could be included for the future, models retaining the term would result in predictions outside of the range of data from the present going forward, which is of concern.Thus, DIM, and year of publication were excluded a priori from the global models.
Absorbed EAA Expressed in Grams per Day.The 1st global equation form utilized AA inputs expressed as grams of absorbed AA/d and considered both linear and quadratic effects for each EAA.InfEAAIn Bld was added to allow for potential differences in responses to AA absorbed from the intestine versus jugular or arterial infusions.There were 6 studies that utilized vascular infusions; all were delivered post-splanchnic.In this case, first-pass use by the gut and the liver during absorption is bypassed.This may result in a difference in the post-splanchnic supply of 5 to 10% (MacRae et al., 1997, Hanigan et al., 2004, Fleming et al., 2019a).Following initial work that demonstrated that most responses to individual absorbed AA were linear, the in-dividual squared AA terms were collapsed into a single summative quadratic term ⅀(EAA 2 (g 2 /d) which was the sum of the squared individual absorbed supplies of each EAA (a = 1 to 10) reflecting the potential of overall saturation of protein synthesis or its regulation by the combination of EAA and assumes a common effect size among the EAA: 1 . [0.4] Absorbed EAA Expressed as a Proportion of Total EAA.Historically, EAA requirements and supply have been expressed as a percentage of MP or of total EAA.This approach implicitly assumes the primary driver of production is MP or total EAA, and the EAA must have an ideal composition to achieve the full animal potential at that level of MP.The global model for that approach using total EAA as a denominator for calculation of percentages was:  [6]

Model Selection
Initial model exploration was conducted by backward elimination using the step function of the lmerTest package (v.3.1-2) of R. Terms were sequentially removed from the global models until all terms had P values less than or equal to 0.05.Because backward elimination can result in biased models (Dunkler et al., 2014), final model selection was conducted using an all-models approach encoded in the dredge function of the MuMin package (v.1.43.17) of R. With this method, all combinations of terms in each global model were assessed, and the models ranked based on statistical performance.This resulted in a maximum of 28 fixed-effect terms (10 EAA, 10 EAA 2 , NEAA (g/d), InfEAA Bld , 4 energy intakes, BW, and Parity) equating to 536 million possible combinations to be evaluated and ranked.Because such a large matrix could not be held and manipulated in memory by R even with 4 Tb of memory, each of the global models were first evaluated to identify fixed terms that could be eliminated based on backward elimination so that the global models could be reduced to no more than 22 terms which equates to 4.2 million combinations.
Final equations were primarily presented in tabular form.Tables were constructed using the tab_model function of the sjPlot package (v. 2.8.4)

RESULTS AND DISCUSSION
A summary of the meta data is provided in Table 1.In general, the data had a large range spanning milk production from 10.2 to 53.8 kg/d with DMI ranging from 9.1 to 31.8 kg/d.There were 2 limitations to the data which were the general lack of characterization of pregnancy state and stage and the limited reporting of BW change.Studies often reported ending BW and compared differences among treatments, but the beginning BW were not generally reported.These deficiencies prevented meaningful assessments of the impact of BW change and gestational nutrient use on production and efficiency.The magnitude of that potential contribution to energy and protein supplies was examined by (Liu et al., 2021).
All dietary nutrient intakes were normally distributed with a minimum range of 5 SD, but correlations among the inputs were generally high (>0.75;Supplemental Figure S1; https: / / doi .org/ 10 .3168/jds .20XX-XXXXX) likely reflecting the scaling of all nutrients with milk production.Much of the range in EAA inputs was driven by trials which fed rumen-protected AA or infused AA (Figure S2 and Figure S3).Although the infusion of proteins helps define the curvature at the top of the production range, the EAA within that infused protein are fixed in ratio and thus such treatments do not help address collinearity problems among individual EAA responses.Methionine absorption had the lowest correlations with other EAA ranging from 0.69 to 0.77 due to the large number of studies which explored the effects of methionine infusions and feeding ruminally protected methionine.The relatively large correlations among all of the other EAA underscores the value of infusion studies making use of single or small groups of EAA and feeding trials utilizing well characterized rumen-protected AA sources.
Digested energy intake and MPIn suffered from the same collinearity problem having a correlation coefficient of 0.83.Digested energy is composed of the energy contributed from each of the nutrients including MPIn, thus the 2 terms are mathematically related.Subtraction of the DE arising from MPIn from the DEIn term reduced the correlation between DEInp and MPIn to 0.69.The reduction was similar for the relationship between DEInp and total absorbed EAA, and between DEInp and each of the individual absorbed EAA (Figure S1).The statistical benefits of removing the subnutrients of interest from the aggregated term was demonstrated for microbial protein flow predictions (Hanigan et al., 2021).Absorbed supplies of each EAA were correlated with milk protein output, but the efficiency of conversion of MPIn to milk protein ranged widely and is very poorly related to the proportion of each EAA relative to total EAA (His, Lys, and Met relationships shown in Figure 1; the remainder not presented).MPIn conversion efficiency is also poorly related to each EAA expressed relative to DEIn (data not shown).Despite tremendous range in milk protein production per unit of MPIn or DEIn, there was no evidence of a break point for DEIn and only weak evidence for MPIn (Figure 1).Such points have classically been used to define a requirement.The apparent quadratic responses to Lys and Met were due entirely to a small number of rumen-protected AA feeding trials and infusions, and the range in milk protein within the mass of observations in the very middle of the input range was much greater than the predicted range based on the regression equations (Figure 1).The lack of a clear relationship was particularly apparent when expressing EAA supplies as a ratio with DEIn.Thus, the data offer little support for modeling milk protein responses as a function of the proportion of individual EAA in total EAA or as a ratio to DEIn.

NRC 2001 Model Performance
From Occam's Razor, models should generally be no more complicated than necessary, and thus an evaluation of the NRC (2001) model and its extension based on the first-limiting nutrient concept (Mitchell and Block, 1946) was undertaken as a reference point.Because the model does not output predictions of milk protein, these were generated by reversing the calculation of MP allowable milk using the observed milk protein concentrations for each observation (Table 2 and Figure 2).The resulting predictions show clear systematic bias as previously reported (White et al., 2016b) with the model generating predictions generally centered in the data for NEL and MP allowable, but underpredicting by 91 g/d when considering the minimum of the 2 and by 202 g/d when considering the minimum of NEL, MP, and all 10 EAA.All of the predictions exhibited negative slope biases ranging from −0.44 g/g to −0.24 g/g indicating the model greatly overpredicted responses to changes in the supply of those nutrients.Although the slope bias marginally improved as more nutrients were considered in the first limiting nutrient framework, the increase in mean bias overwhelmed those gains resulting in predictions with much poorer accuracy and precision as measured by RMSE and CCC.Thus the model has structural or parameter problems that preclude unbiased estimates of milk and milk protein production, and the problems do not derive from the nutrient supply equations as the feed library and the nutrient supply predictions used herein were demonstrated to be largely free of mean and slope bias (White et al., 2016a, White et al., 2017, Fleming et al., 2019a, Fleming et al., 2019b, Daley et al., 2020, Hanigan et al., 2021).Hanigan et al. (2000) observed the same mean bias problem when the 1st limiting nutrient concept was applied at the mammary tissue level, suggesting the concept lacks validity.
The slope bias associated with MP allowable protein indicated the marginal response of milk protein to MP input was approximately 40% rather than the 67% used in the model.Because slope bias was reduced, but not eliminated, when the EAA were added to the calculation of first-limiting production, a portion of the problem may be due to the effect of dietary AA on milk protein production.However, the simultaneous introduction of mean bias as the concept is extended across nutrients clearly indicates that the form of expression is not reflective of the underlying biology.Adding a representation of a mechanism underpinning the aggregated function should never increase the error.
The poor performance of the NRC ( 2001) model should not be surprising.It has been clearly shown that marginal recovery of MPIn decreases as MPIn increases (Doepel et al., 2004, Metcalf et al., 2008).Therefore, usage of a fixed efficiency of utilization of MP or EAA directly leads to this type of slope bias (Lapierre et al., 2020).In addition, evidence began accumulating that milk protein production responded independently to multiple EAA simultaneously (Schwab et al., 1976, Clark et al., 1978).Hanigan et al. (2000) demonstrated the abject failure of the single-limiting EAA approach for predicting milk protein output from EAA removed by the mammary glands, and the clear superiority of representing the effects of each EAA in an additive form.This concept was reiterated by Cant et al. (2003), and a plethora of in vivo and in vitro cell signaling work has demonstrated that milk protein synthesis is highly regulated by at least insulin and IGF-1 concentrations, energy supply, and Ile, Leu, Lys, Met, and Thr with strong trends for several additional EAA (Burgos and Cant, 2010, Burgos et al., 2010, Rius et al., 2010c, Appuhamy et al., 2011b, Appuhamy et al., 2012, Appuhamy et al., 2014a, Arriola Apelo et al., 2014a, Liu et al., 2017, Yoder et al., 2020).However, the regulation has not always been observed (Cant et al., 2001, Doelman et al., 2015a, Doelman et al., 2015b, Nichols et al., 2017),.

Protein and Total EAA Based Predictions
Given the evidence that the current framework is inadequate, work was undertaken to understand the general milk protein response surface with respect to MP and DEIn.Given the complexity of considering the effects of 10 or more individual AA responses, understanding the general response surface was deemed important.
Although CP-based models (NRC, 1971(NRC, , 1978) ) are considered inferior to absorbed protein (NRC, 1989), and MP-based models (NRC, 2001), the improvement in precision across that spectrum is modest (Table 3).Stimulating milk protein output in the absence of a change in the EAA supply necessitates an improvement in post-absorptive EAA efficiency.The meta-analyses of Hanigan et al. (1998Hanigan et al. ( , 2002) ) and experimental work of Bequette et al. (2000) demonstrated that mammary EAA transport activity responds reciprocally to supply, increasing activity when supply is low, and the reverse.Subsequent work by Fleming et al. (2019a) demonstrated that the liver acts in the opposite manner, increasing transport activity when arterial EAA concentrations are high, and the reverse when low.Acting in concert, the 2 tissues can shift more of the EAA supply to mammary use when milk protein production is stimulated by a greater supply of energy providing the needed flexibility in net supply to mammary tissue as DEIn changes (Vanhatalo et al., 2003a, Raggio et al., 2006, Rius et al., 2010a).Although the overall relationship between MPIn and milk protein production appears generally linear (P < 0.01, Figure 1), there were significant nonlinear relationships for all 3 approaches (CP, MP, and EAA; Table 3).This was not due solely to high leverage points at the upper range of the data, as the nonlinear terms cross-evaluated well with repeated, random removal of 15% of the data (discussed below).Representing the relationship between milk protein production and MPIn as a linear relationship resulted in significantly poorer prediction performance (eqn.[3.2] vs eqn.[3.3]), and, perhaps more importantly, a non-biological intercept (419 ± 28 g/d).
The negative value of the intercept should reflect nonproductive net protein (NP) requirements of a fasted animal, and thus the intercept would be expected to be significantly negative for all animals.Inclusion of both the linear and squared terms for MPIn (eqn.14.9 14.9 14.5 14.9 Mean Bias, %  , 2021).Therefore, the efficiency of conversion at maintenance intake would be expected to be less than at production level intakes due to the greater energy density in diets for lactating cows as compared with dry cows.Additional work is needed to understand this inconsistency.Currently, the factorial method is used to represent non-productive use which includes metabolic fecal protein, and a regression equation with a much smaller intercept is used to represent MP use for milk protein.
Adding BW centered to the mean of the data to account for different sized animals (eqn.[3.6]) resulted in a slope coefficient of −0.35 g/kg BW, a 13 g/d decline in the intercept, no real changes in the other slope coefficients, and no improvement in overall fit statistics.The lack of an improvement in fit statistics with the inclusion of BW in the regression would normally indicate the added complexity should be excluded, however, as the model should be as widely applicable as possible, including for small breeds such as Jerseys, and the latter are very poorly represented in the development database, it seemed prudent to retain the term given its significance and general biological consistency.At a BW of 700 kg (Holstein breed average; NASEM, 2021), the non-productive NP requirement estimate from eqn.The linear effects of energy intake and the curvilinear responses to protein are consistent with declining marginal returns to MP supply as supply increases (Whitelaw et al., 1986, Hanigan et al., 1998, Kalscheur et al., 2006, Arriola Apelo et al., 2014b, Daniel et al., 2016), which presumably reflects saturation of the capacity of the mammary tissue to synthesize milk protein at a given energy supply.However, the latter responses are within a given set of animals or cells whereas the overall response captured by a squared MPIn term represents a global response to MPIn across the experimental population of animals.The apex of MPIn responsiveness was predicted to occur at 3,352 g MPIn/d from the linear and squared coefficients for eqn.[3.5].This yields a maximum predicted milk protein of 1,511 g/d.Thus animals with low potential to produce milk protein, such as the cows of Whitelaw et al. (1986), would be predicted to exhibit linear responses to MPIn given low MPIn intakes while very high merit animals would be predicted to operate above the apex predicted with this equation.Although the squared term captures at least part of the nonlinear responses in the cows within this data set, it cannot be expected to be globally applicable to the full genetic range of animals because those of high genetic merit will have greater synthetic capacities than represented by the data herein and the resulting equation.Conversely, the cows of Whitelaw et al. (1986) did not exhibit linear responses because their genetic capacity was low.
We also observed independent quadratic responses to RDP and RUP intakes when replacing MPIn in the equation.The slope coefficients were, however, seemingly greater for both than reported in NRC (2001) at 0.34 for RDP intake g/g and 0.24 g/g for RUP intake, likely because the prior estimates would have been colinear with the DMI term used in the model.At a mean DMI of 20.7 kg/d, these estimates equate to a zenith for milk protein at 15.9% RDP and 11.2% RUP which are both greater than observed by NRC ( 2001), but are consistent with the generally linear response of microbial protein to RDP (NASEM, 2021) and milk protein to MPIn herein.(Orskov, 1992) hypothesized that maximum protein synthesis was related to energy supply to the tissue.As energy is tightly regulated, the partial efficiency of energy use for milk production has been constant over time (Moraes et al., 2015), thus increases in genetic gain in milk production are associated with increased energy intake.If energy supply is related to maximum rates of protein production, this would present as an interaction between energy and protein terms as observed by Hanigan et al. (1998).Unfortunately, we were not able to identify any interactions between DEInp and MPIn including the form used by Hanigan et al. (1998).Any addition of higher order DEInp terms resulted in nonsignificant estimates for those terms, and an increase in AIC.These findings are consistent with the lack of interactions found at the cellular (Clark et al., 1978, Appuhamy et al., 2011b, Appuhamy et al., 2014a) and animal levels (Vanhatalo et al., 2003a, Raggio et al., 2006, Rius et al., 2010c).Although there must be a maximal rate of signaling and protein synthesis that can occur, the energy and hormonal arms of the signaling pathways do not appear to dictate that plateau within normal physiological ranges (Arriola Apelo et al., 2014c, Arriola Apelo et al., 2014d, Castro et al., 2016a), and energy supply cannot be used to scale maximum protein production to accommodate very high levels of production observed on North American dairy farms.It remains to be seen whether mathematical scaling based on pen or herd production will suffice (NASEM, 2021).
Replacing DEInp with DE concentration in DM (data not shown) resulted in much poorer predictions as compared with eqn.[3.5] with a non-significant intercept.Attempts to express MPIn as a ratio to DEIn, as in some swine models (NRC, 2012), also resulted in fits that were poorer than eqn.[3.5] based on AIC.The latter approach requires careful exploration of the ratio to ensure it is normally distributed and it is linearly related to the dependent variable (Lien et al., 2017).Lien et al. (2017) reminds us that if the ratio is linearly related to the dependent variable, its reciprocal must be nonlinearly related, and a nonlinear relationship with the independent variable violates the linear regression assumptions.Use of a ratio also forces a fixed relationship between the 2 variables rather than letting the relationship be discovered via combinations of linear and higher order terms for each.4. As for the MPIn-based models, DEInp was found to be highly significant for all of the models, as were the partitioned energy terms (DEStIn, DENDFIn, DEFAIn, and DErOMIn; Mcal/d).In no case, did the effect of systemic infusions (primarily jugular or mesenteric vein infusions) remain in the model after the selection process suggesting that losses during absorption from dietary sources were not significantly greater than from direct systemic infusions given the errors of measurement.4).The mixed effects model had an AIC of 10,360 and a conditional R 2 of 0.87.The fixed effects portion of the model had a marginal R 2 of 0.51, and RMSE of 131 g/d, which was 14.2% of the observed mean, no mean bias, and only 4% of the MSE associated with a slope bias.As the mixed model was fitted to the data, the latter reflects some structured random effects not captured by the fixed effects terms.This equation form has generally been used in the past as it was thought to approximate the saturation effects of each EAA, and the curves could be linearized to identify a breakpoint which was deemed the requirement (NRC, 2001).However, the negative linear responses for Thr and Val are difficult to explain biologically.The negative effect of Val may be an indication of an antagonism of Ile and Leu, although neither were present in the final model.If Val were truly antagonistic, one would expect a positive slope for at least one of the other 2 branched-chain AA (BCAA).It is also possible that the result is due to co-linearity among variables.The correlations between absorbed Val and Ile and Leu were high at 0.90 and 0.92, although several treatments deviated significantly in the relationship between Leu and Val (Figure S1).There was less variation between Ile and Val, but the low correlations among the linear parameter estimates from the regression model (generally less than |0.10| with only His exceeding that range at −0.38) do not support that argument (data not shown).The response might also be due to lack of independent variation from total EAA with Val having the highest correlation at 0.99 (Figure S1).Thus, the source of the negative response is unclear, but it seems likely to reflect an autocorrelation problem, and thus is biologically suspect.

Absorbed EAA Expressed as Gram per
The presence of a negative linear term for Thr is even more difficult to explain.The zenith of the quadratic occurs at an absorbed Thr supply of 82 g/d with a partial contribution to milk protein of −759 g/d.The partial contribution only reaches 0 at just under the maximum prediction for absorbed Thr, and thus it has a large negative effect throughout the expected dietary range.This is the opposite of that observed by Arriola Apelo et al. (2014d), and thus the response is not supported by direct experimental evidence at the tissue level.But Thr has the 2nd highest correlation with the total absorbed EAA at 0.98, and correlations among parameter estimates were generally greater for Thr with the correlation between Thr and Phe being −0.63.This may explain the relatively large positive slopes for Arg and Phe.Therefore, the negative linear slope may also represent an auto-correlation problem rather than a biological response.Finally, all 4 of the EAA having linear and quadratic terms exhibited very high correlations (<-0.97) between the linear and quadratic terms within an EAA.This is a clear indication that the higher order terms were not independently defined by the data, and thus contributing to the overall model definition problem, and likely the cause of the negative linear term for Thr.Responses to Phe were predicted to plateau at 103 g/d with a partial contribution to milk protein of 275 g/d.The partial contribution of Trp to milk protein was 280 g at the minimum absorbed supply of 9 g/d (51% gross Trp efficiency based on milk protein composition of NASEM ( 2021)) and plateauing at 346 g/d with an absorbed supply of 21 g/d (27% gross Trp efficiency).The partial contribution of absorbed Met to milk protein was predicted to be 27 g/d at the mean absorbed Met value of 46 g/d, yielding a partial response of 0.59 g milk protein/g absorbed Met for a mean gross Met efficiency of 1.8%.Such a low response is not supported by the direct experimental evidence of Pisulewski et al. (1996a) who observed a linear response from 1010 to 1100 g milk protein to infused Met ranging from 0 to 24 g/d.This reduces to a slope of approximately 3.75 g milk protein/g absorbed Met, and equates to a gross absorbed Met efficiency of 11.3%.Plateau values were not observed for His, Lys, or Met offering no information on target requirement levels under the classical approach.From the quadratics, plateaus for Arg, Phe, Thr, and Trp occur at 130, 104, −82, and 21 g absorbed/d, respectively.Obviously the Thr response and plateau are non-biological.Arriola Apelo et al. (2014d) observed that in vivo concentrations fell quite low on the linear portion of the quadratic curves, thus, the lack of quadratic effects for some EAA was likely due to lack of range in the in vivo data.
Of the terms in the reduced model, the linear effects of His, Trp, and Val were marginally defined by the data having CV of 59, 49, and 49% respectively with SD of the CV of 77, 20, and 27%, respectively under exhaustive cross-evaluation (data not presented).Surprisingly, the quadratic terms were all fairly well stable despite the high correlations with the linear terms.Regarding His, Trp, and Val, the large CV are indicative of the parameter estimates being a function of a relatively small portion of the overall data, which minimally suggest more data are required, but it could also be an indication of heterogeneity in the responses for some data sets.
Although eqn.[4.1] performed relatively well on crossevaluation, suggesting that the solution was stable, the large linear slope values for Arg, Phe, and Trp; the negative linear slopes for Thr and Val; the exponential response for Met; and the high correlations among some parameters suggests the model was overparameterized and likely fails to capture the true response, which argues against its use for field application.
Absorbed EAA Expressed as Gram per Day with a Common EAA Quadratic.To address the apparent overparameterization problem of eqn.[4.1], and recognizing that at least mTOR signaling is integrated across AA, energy, and hormonal signals yielding a com-mon regulatory response (Castro et al., 2016a, Sabatini, 2017), we hypothesized that protein synthesis capacity could also be represented with a combined saturation point.Thus the individual quadratic terms for EAA in eqn.[0.3] were replaced with a common quadratic slope term applied to the sum of the individual squared EAA yielding eqn.[0.4] (Table 4).This form has the implicit assumption that the saturation response to individual EAA is additive and the responses to each EAA are similar.It is important to point out that both eqn.[0.3] and [0.4] also imply that AA are interchangeable from a regulatory standpoint, and a single AA at very high concentrations can substitute for low concentrations of other AA.This could lead to mass balance problems; however, the substantial regulation of mammary AA transport mitigates the potential for such a problem (Arriola Apelo et al., 2014b), and based on the observed efficiencies of conversion of individual EAA to milk protein predicted from the system (see below), this does not seem to occur within in vivo ranges.Arriola Apelo et al. (2014d) used a similar set of quadratics for his work and found that the plateaus occurred at EAA concentrations well above physiological ranges.As nonlinearity with respect to DEInp was not observed when exploring the effects of DEInp and MPIn, it was not considered a contributor to the quadratic term.
Correlations and covariance among the fitted parameters was assessed for the eqn.[0.4], and the results are presented in Table S2.The correlations among parameters were surprisingly low having a maximum in absolute terms of 0.53 with most below absolute values of 0.30, indicating the model was generally uniquely definable from the data.
As before, eqn.The NEAA term was also significant suggesting that, although a small effect, the NEAA play some role in driving milk protein output and thus should not be totally ignored.This is consistent with observations of sestrin-2 signaling in response to EAA and NEAA (Luo et al., 2018), which would impact mTOR phosphorylation and protein synthesis.Several NEAA are also actively transported into cells and subsequently exchanged for extracellular EAA via the LAT1, LAT2, Y + LAT1, and Y + LAT2 transporters (Broer and Broer, 2017).Thus, low blood concentrations of NEAA may reduce EAA transport activity (Hruby Weston et al., under review) leading to the reduction in milk protein output.
Excepting Trp, the individual EAA terms were more consistent across EAA.Using the milk protein composition of NASEM ( 2021 The effect of DE from digested FA was presumably mediated via energy supply to the mammary glands as dietary FA do not generally drive ruminal microbial growth (Hanigan et al., 2021) or insulin and IGF-1 concentrations (Becú-Villalobos et al., 2007), although highly saturated FA have been observed to elicit an insulin response (Harvatine and Allen, 2005).
The intercept was more negative than for eqn.The correlations among the linear EAA terms and the quadratic term ranged from a low of −0.12 for Thr to −0.41 for His.Based on these, one can conclude that the slope of the common term was generally independently defined by the data.
The model exhibited slightly greater parameter stability during cross-evaluations (data not presented) with only the NEAA term exhibiting some uncertainty with a CV of 48.2%.The Trp term was more stable with a CV of 33.7%, and exhibited maximal correlations with Lys and Thr terms at modest values of −0.40 and −0.36, respectively.Thus, this simpler model resulted in slightly more stable slope estimates that were consistent with the known biology excepting the Trp slope.Based on these attributes, the similar fixed effects performance favors eqn.[4.2] over eqn.[4.1].

Interactions Among Nutrients
Based on prior work showing that interactions among nutrients was not generally significant (Vanhatalo et al., 2003a, Raggio et al., 2006, Rius et al., 2010c, Appuhamy et al., 2011b, Appuhamy et al., 2012, Appuhamy et al., 2014a, Arriola Apelo et al., 2014d, Yoder et al., 2020), our screening work was conducted without consideration of potential interactions.To test this assumption, we examined all the potential 2-way interactions among the linear energy and EAA terms in eqn.S3), which contained 20 significant, 2-way interactions.Although the model performed slightly better than eqn.[4.2] and [4.1] [AIC = 10,296, marginal R 2 of 0.54, and an RMSE of 128 g/d (13.9%)], an RMSE gain of 4 g/d hardly justifies the addition of 20 parameters and the associ- The linear terms for absorbed Phe and Thr were very unstable during cross-evaluation (data not presented) having CV of 403 and 59%, respectively.This is an indication that the data were not adequate to define the responses given the variation in responses among trials.In the absence of additional observational data, they should be removed from the model.All other terms had CV less than 36% indicating relative stability during cross-evaluation.
The correlation matrix suggested some collinearity problems with negative correlations among the linear and quadratic absorbed EAA terms ranging as great as −0.98.There were also moderate correlations of 0.50 for DEStIn and DENDFIn, −0.48 for Phe and the quadratic Trp term, and 0.47 for absorbed Phe and Trp percentages.The remainder were below absolute values of 0.43 with most less than 0.30.Removal of the Arg, Phe, Thr, and Trp quadratic terms addressed the collinearity problems, but the AIC increased to 10,436.
The intercept was −3,687 g/d, which greatly exceeded the mean MPIn and even CPIn, and thus was well outside of the biologically possible range.The model is clearly centered in the lactating data used to fit the model, but it must transition from the center of the data to a biologically impossible solution as the inputs approach 0, and thus predictions from low inputs would be highly suspect.When the quadratic terms were reduced to address collinearity concerns, the intercept assumed a value of 77 g/d, which was greater than expected, but greatly improved over eqn.[4.2].
Estimates for the effects of the DE components were essentially the same as for the other models, and the linear slope for total EAA was 0.459 g milk protein/g EAA which was less than for eqn.[3.8] indicating the individual quadratic terms for the EAA proportions did not exert as much influence as the summative term in eqn.[3.8].There were positive linear effects for the percentages of Arg, His, Lys, Met, Phe, Thr, and Trp; and negative quadratic effects for Arg, His, Lys, Met, and Trp.The plateaus for each EAA were: Arg = 10.6%,His = 5.0%, Lys = 17.0%,Met = 6.3%, and Trp = 2.2%.With the EAA representing an average of 44.8% of the total AA supply, these can be converted to the equivalent of percentages of MP yielding: Arg = 4.77%, His = 2.25%, Lys = 7.63%, Met = 2.84%, and Trp = 0.98%, respectively.The Lys and Met zeniths compare favorably with the estimates derived by NRC (2001).
Excepting the linear Phe and Thr terms and the collinearity concerns with the quadratic terms, this model was fairly well defined by the data and consistent with prior requirement models, but the form of expression of the individual EAA responses is problematic.Because each term is divided by total EAA, a change in the absorbed supply of any of the individual EAA not only causes a change in the numerator for that EAA, it also causes a change in the denominator for all other EAA.This results in nonlinearity in the relationships among the ratios and the independent variable (a violation of the linear regression assumptions), and heterogeneity in individual EAA slopes with the contribution of each EAA to milk protein varying depending on the summative values of the other EAA.Thus, the effect of the addition of 1 g of absorbed Met will be different when added alone than when added with another EAA or as a component of a protein.In the latter case, the effect will be quite different given the much larger change in the denominator.Such a response is not consistent with existing knowledge where each factor primarily exerts an independent effect (Rius et al., 2010c, Appuhamy et al., 2011b, Appuhamy et al., 2012, Appuhamy et al., 2014a, Arriola Apelo et al., 2014d, Yoder et al., 2020).However, the ability of the mammary glands to alter individual AA transport activity in response to AA supplies (Bequette et al., 2000, Hanigan et al., 2000), and the flexing of liver clearance in the opposite direction yields variable efficiency of conversion of AA supply to milk protein (Fleming et al., 2019a), which is perhaps the source of the interactions identified in Table S3.Such flexing across the tissue beds, which generally reduces mammary affinity and increases splanchnic affinity as absorbed supplies increase, is directionally consistent with the denominator effects within the ra-  4), but the latter are inconsistent with the 9 positive, 2-way interactions among AA identified in eqn.[S3.1], and the ratios are highly unlikely to replicate in vivo response surface given the ratios are fixed mathematically.Finally, the positive interactions identified in Table S3 are directionally and quantitatively inconsistent with the negative effects generated by the encoded ratios.For these several reasons, eqn.[4.3] cannot be recommended for general field use, and was not considered further.
Absorbed EAA Expressed as a Ratio to Nonprotein Digested Energy Intakes.The final equation form that was considered expressed the effect of individual absorbed EAA as a ratio to DEInp (eqn.[0.6]).This form is used in swine models (NRC, 2012).It suffers from the same mathematical deficiencies and violation of linearity assumptions as eqn.[0.5] due to the use of the ratios.Backward elimination yielded eqn.
[4.4], which had 16 fixed terms.There were positive linear terms for the DE components, Arg, Lys, Phe, and Trp and a negative linear term for Thr.It had negative quadratic effects for Arg, Phe, and Trp, and positive quadratic effects for Met and Thr.The AIC was 10,412, which was the same as eqn.[4.2] and greater than the other equations.The marginal R 2 was 0.51, and the RMSE was 129 g/d with no meaningful mean or slope bias, and a relative prediction error of 14% of the observed mean.
Under cross-evaluation, the linear and quadratic Trp terms and the quadratic Arg term were unstable with CV of 96, 49, and 49%, respectively.The quadratic terms for Arg, Phe, Thr, and Trp were highly correlated with the respective linear terms at −0.99 for all indicating inadequate data range to define those terms.The Thr quadratic term was also correlated with the Trp and Phe quadratic terms at −0.61 and −0.60, respectively.There were also correlations among the linear terms for Phe, Thr, Arg, and Trp all exceeding −0.58.There were cross-correlations among the linear EAA terms and other quadratic terms ranging from 0.55 to 0.60.Removal of the quadratics for Arg, Phe, Thr, and Trp resolved all of the correlation problems, but the AIC increased to 10,443 indicating the model was clearly inferior to the other forms.Based on inferior performance, the wide range in collinearity among the parameters, parameter stability problems during cross-evaluation, and the use of ratios in the model, the model cannot be recommended for field application, and was not considered further.

Model Selection by Evaluation of All Possible Models
The screening work established that eqn.[0.4] was generally biologically consistent, excepting the negative linear Trp slope, and the data were adequate to define the energy terms and at least 5 EAA.Because backward and forward elimination algorithms can result in biased model selection (Heinze and Dunkler, 2017), we evaluated all possible combinations of the fixed terms in eqn.[0.4].A summary of the top 5,000 models based on AIC from a total population of 524,288 is provided in Table S4.DEFAIn, DErOMIn, and DEStIn were in all of the top 5,000 models and DENDFIn was in 4,554 of the 5,000 underscoring the importance of DE as a driving factor for milk protein production.Absorbed Lys, Met, and the quadratic EAA term were in 4,535, 4,718, and 4,057 of the 5,000 indicating their relative importance in explaining variation in milk protein production.His and Ile were in more than 3,200 of the models, and the remaining AA including the aggregated NEAA were in roughly half of the top models.Parity and BW were in 2,271 and 3,180 of the top models, respectively.From the CV for the parameter estimates for the 5,000 models, it can be observed that several of the potential terms had extremely variable slope estimates, indicating they were heavily dependent on the presence or absence of small subsets of the data, and thus not well defined.These included absorbed Arg, Leu, Phe, and Val.The remaining terms were more consistent across models indicating they were well defined by the data, and should represent unbiased estimates of the true effects.
Table 5 contains a summary of a selection from the top models including the global model (eqn.[5.1]) and the weighted average model (eqn. [5.2]).The best model (eqn.[5.3]) contained linear terms for absorbed His, Ile, Lys, Met, Thr, Trp, and NEAA; all 4 DE terms; and BW.It had an AIC of 10,389 which was 8 units less than the global model.The difference in AIC is due to the penalty applied to the number of terms in the model.Although this model was the best based on AIC, it lacked some biological coherence given the negative slope for absorbed Trp.The global and weighted average models also had negative terms for Trp.
The first model that did not include a negative linear slope for Trp was rank order 167 (eqn.[5.4]).It had an AIC of 10,397, a marginal R 2 of 0.49, and RMSE of 132 g/d, a relative error of 14.3%, and negligible mean and slope bias.In this regard, it performed almost identically to the top and global models, and had the advantage of biological consistency.The intercept was numerically positive, but not significantly different from 0. The absence of the relatively large negative Trp term resulted in a general reduction in the slope estimates for the other EAA except for Met.Given the low correlations among parameter estimates, the changes should not be due to co-linearity problems.It is more likely the general reductions in slope estimates simply reflect the removal of the significant negative driving term, which algebraically requires an offsetting reduction in one or more other parameters (all positive) to allow the predictions to remain centered in the data.Although the negative slope for Trp is seemingly biologically inconsistent, it clearly is capturing some aspect of the observed variance.
The derived slopes of 0.96 g milk protein/g absorbed Lys and 1.92 g milk protein/g absorbed Met were well below the ranges of 16 to 4 g milk protein/g Met intake and 5.0 to 3.2 g milk protein / g Lys intake reported by Vyas and Erdman (2009).However, Vyas and Erdman (2009) only controlled independently for Lys and Met intake.Any DMI or DE concentration changes that occurred would have altered intakes of DE and the other EAA, and the effects of those would have been assigned to Lys and Met.Lean et al. (2018) reported coefficients of 1.04 g milk protein / g metabolizable Lys and 1.89 g milk protein / g metabolizable Met, which are essentially equivalent to values derived using eqn.[5.4] and other models where negative linear EAA terms were excluded.Plots of residual errors for milk protein production also support these lower estimates given essentially no evidence of bias across the entire data set (Figure 3 and Figure 4), studies using infusions (Figure S2), or studies feeding rumen-protected Lys and Met (Figure S3) with respect to the supply of any of the EAA.The lack of a slope for the residuals within experiments indicates the model correctly captured the response to each of the EAA.Lean et al. (2018) also found that the NEAA were a significant determinant of milk protein production having a slope of 0.23 g milk protein / g metabolizable NEAA, which is almost double the estimate herein.A portion of the difference is explained by the method of estimating NEAA used by Lean et al. (2018).The assumed equivalence of AA mass in the free and protein bound forms, when in fact, each peptide bond is formed by exclusion of a molecule of water results in approximately a 15% loss of mass when bound in protein form.Estimating NEAA as MPIn -EAA as in Lean et al. (2018) results in approximately a 30% underestimate of the NEAA supply which would propagate to the slope estimate.All AA mass interconversions between peptide bound and free forms were appropriately adjusted herein.Such an adjustment of the Lean et al. (2018) estimate yields a slope of 0.161 g/g which is not different than the 0.114 g/g estimate for eqn.[5.4].
Other experimental work (Metcalf et al., 1996, Doepel andLapierre, 2010) does not support the observation of NEAA effects.Doepel and Lapierre (2010) infused 356 g NEAA/d into the abomasum with and without infused EAA and Metcalf et al. (1996) infused 192 g of NEAA into the jugular vein.Neither reported significant changes in milk protein output although Doepel and Lapierre (2010) did observe a numerical effect when the infusion was combined with additional EAA.The expected change in milk protein output based on eqn.
[5.4] was 40 g/d for Doepel and Lapierre (2010) and 22 g/d for Metcalf et al. (1996).Given an SEM of 37 g/d for the first and an SED of 14 g/d for the second, the expected response may not have been detectable.
Leucine is often considered the most potent stimulator of mTOR activity, and thus an important driver of protein synthesis (Lynch, 2001, Escobar et al., 2005).Rank order 172 (eqn.[5.5]) includes an absorbed Leu term, although the slope estimate is relatively low at 0.37 g/g with a SE of 0.26 g/g.The relative error is large, but it should be noted that the SE is half of that for absorbed Ile and a third of that for His and Thr, and removal of Leu impacts the estimates of His, Ile and NEAA.Contrary to the conclusion reached by Lean et al. (2018) that Leu was a potent driver of milk protein synthesis, it seems more likely that Ile is exerting more influence on milk protein production than Leu under normal feeding conditions.
Differential responses to systemically infused EAA as compared with those absorbed from the intestine are likely real, but the magnitude is small enough that it is difficult to discern the effect.The InfTP Art term was only present in half of the top 5,000 models fitted to all of the data, and the slope estimate averaged 0.173 g of milk protein / g of infused EAA indicating a 17.3% first-pass loss of the EAA during absorption.The SE was as large as the estimate indicating the estimate does not significantly differ from 0. From prior work (MacRae et al., 1997, Hanigan et al., 2004, Fleming et al., 2019a), the utilization of EAA on first pass varies among EAA, but generally being 15% or less which is consistent with the estimate herein.Given the variance, much more data would be required to estimate the loss with any degree of precision using milk protein production as a response variable.Given that the presence or absence of the term in the model did not affect the remaining parameter estimates (eqn.[5.5] vs [5.6]), it can be excluded from the model for field application.
Equation [5.7] was rank order 277 and represents the first model that did not include a term for Thr.As compared with eqn.[5.5], its removal did not significantly alter the slope coefficients for Lys and Met, but it did increase the coefficients for His and Ile without greatly altering their SE.The removal also resulted in a  Intercept estimates ranged from −166 to 257 g/d across the 5,000 top models with a mean estimate of 44 g/d (Table S4).However, screening and selection based on fit statistics and biology resulted in models that generally had slightly positive intercepts (Table 5) as compared with an expected value in the −150 to −200 g/d range as discussed above.This may reflect differences among the AA composition of milk protein (the response element defining both the intercept and  the slope) and the composition of proteins representing the non-productive functions.However, if this were the case, one would expect the global model, which contained all of the AA represented in MP, to have a similar intercept as for eqn.[3.5].It is also possible that the common quadratic term failed to fully capture the curvature in several of the EAA, but the intercept for eqn.[4.1] was only slightly negative with most of the individual AA quadratics represented.Given that there were multiple models with the expected intercept among the population of models (Table S4), it seems the problem is related to effects captured in AA that were screened which yielded slopes for the remaining AA that deviate slightly from the expected intercept.Given the apparent discontinuity between the intercepts for equations in Table 5 and factorial estimates (NASEM, 2021), predictions of milk protein below the range of data used herein to derive the equations should be evaluated for validity before use.
[5.6] are presented in Table 6.Stability of parameter estimates generally aligned with the SE of the estimates in the original solves with the Leu, Thr, NEAA, and InfT-PArt terms exhibiting large variation in estimates (CV >50%) across the solves, and mean P-values indicating they were not significantly different from 0. Terms for absorbed His and Ile had moderate CV, but were generally significantly different from 0, and the remainder excepting the intercept were generally stable across solves.Despite moderate variance in the estimates of model parameters, the resulting model evaluations were quite stable.Based on the marginal and conditional R 2 , approximately 2/3 of the remaining variance was associated with random study effects which would generally represent variation among studies in animal genetics, facilities and facility management, environmental and disease stress, nutritional factors not captured as inputs to the model, and random errors of measurement (see Ch. Twenty in NASEM, 2021).Thus, some of this variance could be explained through collection of a wider array of data, but some would remain as truly random, and not biologically explainable.
One challenge with several the EAA was that relatively few studies have been conducted with dietary AA inputs independently manipulated by infusion or with a well-defined rumen-protected AA source.Therefore, the independent variation in inputs of those EAA was not well represented in the data.In selecting 85% of the data for training, several key treatments defining one of those EAA likely would have been excluded from the training set in one or more iterations resulting in poor estimates of those effects.Thus, additional experimental work is needed to gain precision and better evaluate the independent effects of several AA including Arg, Leu, Phe, Thr, Trp, and Val so that confidence in those effect estimates can be improved.
Exploration of Residual Errors.Residual errors were examined to determine whether any additional factors should have been considered in the models.From Figure 3, the residuals offer no evidence of slope bias and are centered around 0. There also are no apparent patterns to the residuals within experiment although there are differences among experiments in the intercept which were adequately handled by the random effect of study as shown in the lower panel of Figure 3.There were also no clear patterns to the residuals when plotted against total EAA divided by DEInp indicating that the independent effects of these 2 terms were properly captured in the model.There were no obvious patterns to the residuals across or within studies when plotted against dietary concentrations of dFA, dNDF, dSt, or rOM or against BW.There was slight evidence for a potential negative slope to the residuals within study relative to increasing milk fat concentration, perhaps indicating that increasing energy use for milk fat reduced the yield of milk protein per unit of input (data not shown).

Consolidation of Eqn. [5.6].
Because the coefficients for DEStIn, DErOMIn and DEFAIn were similar, and an infusion term was not needed for field use, a simplification of eqn.[5.6] was undertaken.In addition to aggregation of the DE terms, the addition of the unused EAA to the NEAA term was explored, but contrary to the earlier NASEM (2021) work, we found that it provided no additional statistical benefit (data not shown), and thus we did not include an Other AA term in place of the NEAA term.The simplified model, eqn.[7.2] is provided in Table 7 along with an evaluation of the published NASEM (2021) equation and the NASEM (2021) equation refitted to the data.The NASEM (2021) equation was derived earlier in the analytical process, and thus represents our exploration of the potential differential effects of energy coming from differing substrates using nutrient concentrations in the diet.We later discovered this could be more cleanly addressed using the DE terms for each which also negated the remaining collinearity problems for the DEInp term.The lesser slope for DENDFIn than for other DE terms is more easily interpreted in this form than when using an aggregated DEInp term and adjustment slopes based on dietary NDF concentrations.The NASEM (2021) equation was also derived before more data screening uncovered additional entry errors that were corrected, hence the refit of the prior model form to the current data to allow like comparisons.Based on RMSE, the refitted NASEM (2021) equation performed the same as the original equation; however, the slope Aggregation of terms yielding eqn.[7.2] increased the AIC very slightly to 10,405, but the RMSE only increased 1 g/d.Thus the consolidated model performed as well as eqn.[5.6], and as well as the equation provided in NASEM (2021).Not only do these models perform equally well, but there are many additional equations that perform similarly.Some were ruled out for biological reasons, but many additional ones are potentially valid.
Given the superiority of the equations based on independent and additive responses to the individual nutrients as compared with those based on first limiting nutrient principles (Table 4) and the extensive mechanistic work (Cant et al., 2003, Anthony et al., 2004, Burgos and Cant, 2010, Appuhamy et al., 2011b, Appuhamy et al., 2012, Appuhamy et al., 2014a, Arriola Apelo et al., 2014c, Arriola Apelo et al., 2014d, Doelman et al., 2015a, Doelman et al., 2015b, Liu et al., 2017, Nichols et al., 2017, Yoder et al., 2020), the first limiting nutrient hypothesis must be dismissed.The disadvantage of the revised representation of nutrient effects is that it is not possible to define requirements for any of the terms in a classical sense.Previous mod-els (NRC, 1966(NRC, , 1971(NRC, , 1978(NRC, , 1989(NRC, , 2001) ) were based on assumptions of fixed efficiencies of nutrient use with no additivity of responses, which is clearly not consistent with the observed data (Figure S4).Under those assumptions, unique requirements for each nutrient could be algebraically determined.Given observed independent, additive, and nearly linear responses throughout the normal feeding range, there are an infinite number of solutions that provide the same level of production, and the best diet will be the one that provides an optimum mix of nutrients that achieves the desired level of performance at a minimum cost.
Although much of the variance is captured in these equations, the random effect of study still accounts for approximately 35% of the overall variation (Table 5).Some of this variation is truly random errors of measurement, but a portion will be due to environmental effects, health challenges, and genetic influences (NASEM, 2021).As these cannot currently be captured in the predictions a priori, the model may exhibit mean bias when applied to individual farms and pens within farms, but the predictions should be representative across a range of farms and pens.The plateau defined by the quadratic term in the equations is also problematic for farm application as it represents an apparent  maximum milk protein production for animals used in studies over the past 30 years which is generally well below that achieved on well managed farms today.As genetics and management at a particular farm change slowly over time, provision should be made for intercept and plateau adjustments by farm as outlined by NASEM (2021) to accommodate farm specific effects until predictions of those effects can be included in the system.

CONCLUSIONS
The existing data were adequate to define the effects of DE, Ile, Lys, and Met on milk protein production with good precision, the effects of His and the NEAA with moderate precision; and the effects of Leu and Thr with poorer precision.Milk protein production equations based on additive and independent effects of the individual AA and energy providing nutrients performed much better than equations based on 1st limiting nutrient concepts.Milk protein responses to supplies of MP and total EAA were curvilinear yielding variable efficiencies of conversion of absorbed EAA to milk EAA.Linear responses to individual EAA were observed, but quadratic responses were not well defined by the data.Expression of EAA as percentages of MP or as a ratio to DEIn resulted in slightly poorer model performance, and given such ratios also violate linearity assumptions relied upon to derive the regressions, those forms cannot be recommended.Several potential interactions among AA and energy providing nutrients were identified, but correlations among the parameters were high indicating some may be artifacts of collinearity problems.Further progress will require completion of additional experiments that isolate the individual effects of Arg, Leu, Phe, Thr, Trp, and Val; explore the effects of individual NEAA; and evaluate potential interactions among the AA and energy providing nutrients.Because many of these nutrients also affect DMI and correlations among nutrients, further modeling work must consider all nutrient inputs to the animal as undertaken herein, or the data must be screened to remove all treatments affecting intake of any of the other macro-nutrients that are not being directly studied to avoid biased estimates of the nutrients of interest.e Derived before additional data cleaning and not re-calculated.
Hanigan et al.: Predicting Milk Protein independent.The complete data set contained 1153 treatment means from 273 experiments.
1] atic for separating the effects of each.Similar problems are introduced for partitioning DEIn effects based on substrate type.Therefore, DEIn was parsed to the individual contributions from starch, NDF, TP, FA, and rOM (DEStIn, DENDFIn, DETPIn, DEFAIn, and DErOMIn, respectively; Mcal/d), and the latter used in regressions.Combinations of these were created ad hoc except for non-protein DEIn (DEInp, Mcal/d) used by (NASEM, 2021) which was precalculated for use in regressions.
Eqn. [3.1] containing linear and quadratic terms for CPIn and a linear DEIn term performed nearly as well, based on AIC, as eqn.[3.4] using linear and quadratic MPIn terms, the latter having a gain in precision of 22 g/d.The total EAA based equation performed the same as the MP-based equation reflecting the very high correlation between MPIn and EAA intakes (Figure S1).Digestible energy intake was a very strong linear, additive determinant of milk protein production regardless of whether the DE from MP was included in the energy term (DEIn; eqn.[3.4]) or excluded to reduce collinearity (DEInp; eqn.[3.5]).Although the latter resulted in a very slight decrease in statistical performance, removal of the collinearity problem likely revealed the true independent effects of DEInp and MPIn.Considering that 22% of DEIn derives from MPIn, an estimate of 10.99 g protein/MCal of DEIn (eqn.[3.4]) equates to a partial slope for MPIn of 0.428 g protein/g MPIn (10.99 × 0.22 / 5.65).Adding the observed linear slope coefficient for MPIn to the partial energy slope results in an overall response to changes in MPIn of 0.755 g protein/g MPIn (0.428 + 0.327) as compared with the derived slope of 0.409 g milk protein/g MPIn from eqn. [3.5].A slope estimate eof 0.755 g/g is more than that used by NRC (2001), which resulted in a very large slope bias.Therefore, one can conclude the slopes from eqn. [3.4] are non-biological presumably due to collinearity problems.There was no evidence of an interaction between DEInp and protein terms (discussed further below).At a mean DEInp of 48.0 Mcal/d, the DEInp term in eqns.[3.5] and [3.8] represented 530 g/d and 532 g/d of milk protein, respectively versus 794 g and 812 g/d, respectively associated with the CPIn and MPIn terms.Thus, if representative, approximately 40% of the driving force for milk protein production is from the nonprotein energy supply.
[3.3]   yielded an intercept that did not differ from 0. When DEIn was also included in the regression, the intercept assumed a value of −150 g/d (eqn.[3.4]).Use of DEInp to remove the effects of MP from the energy term yielded an intercept of −165 g/d NP/d (eqn.[3.5]).This was less than estimates of 223 g/d predicted by NASEM (2021) for scurf and endogenous urinary NP losses by a fasted, non-productive, 650 kg animal, but the difference was less than the SE of the intercept estimate.From eqn.[3.3], the efficiency of MP conversion to NP in a fasted state should be approximately 0.65, which is consistent with NRC (2001) estimates.Because the slope presents the fractional capture of MP in milk NP,Hanigan et al.: Predicting Milk Protein

Figure 2 .
Figure 2. Observed (black) and residual errors (red) for milk protein production versus predicted based on NRC (2001) predictions of NEL allowable, MP allowable, the minimum of NEL and MP allowable, and the minimum allowable predictions from NEL, MP, and each of the 10 EAA.The EAA allowable predictions were calculated using the same framework as for MP but substituting each of the individual EAA to calculate supply, requirements, and thus allowable milk protein production.
= mean random effect variance, τ00 = between group variance, AIC = Akaike's Information Criterion, BW = body weight centered to the mean value of 612 kg, CCC = concordance correlation coefficient, Conditional R 2 represents fixed and random effects variance, DEIn = digested energy intake, DEInp = DEI minus DE from MPIn, EAA = essential AA intake, ICC = interclass correlation coefficient, Marginal R 2 represents the fixed effects variance, MPIn = MP intake, Obs = observed; Pred = predicted; RMSE = root mean squared error (g/d or % of Obs mean); Mean and Slope bias % expressed as a % of the MSE.1 minus the slope reflects the variable cost of production including metabolic fecal protein and inefficiencies of protein capture resulting in protein catabolism.At an efficiency of 0.65, the MP intake required to meet endogenous urinary and scurf protein use would be 254 g/d.The difference of 89 g/d is less than the 127 g NP/d requirement for metabolic fecal protein predicted by NASEM (2021) for an animal at maintenance intake levels, but as DMI increases with increasing milk production, the estimated non-productive NP use does align with the expected losses based on the inefficiency of conversion of MP to milk NP.Thus, the factorial and regression approaches do not perfectly align with respect to non-productive protein requirements at maintenance intakes and likely at very low levels of production.Possible reasons for the misalignment include extrapolation errors given a minimum milk protein of approximately 300 g/d in the data set, and use of a milk protein efficiency for estimation of non-productive functions.If the independent effects of DEInp and MPIn in eqn.[3.5] are accepted, the efficiency of conversion of MP to NP increases as DEInp increases (NASEM [3.6] is 209 g NP/d for a fasting, non-pregnant, nonlactating cow as compared with an estimate of 241 g/d for the same cow by the NASEM (2021) model.Eqn.[3.6] predicts a requirement for a 520 kg Jersey in the same state as 146 g/d, and the NASEM (2021) model as 180 g/d.The differences between the high and low BW cows calculated by the 2 models are almost identical at 32 and 34 g/d, respectively, indicating the BW term of eqn.[3.6] was justified and appropriately represented.
[0.4] was reduced by backward elimination resulting in eqn.[4.2].The reduced equation had question of whether the data are adequate to uniquely define Trp effects.The quadratic term, applied to the sum of the individual squares of each EAA, was highly significant, further supporting the need for some form of nonlinearity in the model.From the quadratics, absorbed EAA inputs (g/d) at response plateaus were: His = 710, Ile = 512, Lys = 385, Met = 503, Thr = 769, and Trp = −1997.These are 1480, 440, 250, 1090, 590, and −7680% of the mean inputs of each (Table 1), respectively.Comparison of the casein synthesis plateaus observed by Arriola Apelo et al. (2014d) with mean plasma AA concentrations reported by Swanepoel et al. (2016) yields approximate plateaus relative to mean concentration of 300% for Ile, 200% for Leu, 800% for Met, and 300% for Thr.Thus both observations suggest the plateaus for individual EAA are far above observed plasma concentrations yielding mostly linear responses in the normal in vivo range which explains the difficulty in deriving individual quadratic terms for each EAA.
), the gross efficiencies of conversion of absorbed EAA to milk protein EAA based solely on the linear terms were: His = 7.0%, Ile = 10.7%,Lys = 11.5%,Met = 5.2%, and Thr = 12.0%.Unlike the EAA terms, the DE terms in both equation forms were fairly consistent with the DE from starch, rOM, and FA all having similar coefficients of approximately 10 g milk protein/Mcal of DE, and the DE from NDF having a coefficient of 5 g milk protein/ Mcal.The lesser effect of DENDFIn likely reflects the lower potentiation of insulin and IGF-1 secretion as compared with digested starch and rOM (Garnsworthy et al., 2009), and it is consistent with the negative digested NDF term in the milk protein equation of NASEM (2021).The NASEM (2021) equation used an aggregated DEInp term for energy where starch, rOM, FA, and NDF were given equal weight.As digested NDF apparently has half the effect, the regression solved for a negative term for digested NDF concentrations to offset a portion of the digested NDF response within the DEInp term.
[4.1] and similar to that derived for MPIn (eqn.[3.5]), although less than expected based on factorial summations of NP requirements summarized by NASEM (2021).If we assume the factorial estimates are unbiased, an intercept of −138 g NP/d may indicate a slight potential for slope bias in the milk protein predictions at very low EAA inputs.
[0.4].Reduction of the modified equation by backward elimination resulted in eqn.[S3.1] (Table Hanigan et al.: Predicting Milk Protein ated possibilities for co-linearity problems and spurious estimates driven by outliers.Additional discussion is provided in the Supplement (https: / / doi .org/ 10 .3168/jds .20XX-XXXXX), and no further consideration was given to the approach.Absorbed EAA as a Proportion of Total EAA.Eqn.[0.5] represents an equation form with linear terms for DE and total EAA supplies, and linear and quadratic terms for the percentages of total EAA represented by each EAA.The subcomponents of DEInp and EAA supply captures the independent effects of energy and protein and ensures mass balance.The remaining AA terms reflect the effects of the EAA composition of the total EAA supply on milk protein production, and as for eqn.[0.3], independent plateaus are implied by the individual quadratics for each EAA.As before, the model was reduced by backward elimination resulting in eqn.[4.3], which performed similar to eqn.[4.1] and slightly better than eqn.[4.2] based on an AIC of 10,363.The model had 19 fixed terms as compared with 14 terms for eqn.[4.2].The RMSE for the fixed components was 129 g/d which was marginally better than eqn.[4.1] and [4.2] by 3 g/d.The relative prediction error was 14.0% of the observed mean, with no mean bias and minimal slope bias.
Hanigan et al.:  Predicting Milk Protein tios of eqn.[0.5] and its reduced forms (Table 's information criterion, BW = body weight (kg); Conditional R 2 represents fixed and random effects variance, DEFAIn = digested energy intake from fatty acids (Mcal/d), DENDFIn = digested energy intake from NDF (Mcal/d), DErOMIn = digested energy intake from residual OMs (Mcal/d), DEStIn = digested energy intake from starch (Mcal/d), InfTP Art = the effect of systemic infusion of AA versus intestinal absorption, ICC = interclass correlation coefficient, Marginal R 2 considers only the fixed effects variance,, RMSE = root mean squared error (g/d or % of observed mean).σ 2 = mean random effect variance, τ00 = between group variance.1 unit increase in AIC, and a 1 g/d increase in RMSE.The changes in the His and Ile coefficients may indicate a lack of orthogonality, but the correlations among the slopes for each versus Thr are low suggesting that is not the case.It is unclear which estimates for His and Ile are more representative, but given the low correlations among the EAA parameter estimates, it is likely to be those from the more complicated model.

Figure 3 .
Figure 3. Residual errors for predictions of milk protein by eqn.[5.6].Data in the top panel were uncorrected for the random effects of study, and data in the bottom panel were corrected.

Figure 4 .
Figure 4. Milk protein residual errors for eqn.[5.6] versus predicted absorbed supplies of EAA.Predictions and residuals were corrected for the random effect of study.
Hanigan et al.: Predicting Milk Protein   estimates for absorbed His and Leu declined and the estimate for absorbed Ile increased.
energy intake from fatty acids (Mcal/d), DENDFIn = digested energy intake from NDF (Mcal/d), DErOMIn = digested energy intake from residual OM (Mcal/d), DEStIn = digested energy intake from starch (Mcal/d), InfTP Art = the effect of systemic infusion of AA versus intestinal absorption, Obs = observed, Pred = predicted, Mean and Slope bias % expressed as a % of the MSE, RMSE = root mean squared error (g/d or % of Obs mean).

a
DEInp = digested energy intake minus the DE associated with absorbed protein (Mcal/d), DEFAIn = digested energy intake from fatty acids (Mcal/d), DENDFIn = digested energy intake from NDF (Mcal/d), DErOMIn = digested energy intake from residual OM (Mcal/d), DEStIn = digested energy intake from starch (Mcal/d), DigNDF = digested NDF (% of DM), dSt = digested starch (% of DM), RMSE = root mean squared error (g/d or % of Obs mean).b NASEM 2021 model was derived using a prior version of the data set before additional cleaning, thus the parameter estimates reflect the earlier data, and the AIC and random effect statistics are not comparable to those derived from the current data.The RMSE analyses reported here reflect performance of that equation using the revised data set herein.c summation of other EAA not explicitly represented plus the NEAA.d summation of the squared individual, absorbed EAA explicitly represented in each equation, e.g., His, Ile, Leu, Lys, Met, and Thr (g 2 /d).
Hanigan et al.:Predicting Milk Protein intake per day, the ratios of individual EAA to DEInp algebraically also reflect mass quantities of absorbed EAA for each negating the need to specify total EAA in the model to maintain conservation of mass: of R. Equation numbers are of the form (Table. Equation) to aid in finding equations within the manuscript.The number to the left of the decimal point denotes the table num-ber and to the right the equation within the table.A table number of 0 denotes equations presented in the body of the text.

Table 1 .
Descriptive statistics of the data set used to derive model parameters and test the models.n = 905 a DEIn = digested energy intake, DEInp = DEIn minus DEItp, DEItp = DEIn from digested true protein.

Table 2 .
Hanigan et al.: Predicting Milk Protein   Residual error analyses of NRC (2001)predictions of NEL allowable, MP allowable, the minimum of NEL and MP allowable, and the minimum of NEL, MP, and each of the 10 EAA allowable predictions.The EAA allowable predictions were calculated using the same framework as for MP but substituting each of the individual EAA to calculate supply, requirements, and thus allowable milk protein production a CCC = concordance correlation coefficient, RMSE = root mean squared error (g/d or % of Obs mean), n = number of treatment means, RMSE = root mean squared error (g/d or % of Obs mean.

Table 4 .
Hanigan et al.:Predicting Milk Protein Evaluation of the form of expression of the individual EAA.Models were derived by backward elimination from the global models (eqn.[0.3] to [0.6] until all terms were significant at P < 0.05.Data (n = 905) were weighted by the square root of the number of experimental units per treatment and Experiment ID was included as a random variable

Table 5 .
Hanigan et al.: Predicting Milk Protein Selected results from evaluation of all possible models derived from eqn. [0.4].Rank order was based on AIC.Data (n = 905) were weighted by the square root of N and Publication ID was included as a random variable

Table 6 .
Hanigan et al.: Predicting Milk ProteinResults from cross evaluation of eqn.[5.6] using 85% of the data for training and 15% for independent evaluations.Data for each were randomly selected from the full data set.Selection and evaluations were repeated 500 times

Table 7 .
Evaluation of the consolidation of terms in eqn.[5.6] and a comparison to the equation used by NASEM, 2021