An Integrative Biology Approach to Quantify the Biodistribution of Azidohomoalanine In Vivo

Identification and quantitation of newly synthesized proteins (NSPs) are critical to understanding protein dynamics in development and disease. Probing the nascent proteome can be achieved using non-canonical amino acids (ncAAs) to selectively label the NSPs utilizing endogenous translation machinery, which can then be quantitated with mass spectrometry. We have previously demonstrated that labeling the in vivo murine proteome is feasible via injection of azidohomoalanine (Aha), an ncAA and methionine (Met) analog, without the need for Met depletion. Aha labeling can address biological questions wherein temporal protein dynamics are significant. However, accessing this temporal resolution requires a more complete understanding of Aha distribution kinetics in tissues. To address these gaps, we created a deterministic, compartmental model of the kinetic transport and incorporation of Aha in mice. Model results demonstrate the ability to predict Aha distribution and protein labeling in a variety of tissues and dosing paradigms. To establish the suitability of the method for in vivo studies, we investigated the impact of Aha administration on normal physiology by analyzing plasma and liver metabolomes following various Aha dosing regimens. We show that Aha administration induces minimal metabolic alterations in mice. Our results demonstrate that we can reproducibly predict protein labeling and that the administration of this analog does not significantly alter in vivo physiology over the course of our experimental study. We expect this model to be a useful tool to guide future experiments utilizing this technique to study proteomic responses to stimuli.


Introduction
The use of non-canonical amino acid (ncAA) labeling for selective identification of newly synthesized proteins (NSPs) in mammalian cells was first introduced by Dieterich et al. 1 and has since been applied to study several biological systems (see reviews 2,3 ). In this technique, an ncAA, typically a methionine (Met) analog, is introduced to the biological system of interest and incorporated into newly synthesized polypeptide chains using endogenous or engineered cellular translational machinery. Distinction of nascent proteins from the constituent proteome is enabled by reactive chemical groups, such as azides and alkynes, which can be covalently modified via azidealkyne cycloaddition (a click chemistry reaction). 2 As such, ncAA-labeled NSPs can be selectively conjugated to affinity or fluorescent tags for identification or visualization, respectively. 1,4 This technique has been successfully employed to probe protein dynamics in a variety of bacterial [5][6][7] and mammalian cells in vitro, [8][9][10] as well as model organisms in vivo, including zebrafish 11 and Xenopus. 12 More recently, ncAA labeling has also been shown to be effective in identifying NSPs in rodents. [13][14][15] Unlike similar techniques such as stable isotope labeling, wherein proteins in low abundance are often undetected due to the absence of an enrichment method for selective isolation of newly synthesized labeled proteins from the pool of more abundant unlabeled proteins, ncAAs offer chemical handles that can be used for enriching NSPs. 5 The expanding applications of ncAA labeling will therefore enable previously inaccessible biological questions, wherein understanding the temporal dynamics of protein synthesis and turnover is critical, to be addressed.
Dietary administration of ncAA, typically enhanced with a Met-free diet, has been shown to achieve adequate labeling efficiency in rodents. 13,16 However, as an essential amino acid, Met deprivation may affect normal physiology, particularly over longer labeling periods. Notably, the presence of Met in mammalian diet is critical for normal embryonic development, [17][18][19] which constrains this method to studies of adult animals. Our group has previously demonstrated that labeling the adult and embryonic murine proteomes can instead be achieved via systemic injection of ncAAs without the need for Met depletion. 14,20 Compared to feeding with an ncAA-enriched diet, the injection method achieves global proteome labeling in a shorter period of time, which enables the detection of proteins synthesized shortly after injection and proteins with high turnover rates. 20 Although direct comparisons of Aha administration techniques have not been previously studied, parenteral administration also reduces inherent variability, bypassing hepatic metabolism and eliminating confounding fluctuations in feeding patterns and intestinal absorption. 21 Nevertheless, the injection technique requires technical expertise and may cause distress to animals.
Despite the application of ncAA labeling in a number of studies to decipher complex cellular processes in animal models, 13,15,16,22 understanding of the kinetics of ncAA distribution in tissues, especially as it pertains to rates of protein incorporation and loss by degradation, is lacking. Determination of the timescale of ncAA uptake by tissues following administration and the lag time before maximum protein labeling are critical information for the design of robust temporal experiments to study the nascent proteome. Predicting ncAA kinetics in murine models will also enable optimizing the dosing regimen to attain the ideal concentrations to achieve sufficient protein labeling in the desired tissue over the course of the study.
The aim of this study was to characterize the distribution kinetics of azidohomoalanine (Aha), a widely used Met analog, in mice following subcutaneous injection, and to investigate the impact on normal physiology. To study the biodistribution of Aha, we measured the concentration in the plasma, liver, kidney, brain and hindlimb skeletal muscle using liquid chromatography-tandem mass spectrometry (LC-MS/MS) over a period of 24 h. This dataset was used to develop a deterministic compartment model of small molecule kinetics that characterizes the movement of freely diffusive Aha (fAha) throughout the mouse circulatory system and into tissues. In addition, we used fluorescent western blotting to measure protein labeling in these tissues during the same period. This second dataset was used to inform a model of relative protein labeling (pAha) as a function of fAha availability to characterize both the predicted labeling profile for a given experimental treatment of Aha and the relative synthesis and degradation rates of Aha-labeled proteins. We demonstrate that this model can be used to characterize nascent protein synthesis and turnover within distinct tissues. Furthermore, we validated the capability of this model to predict NSP labeling under more complex, multiple injection dosing paradigms, which can be a tool to guide the design of future experiments utilizing Aha labeling.
In addition to the kinetic analysis, we compared the effect of the various Aha dosing regimens on the metabolome of Aha-injected mice with that of non-injected controls. LC-MS analyses of plasma and liver metabolomes indicated few statistically significant metabolic changes occur in response to Aha labeling. Taken together, these results provide a fundamental understanding of the interrelation between the distribution kinetics of the ncAA into murine tissues and the associated degree of protein labeling, as well as the potential impact of ncAA injection on normal physiology.

Kinetic Model of Aha Biodistribution Describes Transport and Exchange
To capture the dynamics of free Aha (fAha) distribution in vivo, an in silico model system of ordinary differential equations (ODEs) was generated describing the physiological processes of small molecule transport. All model and code files are accessible at the Kinzer-Ursem lab GitHub repository (Additional File 1: S1). Descriptions of the model, equations, and parameter values are also provided in Additional File 1 (S2-S5). To mimic the injection site of our subcutaneous dosing paradigm, fAha was introduced into the model at a non-localized reservoir. From this reservoir, fAha enters the murine circulatory system at a rate that is a function of reservoir concentration and transport kinetics and is distributed to distinct tissue compartments (Fig. 1).
Within each compartment, the time-dependent rate of change of the fAha plasma concentration ([fAha p ]) available for exchange with each tissue can be described as a mass balance with two stages: transport and exchange. The transport stage is governed by circulatory blood flow.
where Q x is the blood flow rate between tissue 'x' and a systemic venous reservoir (sysrv), and V x is the corresponding volume of plasma relevant to each tissue (Q/V represented as a lumped constant qb, Additional File 1: S2, S3). All kinetic parameters for circulatory transport were normalized by tissue mass to compare relative perfusion rates between tissue compartments of differing size. Once localized to a tissue, fAha in the plasma can also be exchanged across the cell membrane with the intracellular Aha concentration ([fAha t ]) and is described as follows: where k i,x and k e,x are the tissue specific import and export rates (min −1 ⋅mg −1 ) for fAha across the cell membrane. The liver tissue, kidney plasma, and gut lumen compartments were assigned additional system removal terms (k r,x , Additional File 1: S2, S3) accounting for excretion and Biodistribution of free Aha (fAha) via transport and exchange. Introduced at a distinct injection site, fAha is allowed to enter the circulation at a systemic venous reservoir. Driven by circulation, fAha is passed through the arterial system (red) into distinct tissue compartments where exchange occurs at tissue specific rates. Arrows indicate directional movement of fAha. Each tissue compartment consists of two sub-compartments: plasma available for surface exchange (red to blue gradient) and an intracellular volume (illustrated here as the bottom compartments: from left to right muscle, liver, brain, and kidney). fAha is incorporated into protein (pAha) in the intracellular volume. A mechanism for oral dosing (via plasma exchange with the gastrointestinal lumen) is illustrated and was included in the model, but not utilized in this study. metabolization of fAha. The transport and exchange equations describing fAha distribution were combined into a single system of ODEs parameterized and bound within reasonable ranges for a model of small molecule pharmacokinetics (Additional File 1: S2-S5). [23][24][25][26][27][28]

Kinetic Model of Protein Labeling Captures Aha Incorporation
Within each tissue compartment, fAha t is incorporated into proteins (pAha) via protein synthesis. As a Met analog, Aha is able to bind to methionyl-tRNA synthase, albeit at a much slower rate (k cat K m −1 Aha: 1.42E-3, Met: 5.47E-1 μM −1 s −1 ). 29 Therefore, at concentrations comparable to normal physiological [Met] for mice fed a methionine adequate diet (MAD) ([free Met] ≈ 4.4 ug g −114 ), Aha binding to the methionyl-tRNA synthase is much slower than Met, yielding Aha incorporation into NSPs of less than 10% by previous estimates. 14 We therefore made a modeling assumption that [fAha t ] is negligibly depleted by incorporation into protein (further justification, Additional File 1: S2) and that recycling of Aha (as a non-canonical analog of an essential amino acid) is also negligible at the scale of plasma clearance.
where k s,x and k d,x are the tissue-specific (denoted for tissue 'x' as above) rate constants of incorporation of Aha into proteins due to synthesis and loss of Aha-labeled proteins due to degradation. These equations describing Aha labeling of proteins in each tissue were added to the biodistribution model establishing a time resolved predictive model of tissue-specific protein labeling given a variety of input dosing paradigms.

Experimental LC-MS/MS and Western Blotting Data Enable Parameter Fitting
Model parameters were initialized and bound within reasonable ranges, informed from literature and experimental measurements as described in the methods (Additional File 1: S3-S5), 23-28 then underwent least squares regression to match experimentally measured data of Aha concentration and labeling. To inform fitting, fAha concentration profiles in plasma, obtained from blood collected by cardiac puncture, and liver, kidney, brain and hindlimb skeletal muscle tissues were determined by injecting Aha subcutaneously into mice at 0.1 mg g −1 total body weight and sacrificing 0.5-24 h post injection (hpi). Accurate identification of fAha in each tissue was performed using LC-MS/MS multiple reaction monitoring (Figs. 2a, 2b). Additionally, the kinetics of Aha incorporation into tissue proteins were described by examining the degree of protein labeling within each tissue over the duration of the study. To this end, tissue homogenates were reacted with biotin-alkyne via copper-catalyzed click reaction, analyzed by semi-quantitative western blotting using a fluorescent streptavidin conjugate and the change in fluorescence intensity relative to the noninjected controls was measured as an analog for pAha labeling (Figs. 2c, 2d). Equivalent loading for each western blot was confirmed with Ponceau S staining. While not as precise as MS quantification, western blot studies allowed direct comparison of labeling across the measured time interval without performing Aha-labeled protein enrichment experiments that are required prior to MS analysis. 20 This allows for relative turnover rates of newly synthesized proteins to be fit (Fig. 2). The degree of fluorescent signal normalized relative to the background (rF) measured in relative fluorescence units (RFU) by semi-quantitative western blotting was assumed to be linearly proportional to the concentration of pAha using a fluorescent labeling factor, k f . For each tissue (denoted 'x')

Aha Labeling Captures Relative Protein Synthesis and Turnover Dynamics in Murine Tissues
The highest plasma concentration of fAha was recorded at 0.5 hpi (Fig. 2a), indicating a rapid absorption and bioavailability from the injection site. Also notable, volume of distribution (V D ) was calculated to be 1.6 mL g −1 body weight, a value greater than total body volume, indicating rapid sequestration of fAha into body tissues. Tissue profiles of fAha also peaked early, between 0.5 and 1 hpi, with the liver having the earliest peak compared to the other tissues (Figs. 3a-3d). This early peak can be attributed to the high blood perfusion of the liver, 30 which likely results in faster distribution equilibrium of Aha into the liver compared to other tissues. Clearance of fAha was also observed to be rapid, with a plasma half-life (t 1/2 ) of only 1.9 h and a plasma clearance (CL) of 12.1 mL h −1 .
In all tissues, maximum protein labeling was observed around 6 hpi (Figs. 3e-3h), delayed from that observed for [fAha]. Importantly, across tissues, [fAha] is directly proportional to the first derivative of [pAha] in the early time domain, implying that protein synthesis dominates the mass balance at low [pAha], a finding supportive of our model's first-order kinetics. However, the degree of labeling, represented by the maximum fold increase in fluorescence intensity compared to an internal control, and the kinetics of protein incorporation and turnover varied considerably between tissues.
The observed differences in fluorescence intensities are in agreement with previous stable isotope labeling studies that showed faster protein turnover rates in liver and kidney compared to brain and skeletal muscle. [31][32][33] Our model estimated a protein lifetime in the brain that is 2.7 times higher than the liver (Table 1), a ratio in close agreement with previous findings (Price et al. observed a median peptide half-life of 9 and 3 days for brain and liver, respectively 32 ). These discrepancies between the magnitude of previously reported values and the half-lives estimated here may be attributed to the shorter timescale of our experiment as compared to stable isotope labeling techniques and the use of semi-quantitative western blotting measurements rather than MS. A more robust quantitation of pAha concentrations using MS will be required for a more precise determination of absolute protein kinetics in different murine tissues.

Descriptive Statistics Support Model of Biodistribution and Labeling
Model validity was examined using the following metrics to investigate parameter stability and goodness of fit. First, a prediction interval of residuals was generated for each tissue studied, for both the biodistribution and protein labeling models. A 95% prediction interval (PI α=0.05 ) of residuals was calculated using all experimental replicates (Fig. 3). A second, tighter prediction interval was calculated using the mean for each time point (Fig. 3). For each tissue, the width of PI α=0.05 from the average line of best fit ( ŷ ) can be approximated using a naively informed forecast interval that assumes a normal distribution of residual error. 34 where n is the total number of observations, (t i , y i ) are the coordinates for each observation, t is the time point of the predicted residual, t is the average time of all experimental observations. These prediction intervals demonstrate a narrow range of residual error and capture all experimental means and > 95% of experimental values, indicating an accurate predictive model. Second, the covariance matrix of least squares regression was used to inform a standard error of fitting (SE f ) for all fitted parameters in the biodistribution model (Additional File 1: S3-S5). SE f describes variability of each parameter, but also reflects upon the definition of the model. Parameters with best fit values near to the constraints are less predictable; wide error would be indicative of a poorly constrained system. Error also increases with the number of fitted parameters, is inversely related to the quantity of data points available for fitting, and is weighted by the metric used for optimization. Among all fit parameters, we found relatively few SE f values greater than best fit values by more than an order of magnitude (Fig. 3, Additional File 1: S3-S5). The widest error ranges were associated with parameters rates describing tissue import and export of fAha, which are the least well characterized in the literature and therefore the most naively bounded in our model. Low SE f values, particularly  among parameters fit within reported literature values, support a well characterized model, although this metric is not fully sufficient to describe a complex non-linear regression. A more thorough analysis was achieved using Monte-Carlo Latin hypercube sampling (LHS) to perform efficient sampling of the input parameter space and correlation with partial rank correlation coefficients (PRCCs). This global sensitivity analysis was performed against several metrics to probe the influence of variation in parameter values on (1) model fitness and (2) predicted Aha labeling levels in each tissue (rF). To address model fitness, PRCC values, which vary between 1 (perfect positive correlation) and − 1 (perfect negative correlation), were calculated comparing variation in the input parameter values against the sum of square errors for each fitting model. This analysis generates a PRCC value for each parameter that characterizes its relative effect on the fitting metric (PRCC f ), and therefore the fitting process. PRCC f values higher in absolute magnitude indicate parameters that, when varied independently, have the greatest influence on the fitting metric with values > 0.8 considered as indications of regions of instability in the model. 35,36 All model parameters fell safely below a PRCC f value of 0.3. Parameters found to have the highest PRCC f values (> 0.2) were those that influence the system removal of fAha (e.g. elimination and import rates into metabolic tissues), likely due to the rapid metabolic profile exhibited following subcutaneous injection (Fig. 3).
Taken together, PRCC f and SE f values demonstrate that influential parameters are not those with the widest error ranges (Fig. 4, Supplements S5-S7). Furthermore, when comparing these values across the parameter space, there is no one tissue that exhibits over-sensitive behavior by either metric (apart from perhaps systemic elimination). These findings provide support for the model definition, boundary constraints, and biological relevance of best fit parameters.
Further application of PRCC analysis allowed examination of the influence of parameter values upon the predicted degree of Aha labeling in each tissue over time. All model parameters were varied within their estimated fitting range (Additional File 1: S5-S7) and sampled via LHS (n = 10,000). An rF value in each tissue was predicted for each timepoint from 0 to 24 hpi and the resulting output was correlated to the input parameter variation using PRCC analysis. Resulting traces elucidate kinetic parameters and physiological mechanisms that most influence rF during different time domains (Fig. 5). For example, immediately after injection (0-4 hpi), the rate of absorption of Aha from initial injection site (kaSubcu) is a critical mechanism driving rapid labeling efficiency in all four studied tissues. However, the importance of this parameter declines over time, just as the influence of other parameters increases, particularly those representing Aha degradation and elimination. The influence of liver transport is prominent in all four tissues indicating the liver plays a major role in Aha kinetics. This suggests Aha may be primarily metabolized and eliminated by the liver, rather than by eliminated via excretion. This possibility is consistent with the known liver functions of toxin removal and amino acid metabolism. 37

Predictive Simulations Accurately Capture Alternate Dosing Paradigms
Attaining sufficient protein labeling is critical for accurate identification and quantitation of Aha-labeled proteins using LC-MS/MS. For instance, note the near negligible labeling of skeletal muscle in Fig. 3. If muscle labeling is desired, a much higher Aha dose may be required to attain sufficient labeling. This underlines the importance of tailoring the dosing regimen of Aha (i.e. amount per dose and dosing frequency) to the tissue of interest and the biological question under investigation. Using the model described above, fAha biodistribution and tissue protein labeling can be predicted for alternative dosing regimens to aid future experimental design and predict labeling efficiency depending on the conditions of a study. However, the model was based on the concentration of Aha in plasma and tissues over 24 h after a single subcutaneous injection. While this data was sufficiently robust to develop a well-informed model of Aha patterning after one subcutaneous dose, the use of this model to predict Aha content for longer time scales or multiple dose paradigms would generate naïve forecasts. To address this gap, we validated the predictions of our model experimentally against a data set that spans a longer time period and multiple injected doses. The model was used to predict the rF labeling of brain and liver tissues for two alternative dosing paradigms with either 12 h repeated doses (hrd) or 24 hrd over a 32 h period (Fig. S7). As an internal control for western blotting variation, these new experimental values were normalized by a shared time point with the previous study (6 h post initial injection) for each tissue (Fig. 6).
To determine the ability of the original model to predict Aha incorporation into proteins in various tissues, the data from each repeated dose study was used to refit the relative pAha synthesis and degradation rates for each tissue under each repeated dose paradigm. Relative to the parameter fit with the original experimental data, there was only a slight reduction in the standard error of regression (SE reg ), a goodness-of-fit metric, between the original and refit parameters in each tissue (Table 2). Additionally, among all refit parameters, a single parameter was adjusted beyond a single standard error of fit (SE f ) from the original best fit value (12 hrd liver Δ(k f k s ) ≈ +1.97SE), and only the degradation rate in the brain changed by > 20% (Table 2).

Aha Administration Induces Metabolic Adaptation in Mice
In addition to the characterization of Aha distribution kinetics in mice, it is important to understand the impact of Aha administration on normal physiology. Since metabolites are the end products of cellular biological processes, we sought to evaluate the physiological impact of Aha incorporation into NSPs using untargeted metabolomic analysis. To this end, plasma and liver metabolomes were analyzed following different time points (6,18, and 32 hpi) of two alternative Aha dosing regimens (12 hrd and 24 hrd) to examine the effect of Aha administration on the metabolome. LC-MS analysis of plasma and liver metabolomes (Additional Files 2 and 3, respectively) identified a total of 1339 and 1041 mass features (i.e. metabolites), respectively across all sample groups. The peak area of each mass feature is proportional to the amount of the corresponding ion in the sample and was used as a measurement for its relative abundance across samples. Principal component analysis (PCA) was used to examine the variation in metabolite abundance between the different treatment groups. Plasma PCA analysis revealed no distinct segregation between the non-injected control mice and the Aha 6 hpi mice, indicating that there were no global differences in the metabolome between the two groups (Fig. 7a). Greater separation was observed between the control and the other Aha treatment groups; however there was less variation among each of these Aha groups and one another.
In addition to PCA, unsupervised hierarchical clustering analysis (HCA) was conducted and a heatmap was generated to examine variations in metabolic patterns between the Aha and control groups. HCA and heatmap visualization of the identified plasma mass features showed that the biological replicates of each group did not cluster with each other and that there were no distinct differential abundance patterns among the different groups (Fig. 7b). These results suggest that minor metabolic alterations occur after 6 h of Aha injection and become more prominent following longer post injection periods and multiple Aha doses. Additionally, while Aha does induce metabolic adaptation in injected mice, no substantial metabolic changes exist between the various groups.
Similar trends were observed for the liver metabolome. PCA analysis indicated less variation between the control  and the Aha 6 hpi groups as well as between all other Aha treatment groups, whereas greater separation was shown between these treatment groups and the control (Fig. 7c). Additionally, a heatmap with HCA of liver mass features illustrated a lack of group clustering based on metabolite abundance (Fig. 7d). As such, the LC-MS analysis of the liver metabolome further establishes that the control and the Aha 6 hpi mice have similar metabolite profile and that there are no significant differences between 12 and 24 hrd over a 32 h period. Following global analysis using PCA and HCA, Student's t-tests were employed to compare between the control and each of the Aha treatment groups. Using a raw p-value of 0.05 and a fold change of ≥|2| as cut-offs, the percentage of differentially abundant mass features ranged between 3.66% and 10.90% for plasma, and 1.06% and 18.73% for liver (Table 3). These values dropped to 0% and 4.48% for plasma and 0% and 11.82% for liver using a FDR-adjusted p-value, which has been suggested for untargeted metabolomics data to control for false positives 38 . Notably, these minor metabolic alterations are in agreement with a recent study that investigated the metabolic effect of growing E.coli in media supplemented with ncAAs 39 .
Collectively, the LC-MS analyses of plasma and liver metabolomes demonstrate that Aha administration causes metabolic adaptation in mice. While 6 h following a single Aha injection is a relatively short time to detect a difference, more distinct metabolic changes occur after prolonged post injection periods and repeated Aha injections. PCA, HCA and t-test results refer to a limited metabolic perturbation occurring in response to Aha administration. Notably, our labeling technique does not involve Met restriction or depletion as the case with other labeling strategies that use an Aha-supplemented Met-free diet. Met dietary restriction has been shown to alter the metabolism in mouse models and in humans. [40][41][42] Yet, subsequent MS/MS identification studies of the differentially abundant metabolites are needed to understand the implication of these small changes and the potentially impacted metabolic pathways.

Discussion and Conclusions
To optimize the applications of ncAAs, effective protein labeling is critical to enrich labeled protein with the high signal-to-noise ratio required for accurate quantitative MS measurements and identification of newly synthesized proteins. Here, we reported for the first time the biodistribution of the widely used Met analog, Aha, in murine tissues, as well as the associated relative rates of incorporation of Aha into protein via protein synthesis and protein turnover.
These results showed that liver and kidney have faster protein synthesis and turnover rates compared to brain and skeletal muscle, which is consistent with previous studies that utilized stable isotope labeling. 32 In previous work, we found that protein labeling with Aha was dose dependent and the optimal dose via intraperitoneal injection (IP) of Aha was 0.1 mg g −1 mouse. 14 In this current study, we utilized the same dose with SC injection at various timepoints. This indicates that protein labeling with Aha is robust to this dosing level, dosing schedule, and to injection route. Further, we demonstrated that under a subcutaneous injection dosing paradigm protein labeling reaches a local maximum within a relatively short window in all studied tissues (~ 6 h). We expect this time resolution to enable study of proteins with shorter half-lives. In contrast, alternative dosing of ncAA (e.g. oral feeding over a period of days to achieve desired Aha labeling 13 ) or isotope-labeled amino acids 32,33 are typically administered over longer time period. Other ncAAs such as homopropargylglycine (Hpg) and azidonorleucine can be used for in vivo labeling. 2 The advantage of using Aha and Hpg is that they can be incorporated into proteins in native systems using endogenous translation machinery. Our group and others have previously shown that Aha incorporates into proteins more efficiently than Hpg. 14,29 In order to label proteins with azidonorleucine and many other ncAAs, a mutant tRNA synthetase recognizing the ncAA is required. 14,15 While Aha (and other ncAA) labeling shows great promise for future study of the nascent proteome, the biodistribution of Aha makes it clear that depending on the tissue type and biological processes to be studied, multiple injections of Aha may be required. Therefore, optimizing the dose and frequency of Aha injections is critical for the appropriate design of labeling studies. In this regard, we also presented an experimentally informed deterministic model. Our model was developed to describe: (1) the transport of fAha within the plasma, (2) the circulatory exchange of fAha into tissues of adult mice and (3) the degree of Aha incorporation into proteins in specific tissues. We utilized a multi-compartment system characteristic of pharmacokinetic models of small-molecule and AA pharmacokinetics. 23,27,28 Such models have been extensively used to study protein turnover, typically in the context of stable isotope labeling, and have been used to estimate kinetic parameters such as total protein degradation in isolated tissues (e.g. Guan et al., employed isotope labeling and a compartment model of similar structure to estimate kdLiver ≈ 0.03 d −1 , 28 in close agreement with our kdLiver ≈ 0.02 d −1 ). Our model was empirically parameterized with time resolved measurements of Aha concentrations and fluorescent labeling and used to extract the distribution kinetics of Aha in murine tissues and their relation to the degree of protein labeling as well as to compute the relative rates of protein synthesis and turnover. We Table 3 The number of differentially abundant plasma and liver mass features for each Aha treatment group compared to control using either a raw p-value (0.05) or FDR-adjusted p-value (0.05) and a fold change of ≥|2|.

Group
Plasma Liver Raw p-value FDR Raw p-value FDR   6 hpi  49  2  11  0  12 hrd_18 hpi  88  0  97  22  12 hrd_32 hpi  157  60  159  68  24 hrd_18 hpi  145  43  154  60  24 hrd_32 hpi  146  17  195  123 further validated this framework for predictive modeling of Aha labeling against an experimental dataset including two different repeated injection dosing paradigms to demonstrate its efficacy as a tool for future experimental design. Of note, semi-quantitative fluorescent gel imaging was used to inform the relative change in the amount of NSPs in each tissue compartment, as has been previously demonstrated. 43 For absolute measurement of individual protein synthesis rates, future studies will need to employ quantitative proteomic analysis.
With respect to the domain of validity for our model, we note that our study only probed subcutaneous dosing at a single dose concentration optimized in previous work (0.1 mg g −114 ) and that Aha concentrations are recorded on a limited timescale (up to 32 hpi). To extrapolate beyond these limitations, we can consider Aha from a pharmacokinetic perspective. At low [fAha], adequate homeostatic regulation should compensate for an increased rate of administration with a proportional increase in the rate of elimination. However, as concentrations and dosing rates continue to rise a further domain might exist where the [fAha] exceeds regulation capacity and begins to demonstrate saturation of elimination mechanisms. 44 If fAha dosing enters this saturation domain, we expect a deviation from the first-order dynamics towards a zero-order system. However, at our highest dosing frequency (12 hrd), we observed no evidence of a deviation from the first-order kinetics. The model predicts both fAha and protein labeling without variation in the rates of absorption or elimination. Therefore, we expect this model to be valid at dosing rates up to the maximal observed clearance of fAha, ~ 30 μg h −1 per g body mass. At our dose concentration, this would correspond to injections administered at ~ 3 hrd intervals, a rate that approaches the limit of feasibility when weighed against the stress burden of multiple injections in a rodent model. Furthermore, to adapt the model to capture an oral dosing study, we suggest refitting the distribution model with plasma [fAha] from at least 2 time points using feed/water supplemented with Aha. Refitting can reasonably be confined to relevant parameters (such as gastrointestinal absorption (kiGasin), excretion (keGasin), blood flow (qbSmint) and all systemic elimination parameters) and should employ a realistic forcing function for oral dosing. Lastly, for studies on a longer timescale than 48 h, we suggest refitting the protein labeling model for desired compartments with at least 2 time points on the desired time scale using the best fit values for each parameter as initial estimates.
Our model demonstrated that Aha exhibits high rates of plasma clearance relative to a high volume of distribution (Fig. 2a). This supports the hypothesis that elimination of Aha predominantly occurs within metabolically active tissues, such as the liver. This is also supported by the prominence of liver parameters in the model analysis (Fig. 5). Further, the liver showed the highest degree of labeling (Fig. 3f) as well as the highest relative rates of Aha incorporation and protein turnover (Table 1), whereas skeletal muscle had the lowest degree of labeling (Fig. 3e) and the slowest relative rates of incorporation and turnover (Table 1). Interestingly, the amount of fAha (μg) per unit mass tissue (g) was higher in skeletal muscle than in liver (Figs. 3a, 3b), indicating that the low degree of labeling observed in skeletal muscle is predominantly due to a slow rate of muscle Aha incorporation, perhaps implying a lower rate of overall protein synthesis on this timescale. Finally, we investigated the impact of Aha administration on plasma and liver metabolomes. Despite the growing use of ncAAs, and Aha in particular, evaluation of the physiological impact of Aha administration, has previously been limited to examining changes in animal gross behavior, physical appearance and body weight. [13][14][15] A more robust analysis of the effect of Aha administration on mice fed a MAD and the corresponding implications for cellular function, is still required to confirm the suitability of the method for longer-term in vivo studies. Here we have demonstrated that Aha incorporation into cellular proteins induces minimal metabolic perturbations in injected mice. This observation further confirms previous results from our group that demonstrated that ncAAs do not affect the gross behavior nor the physical appearance of treated mice. 14 Taken together, we expect the biodistribution model to be a useful tool to guide future experiments for maximizing Aha labeling of newly synthesized proteins. This also provides a platform for investigation of other biodynamics of other ncAAs and other delivery methods (oral or biodegradable depots).

Animal Model
Animals used in these studies were derived from female age-matched (9-11 weeks) wild-type adult C57BL/6 mice (Mus musculus) purchased from The Jackson Laboratory. All experimental protocols were performed in compliance with established guidelines and all methods were approved by Purdue Animal Care and Use Committee (PACUC, pro-tocols# 1209000723 and 1801001682). PACUC requires that all animal programs, procedures, and facilities at Purdue University abide by the policies, recommendations, guidelines, and regulations of the United States Department of Agriculture (USDA) and the United States Public Health Service (USPHS) in accordance with the Animal Welfare Act and Purdue's Animal Welfare Assurance.

Aha Injection, and Plasma and Tissue Collection
l-azidohomoalanine (Aha; Click Chemistry Tools) was resuspended in 1 × phosphate buffered saline (PBS) to 10 mg mL −1 , adjusted to pH 7.4, sterile filtered and stored at − 20 °C. All Aha injections were administered subcutaneously at 0.1 mg g −1 total mouse weight. Mice (n = 3, biological replicates) were euthanized 0.5, 1, 2, 4, 6, 12 and 24 h post injection (hpi). Blood was harvested by cardiac puncture, collected in EDTA-treated tubes and centrifuged at 1500×g for 10 min at 4 °C. The supernatant (plasma) was transferred into a new tube, snap frozen in liquid nitrogen and stored at − 80 °C. Liver, brain, kidney and hindlimb skeletal muscle tissues were dissected at each time point, snap frozen in liquid nitrogen and stored at − 80 °C. Control plasma and tissues were collected as described above from non-injected mice (n = 3 biological replicates). For the validation of model predictive ability, two Aha dosing regimens were used: (1) 12 h repeated doses (hrd) and (2) 24 hrd. Liver and brain tissues (n = 3 biological replicates) were dissected as described above at 6, 18, and 32 hpi, snap frozen in liquid nitrogen and stored at − 80 °C.

Sample Preparation for Aha Analysis
For plasma sample preparation, 50 µL of plasma were mixed with 10 µL of 1 × PBS, pH 7.4, and 5 µL of 100 ng µL −1 l-α-aminobutyric acid (α-ABA; Sigma Aldrich) that was used as an internal standard. 12.5 µL of trichloroacetic acid (TCA; Sigma Aldrich) were added to the mixture to precipitate proteins. The mixture was incubated for 10 min at 4 °C and centrifuged at 16,000×g for 10 min at RT. The supernatant was then mixed with 100% acetonitrile (ACN; Fisher Scientific) at a 1:1 ratio (v/v). The mixture was transferred to an HPLC autosampler vial for LC-MS/MS analysis. For calibration curve generation, Aha standards were prepared by mixing 50 µL of non-injected plasma with 10 µL of a known concentration of Aha and 5 µL of α-ABA. Proteins were then precipitated with TCA and prepared for LC-MS/ MS analysis as described above.
For tissue sample preparation, tissues were rinsed with ice-cold 1 × PBS, pH 7.4 to remove residual blood and homogenized in ice-cold 1 × PBS, pH 7.4 using a TissueRuptor (Qiagen). The final homogenate weight was measured and converted to volume by using a homogenate density of 1 g mL −1 . Samples were then prepared for LC-MS/MS analysis as described for plasma by using 50 µL of the tissue homogenate. The remaining plasma samples and tissue homogenates were snap frozen and stored at − 80 °C until use for western blot and untargeted metabolomic analyses as described below.

LC-MS/MS Targeted Analysis of Aha
An Agilent 1260 Rapid Resolution liquid chromatography (LC) system coupled to an Agilent 6470 series QQQ mass spectrometer was used for Aha analysis (Agilent Technologies). An Intrada Amino Acid 2.0 mm × 150 mm, 3.0 µm column (Imtakt Corporation) was used for LC separation.  Table 4. The jet stream ESI interface had a gas temperature of 325 °C, gas flow rate of 9 L min −1 , nebulizer pressure of 35 psi, sheath gas temperature of 250 °C, sheath gas flow rate of 7 L min −1 , capillary voltage of 3500 V in a positive mode, and nozzle voltage of 1000 V. The delta electron multiplier voltage was 300 V. Agilent MassHunter Quantitative Analysis software was used for data analysis (v.8.0).

Kinetic Modeling of Aha Distribution
The model of Aha distribution is shown schematically in Fig. 1. Both injection and oral administration of Aha are included in the model system to enable future studies of ncAA labeling with alternative/combined dosing paradigms. Simulations were performed using custom modeling scripts written in Python 3.6 (code location in Additional File 1: S1). Systems of ordinary differential equations (Additional File 1: S2) were solved using a flexible high order solver from the SciPy python package (v1.4.1 45 ). Simulations were run on a Lenovo Yoga PC with an Intel Core i7-8550U CPU @ 1.8 GHz and 8 GB RAM.
Most parameter values and ranges for fitting were informed from reported literature values or experimental measurements from this study. For parameters related to an experimental output ([fAha] or rF) without a reported literature value, an initial best estimate was selected to produce a single time-step change one order of magnitude lower than the maximum recorded experimental value. These parameters were then allowed to fit within a range of 1.5 orders of magnitude from the initial estimate. Parameters were fit with a least squares minimization algorithm from 'Lmfit', a prebuilt python library (v1.1.0 46 ). All best fit values and boundary conditions can be found in the parameter tables (Additional File 1: S3-5). Other package dependencies to generate data and figures in this manuscript include NumPy (v1. 16.1, 47 ), and matplotlib (v3.3, 48 ).

Parameter Sensitivity Analysis and Model Validation
To effectively sample the input parameter space, Latin hypercube sampling (LHS) was utilized to generate unique parameter sets (n = 10,000), sweeping each parameter value through a range defined by the boundary constraints from literature (Additional File 1:Tables S3-S5) as previously detailed. 35,36 Global sensitivity analysis was performed on the LHS generated parameter sets using partial rank correlation coefficient (PRCC) analysis. This analysis quantifies the sensitivity of an output variable on the variation in input parameter values. 35,36 Here, PRCCs were determined for each of the 19 fitted parameters (Additional File 1: Table S3) and 12 static parameters (Additional File 1: Table S4) in the biodistribution model, as well as for all 8 parameters in the protein incorporation model (Additional File 1: Table S5). PRCCs were used to characterize the influence of each parameter on the sum of square errors (SSE), the optimization metric for non-linear regression. Simulations used to inform PRCCs were performed on the Brown Supercomputing Community Cluster at Purdue University, 49 with each simulation run on a single node with dual 12-core Intel Xeon Gold "Sky Lake" CPUs @ 2.60 GHz and 96 GB of memory.
The standard error of fitting was determined for the 19 fitted parameters (Additional File 1: Table S3) in the biodistribution model and for all 8 fitted parameters (Additional File 1: Table S5) in the protein incorporation model. Standard error values were determined from the covariance matrix during non-linear regression using the built-in functionalities of the 'Lmfit' python library. 46

Plasma and Liver Sample Preparation for Untargeted Metabolomic Analysis
Plasma and liver metabolomes of non-injected control samples (n = 3 biological replicates) and samples collected at different time points post Aha injection (n = 3 biological replicates per time point) were extracted by adding methanol: chloroform: water (1:1:1 v/v) to 80 μL of each plasma sample. Samples were vortexed briefly and centrifuged at 8000×g for 5 min at RT. The upper layer was transferred into a new tube and vacuum-dried overnight. The dried fraction was reconstituted in 75 μL 5% ACN and 0.1% FA. Reconstituted samples were sonicated for 5 min, centrifuged at 16,000×g for 8 min at RT, and the supernatants were transferred to HPLC autosampler vials.

Untargeted LC-MS Metabolomic Analysis
Separations were performed on an Agilent 1290 UPLC system (Agilent Technologies). Metabolites were analyzed using a Waters Acquity HSS T3 column (1.8 µm, 2.1 × 100 mm), with a mobile phase flow rate of 0.45 mL min −1 , where the mobile phase A and B were 0.1% FA in double distilled water and ACN at a 1:1 ratio, respectively. Initial conditions were 100:0 A:B, held for 1 min, followed by a linear gradient to 20:80 at 16 min, then 5:95 at 22.5 min. Column re-equilibration was performed by returning to 100:0 A:B at 23.5 min and holding until 28.5 min.
Mass analysis was obtained using an Agilent 6545 Quadrupole Time of Flight (Q-TOF) MS with ESI capillary voltage + 3.2 kV, nitrogen gas temperature 325 °C, drying gas flow rate 8.0 L min −1 , nebulizer gas pressure 30 psig, fragmentor voltage 130 V, skimmer 45 V, and OCT RF 750 V. MS data scans (m/z 70-1000) were collected using Agilent MassHunter Acquisition software (v.B.06). Mass accuracy was improved by infusing Agilent Reference Mass Correction Solution (G1969-85001). MS/MS was performed in a data-dependent acquisition mode on composite samples. Raw data for plasma and liver metabolomics is found in supplemental Additional Files 2 and 3, respectively.

Metabolomic Data Statistical Analysis
Peak deconvolution and integration were performed using Agilent ProFinder (v.10.0). Bioinformatic analyses were performed using Agilent Mass Profiler Professional (v.13.1). Chromatographic peaks were aligned across all samples. Missing values (3 in plasma and 2 in liver) were replaced by the lower limit of detection (1/5 of the minimum positive value of each variable). Data were normalized by logtransformation and auto-scaling. Mass features with P < 0.05 and fold change ≥ 2 were considered significant. Principal component analysis (PCA) and hierarchal clustering analysis (HCA) were performed using MetaboAnalyst v.5.0.