A Novel Pulse-Chase SILAC Strategy Measures Changes in Protein Decay and Synthesis Rates Induced by Perturbation of Proteostasis with an Hsp90 Inhibitor

Standard proteomics methods allow the relative quantitation of levels of thousands of proteins in two or more samples. While such methods are invaluable for defining the variations in protein concentrations which follow the perturbation of a biological system, they do not offer information on the mechanisms underlying such changes. Expanding on previous work [1], we developed a pulse-chase (pc) variant of SILAC (stable isotope labeling by amino acids in cell culture). pcSILAC can quantitate in one experiment and for two conditions the relative levels of proteins newly synthesized in a given time as well as the relative levels of remaining preexisting proteins. We validated the method studying the drug-mediated inhibition of the Hsp90 molecular chaperone, which is known to lead to increased synthesis of stress response proteins as well as the increased decay of Hsp90 “clients”. We showed that pcSILAC can give information on changes in global cellular proteostasis induced by treatment with the inhibitor, which are normally not captured by standard relative quantitation techniques. Furthermore, we have developed a mathematical model and computational framework that uses pcSILAC data to determine degradation constants kd and synthesis rates Vs for proteins in both control and drug-treated cells. The results show that Hsp90 inhibition induced a generalized slowdown of protein synthesis and an increase in protein decay. Treatment with the inhibitor also resulted in widespread protein-specific changes in relative synthesis rates, together with variations in protein decay rates. The latter were more restricted to individual proteins or protein families than the variations in synthesis. Our results establish pcSILAC as a viable workflow for the mechanistic dissection of changes in the proteome which follow perturbations. Data are available via ProteomeXchange with identifier PXD000538.


Isotope patterns of doubly labeled peptides
pcSILAC labeling produces isotope patterns more complex than standard SILAC in the case of peptides containing more than one labeled amino acid, i.e. 2 x K-, R+Kor 2 x R-containing peptides resulting from missed trypsin cleavages (including RP, KP sequences).
"Pure" double-labeled (2 x K, 2 x R) peptides were used for pcSILAC quantification after we manually inspected their isotope patterns at t=6h,12h in both samples before and after mixing. Indeed, such peptides could be showing partial labeling peaks, i.e. peptides carrying one heavy and one light Lys, deriving from label mixing ("old" label with "new" label) occurring especially in the first hours after medium exchange (Supplementary Figures S2-3 and S5-6). Analysis of doublelabeled (2 x K and 2 x R) peptides has been previously used to quantify the amount of label mixing [1] [2]. We found that, in general, partial labeling peaks were of negligible intensity (<5% of base peak, Supplementary Figures S5-S6), indicating that the replacement of the intracellular label, when it occurred, was nearly complete and thus label mixing did not influence the quantitation obtained at the time points considered.  Figure S3). We verified MS signals for 20 distinct peptides from various proteins at 6h and 20h and observed that their isotope patterns were indeed complex and did fit the predictions (Supplementary Figure S5). Like for pure labeled peptides, signals derived from label mixing (old+new label) were of low or very low intensity. Interestingly, such mixed KR-or RK-containing peptides could theoretically be used for quantitation purposes since they contain specific information e.g. on total protein amounts in the two samples (each peak corresponds to the four labeling species R6K0, R10K0, R0K4, R0K8). However, their complex isotope patterns are overlapping and difficult to disentangle. To our knowledge, at the moment there is no software tools that can exploit the information contained in these signals for quantitation.

Incorporation ratios of K-vs. R-containing peptides
Since R and K are used independently, we tested if there were differences in the ratios measured on peptides containing one or the other amino acid. Such differences could arise, for example, due to different incorporation kinetics. For this purpose we analyzed pcSILAC cell extracts at the three time points from control and GA-treated cells before mixing, to measure new/old ratios which express label incorporation ratios. Non-normalized (H/L) or (M/L) ratios obtained with R and K showed highly correlated incorporation kinetics (Supplementary Figure S4) in both conditions with linear correlations coefficients R2 above 0.95 and intercepts virtually at 0. However, the slope in both treated and control samples was reproducibly between 0.8 and 0.9 at all time points, meaning that for the same protein the ratios measured on K-containing peptides were systematically lower than those measured on R peptides. This cannot be due to a different size in the initial pool of amino acid which needs to be exchanged, since this would have given rise to a delay in time (a non-zero intercept on the plot) rather than a different slope.
A possible interpretation would be that of a persistent (until t=20h) defect in incorporation of new label for K compared to R, for example through the recycling of K from degradation of preexisting proteins, which would constantly release "old" label. Such a phenomenon has already been postulated by Boisvert et al. [3], who also reported values close to 20%, without analyzing separately incorporation of R and K. However, amino acid recycling was not compatible with our observation that the degree of label mixing assessed on the spectra of multiply-labeled peptides was minimal in our samples. Furthermore, examination of some standard SILAC datasets acquired on the same cells, after an identical sample processing with separation of the values for K,R, revealed that ratios calculated before normalization showed a similar bias, with K-based ratios 10-12% lower than Rbased ratios. Such differences typically disappear in normal SILAC data analysis because the MaxQuant software normalizes peptide ratios separately for R-and Kcontaining peptides. The observation of a similar bias also in standard SILAC samples which are 100% labeled tends to exclude amino acid pools and incorporation mechanisms as a possible cause. The reason(s) for these differences between R and K ratios are presently unknown to us. We can speculate that phenomena related to the isotopic labeled amino acids used may play a role. The satellite peak at -1 a.m.u visible for K4-and K8-labeled peptides (Supplementary Figure S6) could somewhat cause an underestimation of the medium or heavy K isotopes. Alternatively, an isotope effect in trypsin digestion similar to what recently described [4] could be responsible for this difference. Different rates of cleavage for the light, medium and heavy labels of R-and K-peptides were therefore integrated in the kinetic model based on pcSILAC data (see below). Systematic biases linked to the amino acid are occasionally detected in standard SILAC but are systematically removed through normalization, which for this reason is performed separately by MaxQuant ( [5] and J.Cox, personal communication).  S7C). The lower r value for R could be explained by the much smaller spread of the data for (H/M) R , i.e. by the fact that most values of (H/M) R are very close to 1 and therefore the correlation is strongly influenced by technical noise, which tends to be constant. Indeed, examination of the values of ratio (H/L) R showed that for other Rbased ratios correlation between experiments and replicates was good (r >0.85, Supplementary Fig. S7). These observations also suggest that using all six ratios as an input for further calculations could be an advantage as it may permit to correct for measurement errors and increase robustness of fitting.

Mathematical formulation of variables linked to pcSILAC experiments
Variables and indexes were defined as follows (see main text for description): For control (untreated) cells: For simplicity, the indexes and time dependence will be omitted in most equations.
The following relationships can be derived: Similarly, by inserting (2) into (3) or (1) into (4) the same can be obtained for the treated (heavy) sample : Now, continuing to work with these equations it is possible to obtain formulas to calculate the total ratio in the form of two very symmetrical expressions:

Calculation of pseudo-absolute protein concentrations in the control sample (p A ) using iBAQ scores
The use of the integrated peptide ion intensities generated by the MaxQuant software to estimate protein concentrations has been described previously (iBAQ) [6]. The standard output table proteinGroups.txt generated by MaxQuant contains the signal intensities for the three channels L, M and H for every protein, regardless of the sequence of the peptide (K,R peptides are counted together). Given the characteristics of pcSILAC data, the M and H channels are, respectively, related to the total amounts of protein in the control (U+V) and treated (X+Y) samples : Where I M and I H are the intensities of the signal in the corresponding channels. For example, the iBAQ score for a protein i in the sample A (control) is thus Where n det is the number of detectable peptides for the protein i [6]. Since no data for absolute concentrations of internal standards were available in our datasets, we expressed the concentration of proteins in arbitrary units relative to the median of all iBAQ scores at one time point. We verified that the values obtained in this way for the sample A (control) were very similar at the 3 time points (data not shown) after which we obtained a general value for the concentrations p i,A as: It is important to remember that the concentrations in sample A (untreated cells) are assumed to be constant since the system is considered to be at steady-state.

Model and calculation of kinetic parameters
The meaning of the individual variables linked to protein concentrations and their evolution can be summarized as follows : The time of exchange of label was defined as t=0, while t d is the time of start of drug treatment. First order decay and synthesis kinetics were postulated. The four variables u,v,x,y describe the concentrations of proteins synthesized before (u,x), resp. after label exchange (v,y). V i (with i=A,B) is the rate of synthesis, while k i,app (with i=A,B) is the apparent decay rate constant, i.e. the constant determined from the data before correction by the cell growth rate.
Here we present in more detail the derivation of some equations, the explicit solutions of the ODE's and further parameters of the model.
Here we will show how the equation (1) from the main text is derived from the mass balance equations. The total mass balance describes the evolution of the abundance of a given component, here for a protein: where P A is the abundance of the protein in one cell (with index A to refer to the control cell cultures; we emphasize that P A is an abundance and not a concentration, so the units of P A could be, for example, moles), V A c is the volume of the cell, V A is the synthesis rate for this protein species (in mol L  s       ) and k A,d is the protein decay rate constant.
But then, instead of looking at the abundance, we consider the concentration of the protein (p A ), which is directly related to the abundance by Taking the time-derivative of (17) and applying the chain rule gives: where  A is the growth rate (or dilution term).
Combining equations (16) and (18), keeping only the concentration-derivative on the left and dividing by the volume of the cells gives the wished equation: The calculation of the growth rate  A is described below in the "Additional notes".
For all subsequent steps we will thus use the apparent decay rate, k A,app .

Equations of evolution of protein concentrations
Then we can solve the ODE equation (19) and similar ones that represent the evolutions of the proteins with different markers. Considering the evolution of the various variables u,v,x,y as described above and assuming steady state at t < 0 we And for example, for t > 0 Solving the equations for u and v gives u t Now, x and y follow the same evolution than u, v, until t d , so at t= t d : And solving the evolution after t d gives: x t So, equations (23) and (25) describe the evolution of the proteins in both cultures until cell lysis and protein extraction.

Normalization, calculation of mixing ratios and data correction
In standard proteome analyses by SILAC, the peptide ratios obtained are typically normalized on the basis of the median of the global distribution of all ratios, assuming that the overall composition of the two samples is comparable and that when the median of ratios is not 1.0 the difference is caused by imperfect (not 1:1) mixing of the two extracts. This is however an approximation, because the median of the individual protein ratios does not necessarily coincide with the actual total protein mixing ratio, which is based on an experimental measurement of the total protein concentrations of a sample. For example, the total protein mixing ratio can be influenced by the changes in a few very abundant proteins, while the same changes will not significantly affect the median of protein ratios. Nevertheless, in most standard SILAC experiments global normalization is an appropriate way of compensating for an imperfect mixing ratio.
In . Therefore, the mixing ratio  can be calculated as: The mixing ratios so determined were used for example to correct the values of all non-normalized ratios used for the plots in Supplementary Fig. S9. More importantly they were integrated in the calculations of kinetic parameters. We thus introduced the mixing ratio  and we defined the amounts of the proteins in the mixed sample as U, V, X and Y (i.e. in contrast to u,v,x,y, which represent the cellular concentrations). Without loss of generality, we can define the reference amount as the one from the control cells so that we then have in the mixed sample:

Cleavage of proteins into peptides
The next step in the model construction was accounting for the systematic discrepancy between ratios measured with K-vs. R-containing peptides. We assumed that the difference in quantitation originated from a different efficiency in the cleavage of proteins during the enzymatic digestion step, with the heavy and medium isotopes being cleaved less efficiently than the light ones.
We described the evolution of the protein amount during the cleavage into peptides with the following scheme: where U R a represents the amount of protein U that are still available for cleavage, U R represents the peptides containing the amino acid R (this is what is observed in mass spectra), and k L , k M and k H correspond to the rates of cleavage into peptides for the proteins containing Light, Medium or Heavy amino acids markers respectively.
For the initial conditions for this "cleavage evolution" we have the values present when the cells were harvested and lysed, for example: where s is used here for the "time" parameter for the system of protein cleavage (this time is independent from the time t of evolution of the proteins in the cells), and t 1 is the time when the cells where harvested and lysed.
Solving for the evolution of the system from the scheme (28), for example for U R gives: Similarly for the other proteins we get these relationships:

Experimentally measured quantities and their integration into the model
Mass spectrometry data yield the ratio of the abundances of the peptides containing different labels, which are then used to infer the identified proteins and their quantitative ratios. For each protein there were 6 ratios measured: Note how in these equations we have used the indices K and R to indicate that what is really measured in the mass spectra are the peptides and not directly the proteins. The ratios were measured for 3 time-points in the experiment.

Fitting the data to the model
Fitting of the data, included several steps. Equations (32) were measured at 3 times, so for each protein there were 18 experimental data points that could be used to fit our parameters. Only proteins having valid ratios measured at all time points were used for the calculations. The unknown parameters in the model were k A , k B , V A , V B for each protein, as well as the Q L , Q M and Q H that were global coefficients assumed to be the same for all proteins. It is important to note that in the equations, the V A or V B never appear isolated but only as ratios V B V A so it will only be possible to derive the ratio of the synthesis rates in the two conditions.
The flowchart describing the estimation of parameters is shown in Main Figure 3. In detail:

First step of fitting
In the first step of fitting, we fit simultaneously the data for 100 protein species, It is important to note that the corrections for the mixing ratios and for the cleavage rates for K,R peptides were very important to obtain good fitting of the pcSILAC ratios to the model as well as to avoid computational artifacts such as negative decay rates, which arose when these corrections were not done.

Derivation of V A , V B , p B
From the values obtained up to this point, it is already possible to calculate the expected ratio of concentrations of total proteins p B /p A , considering that for the drug treated cells at their steady state: And therefore, if p A , p B are the protein concentrations at steady state in the control or drug cultures, we have: Since V B /V A as well as k A,app and k B,app are known, p B /p A is fully determined. Furthermore, it is also possible to exploit the values of p A determined using iBAQ scores to obtain more values: and then by substitution: And finally It is important to note that p B is the calculated protein concentration in the treated cells at steady state, i.e. with t=infinite. In conclusion, pcSILAC data allow the determination of all kinetic parameters, although in the present case the values of synthesis rates and protein concentrations remain in arbitrary units because of the lack of an internal standard allowing the calibration of intensities to obtain absolute concentrations.

Calculation of fluxes
To assess quantitatively the instantaneous changes in the treated sample we calculated fluxes due to synthesis, degradation and the dilution effect. Let's recall that the mass balance equation is: All the 3 terms on the right side can be considered as fluxes W: and therefore On the other hand, the drug treated system will need some time to approach equilibrium. At each time point it is thus possible to estimate the contribution of the terms of degradation and synthesis to the total change. dp B dt t But since at later time points p B has been influenced by other terms such as dilution and synthesis, the contributions of the individual fluxes cannot be clearly separated anymore.

Correction for growth rate
We have to assume that cells are growing exponentially and that the growth rate does not change during the experiment. In this case, the cell growth can be expressed by: where N(t) is the biomass measured in million cells per ml and t and t 0 are the times at which the cells were counted.
Thus the correction term can be calculated as

Calculation of half-lives
Assuming a first order exponential decay, the half-life of a species can be calculated very simply (for example for u) as follows:

Supplementary Figure S2 : Samples and conditions of standard and pcSILAC experiments
Experimental design for pcSILAC experiments; the yellow arrow indicates treatment with geldanamycin (GA) or DMSO. Correspondingly numbered replicates were mixed (e.g. M1+H1, M2+H2) after total protein quantification. For replicates M1 and H1 the treatment was inverted, i.e. H1 received DMSO while M1 was treated with GA. A) experimental design for pcSILAC experiment 1 ; coloured triangles indicate medium exchange. The samples used as an internal standard SILAC experiment in pcSILAC experiment 1 were derived from the same culture used for pcSILAC labeling, therefore no medium exchange was done and quantitation was performed only on R (H/M ratio). B) same as A) for pcSILAC experiment 2 (no internal standard SILAC replicate was done in this experiment).  Figure S3 : predicted and possible isotope peaks for various classes of tryptic peptides in pcSILAC A) Isotope labelled amino acids used in the study and their monoisotopic mass shifts relative to light amino acids B) General Labelling scheme for pcSILAC with description of the isotopomers present in the media before/after medium exchange (represented as a coloured triangle) at t=0 C) possible isotope peaks of 1 x R-or 1 x K-containing peptides after mixing H+L pcSILAC samples, with expected mass shifts relative to a R0, resp. K0 peptide. D) Same as C) but for peptides containing 2 K residues. The occurrence of mixed-label peptides (*) during the phase of medium exchange was considered. E) Same as D) but for peptides containing 2x R residues. F) Same as D) but for peptides containing both    Fig. S5 : peptides containing both R and K in pcSILAC A) Isotope pattern observed for peptide GHYTEGAELVDSVLDVVRK (TBB2A_HUMAN) at t=6h in the H+L mix of pcSILAC experiment 1, replicate 2. The position of possible peaks resulting from label mixing (old label + new label) are indicated below the axis, with their expected mass shifts (see Supplementary figure 2) B) Same as A) but for peptide GVAINMVTEEDKR (IF4A1_HUMAN) at t=12h in pcSILAC experiment 1, replicate 1. Overall, the predicted isotope patterns for K+R-peptides were observed. Peaks resulting from label mixing (old label + new label) were sometimes detectable but had very low intensity (<3% of base peak, marked below x-axis) compared to homogeneously (old+old, new+new) labelled species. Spectra obtained for samples before mixing were also inspected with similar results.    A) correlation between theoretical S 1 and S 2 values (from pcSILAC data) obtained with two complementary equations (eq.5 and eq.6 in main manuscript ) using different ratios. Data were from replicate 2 of experiment 1, not corrected for mixing ratio inequalities B) Correlation of S (calculated as average of S 1 , S 2 , normalized ) with experimental standard SILAC (net protein) ratios from a parallel experiment performed simultaneously (« internal » control) or C) an independent standard SILAC experiment performed 6 weeks earlier.  Values are after correction for mixing inequality and log2 transformation. Some proteins with apparent multiphasic behavior are shown in red, together with their value of the ResNorm parameter (Blue), which is a measure of the quality of the fitting. Most proteins highlighted here have values of ResNorm well above the median (0.18) of ResNorm for the whole population. It is to be noted however that, since the fitting is performed on 6 series of ratios, deviations due to noise in one series can be compensated by other values. DNAJB1 and GRP78 were quantitated with large numbers of peptides and thus their values should be reliable and reflect a complex temporal dynamic of changes.