Comparison of β-Casein and β-Lactoglobulin Proteolysis by Trypsin: Quantification of Demasking and Hydrolysis of Peptide Bonds

1. Introduction

Enzymatic hydrolysis of proteins by proteases (proteolysis) is widespread in living systems and is also used in the food industry, biotechnology, pharmaceuticals, cosmetology etc. Proteolysis is quite difficult to study, since it is a continuous process of interconversion of various peptide fragments, i.e. the composition of the reaction mixture changes significantly during the hydrolysis of peptide bonds. The concentrations of the final products of proteolysis and of sufficiently long intermediate peptides containing not yet hydrolyzed enzyme-specific bonds change with the time of hydrolysis and it is difficult to predict these concentration dependences [1,2,3].

A practically important task is to optimize the production of biopeptides from proteins using proteolysis. To obtain the final products of proteolysis, proteolysis must be carried out to the end. In contrast, intermediate biopeptides can be obtained by limited proteolysis during a given time interval, after which the process should be stopped to avoid further hydrolysis of peptide bonds. To determine this time interval, or the corresponding degree of hydrolysis of peptide bonds, it is necessary to use proteolysis modeling to assist the experimental path, which requires a large number of experiments [4,5,6,7,8,9,10,11,12].

Using the rate constants of hydrolysis of peptide bonds and other kinetic parameters, the time course of final and intermediate products can be calculated using the methods of chemical kinetics [4,9,10]. To simplify this task, it is proposed to use several approaches that both reduce the number of rate constants and the number of peptide fragments that are taken into consideration when modeling the proteolysis process [9].

Unlike low-molecular protease substrates, many peptide bonds in proteins are initially inaccessible to the enzyme, since they are located inside the protein globule or, for example, in the casein micelle. The destruction of the protein globule or micelle at the first stage of proteolysis initiates the process of opening the internal peptide bonds for a successful enzymatic attack [1,4,13]. This process, called demasking, eliminates steric hindrances for the enzyme so that the active center of the enzyme and the amino acid residues of the peptide substrate can fully interact. After demasking, hydrolysis of various peptide bonds occurs with the different hydrolysis rate constants. At this stage of proteolysis, the peptide bonds are hydrolyzed in accordance with their specificity, determined by the amino acid sequence of the polypeptide chain of the substrate.

Previously, we described methods for determining the demasking rate constants using fluorescence spectroscopy, which allows one to estimate the rate of substrate opening during proteolysis by the red shift of fluorescence accompanying proteolysis [13]. The use of this spectral method in addition to the traditional determination of the rate of hydrolysis of peptide bonds allowed us to analyze the interdependence of the demasking and hydrolysis processes within the framework of the two-step proteolysis model [4,13]. The demasking rate constants obtained by the fluorescence method were used to model the release of peptides during proteolysis of β-lactoglobulin (β-LG) by trypsin [9].

Modeling of peptide release during proteolysis of β-LG by trypsin was recently completed using a proteolysis model that takes into account both one- and two-stage demasking of peptide bonds without their secondary masking [9]. In addition to simplifications regarding the proteolysis mechanism, some simplifications of the calculation procedure were also made. It was assumed that there is no secondary masking [14], the size of the intermediate peptide fragments is limited by two enzyme-specific peptide bonds within them, which can be subjected to further hydrolysis. The formation of some minor peptide fragments, which are formed very slowly, was neglected [9].

It is believed that the hydrolysis of proteins depends on their spatial structure, at least at the initial stage of the process [1,15,16]. Traditionally, proteolysis of globular milk proteins, including β-LG, and structurally disordered caseins, most often β-casein (β-CN), is compared to reveal these dependencies [1,14,17].

The major protein in bovine whey is the small globular protein β-LG (55–60%) [18,19]. The monomeric β-LG is made up of 162 amino acid residues (∼18.3 kDa) and is stabilized by 2 disulfide bonds [19]. Its secondary structure was predicted mainly as β-sheet (50%) [20]. β-LG and its enzymatic hydrolysates have high functional and nutritional properties. The hydrolysates of β-LG contain a lot of biopeptides that can be used as biologically active additives in therapeutic nutrition and cosmetic [21].

Casein micelles are colloidal complexes of proteins and salts, and their main biological function is to transport sparingly soluble calcium phosphate in liquid form to infants [22]. They have four major protein species, termed αs1-, αs2-, β- and κ-caseins [22]. β-CN is the most hydrophobic one within the group of caseins, which constitutes about 45% of the casein of bovine milk [18]. It is a ∼24 kDa single polypeptide chain which consists of 209 amino acid residues [23]. Since β-CN does not contain disulfide bonds, it has no tertiary structure, but some regions of the polypeptide chain have secondary structure. β-CN has a hydrophobic C-terminus and a hydrophilic negative N-terminal region [24]. Caseins are known to be easily hydrolysed by proteases due to their conformational flexibility and the abundance of the enzyme-accessible peptide bonds. β-CN is of interest because of its nutritional importance and utility as a drug delivery vehicle [25,26]. The hydrolysates of β-CN, as well as β-LG, are rich in biologically active peptides [27,28].

Since demasking of peptide bonds predetermines their hydrolysis, it can be assumed that the demasking rate constants can be calculated from the distribution of the resulting peptide fragments. Previously, the sum of the concentrations of the peptides obtained during the hydrolysis of the analyzed peptide bond was used to determine the selectivity parameter of this bond in the hydrolyzed protein [29,30,31]. Herein, we expand this approach and present a method for estimating the demasking rate constants by analyzing changes in the sums of the concentrations of some peptides over time. Grouping of such peptides and calculation of the sum of their concentrations corresponds to the proposed model of proteolysis [9].

Computer search for bioactive peptides in enzymatic hydrolysates of proteins of various origins has recently become increasingly popular, requiring the development of reliable algorithms for in silico proteolysis [32,33,34]. Existing approaches are primarily aimed at predicting peptide fragments that can potentially be released during enzymatic hydrolysis of various protein substrates. Prediction of concentration dependences for peptide fragments and the degree of hydrolysis required to obtain target peptides was not performed. These aspects of in silico proteolysis are of particular interest to us, as the proteolysis model we are developing is capable of solving such problems.

A simpler method for determining demasking rate constants than using fluorescence spectroscopy is described here for a proteolysis model formulated earlier [9]. Previously, demasking rate constants were estimated from the fluorescence shift during proteolysis, but now this is done using experimental data on peptide release. As a proteolysis substrate, along with β-LG, studied earlier [9], β-CN is also considered. In addition to a detailed demonstration of the method for calculating kinetic curves, our goal is to compare the results predicted by the model with experiment for these two substrates with different protein structures.

2. Results

2.1. Proteolysis model

The same scheme of polypeptide chain cleavage was used here as in the previous work devoted to modeling the release of peptides during proteolysis [9]. The main principles of the model are presented in Figure 1. The hydrolysis of masked peptide bonds is preceded by their two-stage demasking, associated with conformational changes in the polypeptide chains and leading to the accessibility of peptide bonds for the enzyme. The illustration shows the peptide fragment ABC formed by the hydrolysis of the most specific peptide bonds at the beginning of proteolysis. It consists of three relatively short peptides A, B and C, which include only peptide bonds that are nonspecific for the enzyme. These peptides are linked by enzyme-specific peptide bonds A-B and B-C, indexed i and j, respectively.

The left side of Figure 1 shows that these bonds are not hydrolyzed until their demasking is complete. The right side of Figure 1 shows how these demasked bonds are hydrolyzed to form the corresponding peptide fragments.

Thus, Figure 1 shows the transformation of the peptide chain -ABC- into a partially demasked fragment -ABC at the first stage of demasking. Then, at the second stage of demasking, a completely demasked peptide fragment ABC is obtained. After demasking, the fragmentation of the ABC trimer occurs as a result of hydrolysis of the peptide bonds i and j with the formation of peptide fragments AB, BC, A, B, and C. The model determines that hydrolysis of the most specific bonds occurs at the beginning of proteolysis together with demasking, and hydrolysis of less specific bonds is possible only after the completion of demasking.

2.2. Application of the fragmentation scheme to the proteolysis of β-CN by trypsin

To specify the general scheme (Figure 1) for the proteolysis of β-CN by trypsin, we collected data that were previously published for this enzyme-substrate pair (Table 1). These parameters include the selectivity and initial rate of hydrolysis of the enzyme-specific peptide bonds Lys-X and Arg-X in β-CN. Both of these characteristics correlate fairly well with each other, and they allow ranking various peptide bonds by their rate of hydrolysis. From these data it follows that peptide bonds with indexes 99, 105, and 169 are the most rapidly hydrolyzed bonds, while bonds 1, 25, 32, 48, and 202 are the most slowly hydrolyzed bonds, the hydrolysis of which is neglected within the framework of the considered model.

The remaining bonds (28/29, 97, 107, 113, 176, 183) are considered initially masked, and the model must quantitatively determine how these peptide bonds are demasked and hydrolyzed. Thus, during proteolysis of β-CN by trypsin, the ABC fragments are f(1-99), f(106-169), and f(170-209), and the fragment f(100-105) is a rapidly forming peptide that does not contain enzyme-specific bonds.

Peptide bonds assigned to the group of the most slowly hydrolyzed bonds are considered non-hydrolyzed in our simplified approximation, therefore the concentrations of peptide fragments obtained during their hydrolysis are added to the concentrations of the peptides from which they originated. For example, the concentrations of fragments 184-202 and 202-209 are also added to the concentration of peptide 184-209, since bond 202 is considered the most slowly hydrolyzed bond (Table 1). The sum obtained in this way is considered the concentration of fragment 184-202.

In the previous analysis of β-LG proteolysis by trypsin, demasking of two ABC fragments occurred by the two-stage demasking mechanism, and of one trimer by the one-stage demasking mechanism [9]. Proteolysis of β-CN by trypsin is a multistep process, including association of intermediate peptides and formation of nanoparticles during proteolysis [35]. Due to the complexity of this process, we consider demasking all three trimers f(1-99), f(106-169) and f(170-209) using the two-stage demasking scheme.

2.3. Estimation of rate constants for demasking and hydrolysis

Like the selectivity parameter [29], $k_d^f$was determined using a simple exponential equation for the most rapidly hydrolyzed peptide bond (Equation (9)).

In this work, a new approach to estimate the demasking rate constants was used, different from the fluorescence spectroscopy method [13]. This approach requires additional calculations using the dependences of the sum of concentrations of some peptide fragments on the hydrolysis time. To calculate the demasking rate constant k_d the equation following from the material balance equations was used. This Equation (10) is given in Section 4 for the sum of the concentrations [ABC], [AB], and [A], corresponding to balance of A units. To track balance of C units, it is necessary to calculate the sum [ABC]+[BC]+[C]. As a result of fitting according to Equation (10), k_d value was obtained for each of ABC trimer. The found values of k_d then allowed us to determine the rate constants of hydrolysis of peptide bonds using the Equation (11) we had previously derived for the concentrations of peptide bonds [1,14].

A new approach to estimate the rate constants required to model the peptide release according to Scheme 1 was applied to the proteolysis of β-CN by trypsin (Table 2). For comparison with the previous result on peptide release [9], this method was also applied to the proteolysis of β-LG by trypsin (Table 2).

2.4. Simulation of peptide release during proteolysis by trypsin

Figure 2a shows an example of the calculated dependences of peptide concentrations on the hydrolysis time for the intermediate and final peptide fragments. The intermediate peptide products (ABC, AB, and BC) are first formed and then disappear due to the hydrolysis of internal enzyme-specific peptide bonds (Figure 2a). The final products (A, B, and C) only accumulate, since they do not contain internal enzyme-specific peptide bonds. Demasking of the peptide bonds leads to the fact that the concentrations of proteolysis products may increase not immediately with the onset of proteolysis, but after a lag phase [1].

The concentrations of peptide fragments can also be represented as functions of the degree of hydrolysis achieved by a given time (Figure 2b), since the dependence of the degree of hydrolysis on time is known. Representation of concentration curves as functions of the degree of hydrolysis is more convenient for determining the mechanisms by which peptides are released during proteolysis. The transition from time to the degree of hydrolysis does not change the shape of the curves for intermediate products, but changes the dependences for the final peptides. For them, the curves for rapidly released peptides remain convex, and the curves for slowly released peptides become concave.

The parameter d_r, the degree of hydrolysis at which the major part of a given peptide is released, was calculated for each of the intermediate products using the following simple equation:

$$d_r={\displaystyle \sum _{i=1}^6{d}_i}\times C(d_i)/{\displaystyle \sum _{i=1}^{6}C(d_i})$$

(1)

where C(d_i) are the concentrations of peptide fragments determined at six hydrolysis degrees d_i (0, 1.5, 3, 4.5, 6, 7.9% for β-LG and 0, 1, 2, 3, 4, 4.3 for β-CN). Thus, d_r is defined as the weighted average of 6 degrees of hydrolysis, and the relative concentrations of intermediate peptides corresponding to these degrees of hydrolysis are the weighting factors.

For the final peptides, the parameter n, which characterizes their concentration dependence on the degree of hydrolysis C(d_i), was calculated using the following equation:

$$C(d_i)=a{\left(d_i/d_{\max }\right)}^n$$

(2)

where a is a constant factor, d_max is the maximal hydrolysis degree, and n is the exponent of the power function. Thus, we compared the kinetic curves for the final products simply by comparing the parameter n, which reflects the curvature of the curves.

Examples of the calculated and experimental concentration dependences for intermediate (Figure 3a) and final (Figure 3b) peptide fragments are given. Figure 3a shows the release of the intermediate trimeric peptide ABC during proteolysis and the subsequent formation of the intermediate dimeric peptide AB from it. Theoretically, one would expect that the d_r for the ABC fragment is less than that for the AB or BC fragment. Figure 3b shows the concentration curves for the final peptides. For a rapidly forming final peptides the curve is convex and n is less than 1, while for a slowly forming final peptides the curve is concave and n should be greater than 1.

Table 3 presents the parameters of the concentration curves for all studied peptide fragments, namely d_r for the intermediate peptides and n for the final peptides. The parameters d_r and n were calculated using Equations (1, 2) from the simulated curves obtained using the rate constants from Table 2 and the experimental concentration curves. The published peptide concentrations were for tryptic proteolysis of β-LG and β-CN, and the maximal hydrolysis degree (d_max) was 7.9 and 4.3% for β-LG and β-CN, respectively [7]. The concentrations of the peptide fragments were calculated using Equations (3-8) for the same degrees of hydrolysis as in the experiment.

A comparison of the n parameter for the simulated and experimental concentration curves obtained for both β-LG and β-CN proteolysis is shown in Figure 4a. The coefficient of proportionality between the experimental and calculated n values was 0.999 ± 0.038 (r²=0.9998) for β-LG (Figure 4a) with the expected coefficient of 1. Thus, the agreement between the simulation results and the experiment was very good when using the n parameter for comparison. The proportionality coefficient for the d_r values was 0.68 ± 0.07 (r²=0.87), indicating an underestimation of this parameter in the simulation for intermediate peptide fragments with high d_r.

For the parameter n, the proportionality coefficient was 0.77±0.10 (r²=0.892) with previous quantification of rate constants [9], while in this study it is 0.96±0.04 (r²=0.989) for all points in Figure 4a. Thus, the use of a new method for estimating the rate constants of demasking and hydrolysis has led to a more accurate prediction of the release of proteolysis products using the same proteolysis model that we proposed earlier [9]. This confirms the correctness of the main ideas of the model based on taking into account the processes of the peptide bond demasking before their hydrolysis and the difference in hydrolysis rate for the various peptide bonds.

Despite the differences in the protein structure of the substrates, calculations of the concentration dependences and parameters of d_r and n were performed using the same algorithm and yielded results consistent with the experiment.

2.5. Prospects for in silico proteolysis

The possibility of predicting the hydrolysis constants kⁱ for various protein substrates based on the analysis of their amino acid sequences is illustrated by Figure 4b. The total number of charged amino acid residues in the positions P₂ and P_-2 was chosen as the determining parameter. It was previously shown that the presence of charged residues (E, D, R, K) in these positions significantly reduces the selectivity parameters for proteolysis with trypsin [7]. Figure 4b shows that lg(kⁱ) decreases with the number of charged groups in the positions P₂ and P_-2. To construct more reliable correlation equations, a larger volume of experimental data obtained for a number of protein substrates is required. It is also possible to take into account the amino acid residues in the positions P₁ and P_-1, as was done to predict the Michaelis parameters for the hydrolysis of peptide bonds in low-molecular substrates. Within the framework of our approach, the structure-reactivity correlation equations must be sought specifically for the hydrolysis rate constants kⁱ, since these parameters reflect the hydrolysis step of proteolysis. According to the model, all demasked regions of the polypeptide chain and their corresponding kⁱ, regardless of their affiliation with specific protein substrates, must be considered together to find the correlation equation corresponding to a given proteolytic enzyme. Therefore, we plotted the data for β-LG and β-CN on a single graph (Fig. 4b).

The demasking rate constant $k_d^f$=0.46 for β-LG, determined by the opening of protein globule, is expectedly lower than that of β-CN ($k_d^f$=3.5), which is a disordered coil of the polypeptide chain with more accessible peptide bonds. Consequently, the ratio of the demasking rate constants of β-CN and β-LG determined here is 7.6 under the proteolysis conditions published in [7]. It makes sense to compare proteolysis processes carried out at the same concentrations of substrate and enzyme, if the temperature and other environmental conditions remain unchanged. Previously, the ratio of the demasking rate constants during proteolysis of these two substrates by trypsin was determined from the fluorescence shift [13], and the value kII(β-CN)/kII(β-LG) = 0.0244/0.0189/0.1 = 12.9 was obtained [13]. The proteolysis conditions in that study differed from those used in the present calculation, which may explain the difference in the estimates of this ratio. The secondary demasking rate constants k_d are apparently more complexly related to the structure of the polypeptide fragments. For β-LG proteolysis, the lowest k_d = 0.32 was observed for the f(76-100/101) fragment, which participates in the formation of the core resistant to trypsin hydrolysis [36]. For β-CN proteolysis, the lowest k_d = 1 was determined for the central part and the hydrophobic C-terminus of β-CN f(106-209). A significant contribution of the hydrophobic C-terminal fragment to the demasking process was previously noted, explaining the formation of C-terminal peptides of β-CN with a lag [4].

The dependence of the model parameters on enzyme and substrate concentrations is essential for the comparison and correct interpretation of experiments performed under various kinetic conditions. We have demonstrated how the contribution of secondary masking increases during proteolysis and how the concentration dependences change with decreasing enzyme concentration at a constant substrate concentration [14]. However, experimental data are not yet sufficient to accurately determine the dependence of the model parameters on the conditions of proteolysis, mainly on the substrate concentration.

3. Discussion

The algorithm for calculating the release of peptides during proteolysis that we used is a new method that considers the need to demask peptide bonds before their hydrolysis [1,9]. The phenomenon of peptide bond demasking was first discovered based on the analysis of non-monotonic changes in the rate of hydrolysis of peptide bonds during enzymatic hydrolysis of milk protein substrates [4]. Then, the kinetics of demasking was studied by monitoring the fluorescence of tryptophan residues in the comparative analysis of β-LG and β-CN proteolysis by trypsin [13,14,17]. The importance of considering demasking processes in studying proteolysis was shown by other analytical and physicochemical methods [37,38].

It was shown previously that the precise determination of the demasking parameters may be even more important than the precise determination of the rate constants of peptide bond hydrolysis [9]. In the previous work, the demasking rate constants were determined using the spectral method, which itself could yield additional errors. In order to exclude possible inaccuracies in determining the demasking rate constants, in this work the demasking parameters were determined precisely using the published experimental data [7] on the distribution of peptide fragments during proteolysis. Here, more precise values of k_d were used for each of the ABC trimer under consideration, but the calculation algorithm of fragment concentrations was the same. The results of modeling the release of peptides gave better agreement with the experiment than in the previous work [9].

The focus here is on the calculation algorithm. It turned out that the simulation results using the same equations for the proteolysis of β-LG and β-CN proteins provide good agreement with experiment. Thus, the proposed calculation algorithm with more precise parameters correctly accounts for the key features of proteolysis even for structurally different protein substrates.

According to AFM, FTIR spectroscopy and static light scattering data, proteolysis of β-CN by trypsin is accompanied not only by degradation of the initial β-CN micelles, but also by parallel formation of new nanoparticles, which then degrade to form peptide proteolysis products [1,35,39]. When viewed at the macromolecular level, there is a parallel degradation and association of intermediate polypeptide products with degradation gradually becoming predominant as the process progresses. New nanoparticles formed during proteolysis are denser than the initial micelles, and the content of secondary protein α-, β-structures in them does not decrease during proteolysis, but is comparable to that in the initial β-CN micelles. In this regard, during proteolysis at low enzyme concentrations, the content of secondary protein structures changes non-monotonically. It first decreases due to hydrolysis of the original micelles, then increases due to the formation of these new nanoparticles, and finally slowly decreases as they are hydrolyzed [35,39]. The kinetics of degradation of β-CN micelles by trypsin was described by analytical functions in the case of modeling the proteolytic system using linear differential equations [35]. This process was also described by numerical modeling methods using nonlinear differential equations [1]. In this case, the original micelles, nanoparticles of different sizes and their aggregates, as well as the change in the concentration of the active enzyme during proteolysis were considered [1].

The macromolecular description of β-CN proteolysis by trypsin [35,39] and the model considered here represent sequential three-stage processes. A fairly rapid hydrolysis of micelles is the first stage in the macromolecular model and it is initiated by the hydrolysis of the most specific peptide bonds [35]. Hydrolysis of the most specific bonds also initiates the process of the first stage of demasking with rate constant $k_d^f$. The ratio of the rate constant of micelle hydrolysis to the rate constant of peptide product formation in the macromolecular model was 33:1 [35]. A close value of 35 was obtained for the ratio of $k_d^f$to k^j for the peptide bonds with j=107 and 113.

Relatively small β-CN and β-LG are used as model reactions of proteolysis with various proteases, allowing the kinetic features of proteolysis to be studied [1]. The number of peptide fragments is not very large (several dozen) and they are easily identified and quantified by HPLC-MS [40,41,42,43,44,45]. Among these studies, those in which the concentrations of peptide fragments are determined at several time points during proteolysis are especially important. Such data are suitable for verifying proteolysis models and determining the numerical values of kinetic parameters. In particular, such studies have been carried out for trypsin proteolysis of β-CN [7] and β-LG [7,8,10]. Accurate determinations of the concentrations of individual peptides made it possible to determine the change in the concentrations of peptide bonds over the course of proteolysis and then determine the selectivity parameters for various peptide bonds [7]. In this study, summing the concentrations of individual peptides ABC, AB, and A was used to determine the demasking rate constants.

Although spectroscopic methods were not used in this work, these methods have great potential for the monitoring of protein opening during proteolysis and hence for the monitoring of demasking. The possibility of the monitoring of protein structure changes in course of proteolysis was demonstrated using fluorescence [13], infrared (FTIR) [35], ultrasound [46,47,48] and other physicochemical methods. Undoubtedly, further development of these methods is of great importance for successful modeling of proteolysis.

In practical-oriented proteolysis models, the peptide bonds and their hydrolysis rate constants are assumed to be identical, and the main model variables are the degree of peptide bond hydrolysis and the total hydrolysis rate [49,50,51,52,53,54,55]. Such models can be called total proteolysis models, since they do not consider the hydrolysis processes of different hydrolysis sites (different peptide bonds). The model we are developing is more accurate because it considers the differences in the hydrolysis of various peptide bonds in proteins, as well as the limited accessibility of the enzyme to some temporarily masked regions of the polypeptide chain [1,9]. To optimize the production of protein hydrolysate with the required degree of hydrolysis, it is necessary to use total proteolysis models, while to optimize the yield of target biopeptides, it is necessary to use our approach.

Although the relationship between the structural features of β-LG and β-CN and the demasking parameters is demonstrated in this study, it is still difficult to formulate specific correlation dependences that could be used for the quantitative prediction. Since in silico proteolysis of various protein substrates is impossible without successful prediction of model parameters, a hybrid approach can be proposed in which some parameters are determined experimentally. Thus, in the calculations presented here, the demasking parameter d_r was determined by analyzing experimental data with Equation (10), and the identification of rapidly hydrolyzed bonds was also based on the experimental data (Table 1). Apparently, the success of hybrid approach is possible if the number of experimentally determined parameters is minimized. We hope that the use of our approach will bridge the gap between the experimental study of specific proteolytic processes and the development of a general model of proteolysis.

4. Materials and Methods

4.1. Quantitative modelling of proteolysis

The solutions of the differential equations corresponding to the proteolysis scheme (Fig. 1) were found in the form of analytical functions of time t, min. The relative concentrations C(t) of the peptide fragments ABC, AB, BC, A, B, and C are:

$$[ABC]={\displaystyle \frac{k_d^f{k}_d}{\left(k_d-k_d^f\right)\left(k^i+k^j-k_d^f\right)}}{e}^{-k_d^{f}t}+{\displaystyle \frac{k_d^f{k}_d}{\left(k^i+k^j-k_d^f\right)\left(k^i+k^j-k_d\right)}}{e}^{-(k^i+k^j)t}-{\displaystyle \frac{k_d^f{k}_d}{\left(k_d-k_d^f\right)\left(k^i+k^j-k_d\right)}}{e}^{-k_{d}t}$$

(3)

$$\begin{array}{l}[AB]={\displaystyle \frac{k_d^f{k}_d{k}^j}{\left(k_d-k_d^f\right)\left(k^i+k^j-k_d^f\right)\left(k^i-k_d^f\right)}}{e}^{-k_d^{f}t}-{\displaystyle \frac{k_d^f{k}_d}{\left(k^i+k^j-k_d^f\right)\left(k^i+k^j-k_d\right)}}{e}^{-\left(k^i+k^j\right)t}-{\displaystyle \frac{k_d^f{k}_d{k}^j}{\left(k_d-k_d^f\right)\left(k^i+k^j-k_d\right)\left(k^i-k_d\right)}}{e}^{-k_{d}t}+\\ {\displaystyle \frac{k_d^f{k}_d\left[\left(k^i+k^j\right)\left(k^i+k^j-k_d^f-k_d\right)+k_d^f{k}_d\right]}{\left(k^i-k_d^f\right)\left(k^i-k_d\right)\left(k^i+k^j-k_d^f\right)\left(k^i+k^j-k_d\right)}}{e}^{-k^{i}t}\end{array}$$

(4)

$$\begin{array}{l}[BC]={\displaystyle \frac{k_d^f{k}_d{k}^i}{\left(k_d-k_d^f\right)\left(k^i+k^j-k_d^f\right)\left(k^j-k_d^f\right)}}{e}^{-k_d^{f}t}-{\displaystyle \frac{k_d^f{k}_d}{\left(k^i+k^j-k_d^f\right)\left(k^i+k^j-k_d\right)}}{e}^{-\left(k^i+k^j\right)t}-{\displaystyle \frac{k_d^f{k}_d{k}^i}{\left(k_d-k_d^f\right)\left(k^i+k^j-k_d\right)\left(k^j-k_d\right)}}{e}^{-k^{j}t}\\ +{\displaystyle \frac{k_d^f{k}_d\left[\left(k^i+k^j\right)\left(k^i+k^j-k_d^f-k_d\right)+k_d^f{k}_d\right]}{\left(k^j-k_d^f\right)\left(k^j-k_d\right)\left(k^i+k^j-k_d^f\right)\left(k^i+k^j-k_d\right)}}{e}^{-k^{j}t}\end{array}$$

(5)

$$[A]=1-{\displaystyle \frac{k^i{k}_d}{\left(k_d-k_d^f\right)\left(k^i-k_d^f\right)}}{e}^{-k_d^{f}t}+{\displaystyle \frac{k_d^f{k}^i}{\left(k_d-k_d^f\right)\left(k^i-k_d\right)}}{e}^{-k_{d}t}-{\displaystyle \frac{k_d^f{k}_d\left[\left(k^i+k^j\right)\left(k^i+k^j-k_d^f-k_d\right)+k_d^f{k}_d\right]}{\left(k^i-k_d^f\right)\left(k^i-k_d\right)\left(k^i+k^j-k_d^f\right)\left(k^i+k^j-k_d\right)}}{e}^{-k^{i}t}$$

(6)

$$\begin{array}{l}[B]=1-{\displaystyle \frac{k_d{k}^i{k}^j(k^i+k^j-2k_d^f)}{\left(k_d-k_d^f\right)\left(k^i+k^j-k_d^f\right)\left(k^j-k_d^f\right)\left(k^i-k_d^f\right)}}{e}^{-k_d^{f}t}+{\displaystyle \frac{k_d^f{k}_d}{\left(k^i+k^j-k_d^f\right)\left(k^i+k^j-k_d\right)}}{e}^{-\left(k^i+k^j\right)t}\\ +{\displaystyle \frac{k_d^f{k}^i{k}^j(k^i+k^j-2k_d)}{\left(k_d-k_d^f\right)\left(k^i+k^j-k_d\right)\left(k^j-k_d\right)\left(k^i-k_d\right)}}{e}^{-k_{d}t}-{\displaystyle \frac{k_d^f{k}_d\left[\left(k^i+k^j\right)\left(k^i+k^j-k_d^f-k_d\right)+k_d^f{k}_d\right]}{\left(k^i-k_d^f\right)\left(k^i-k_d\right)\left(k^i+k^j-k_d^f\right)\left(k^i+k^j-k_d\right)}}{e}^{-k^{i}t}\\ -{\displaystyle \frac{k_d^f{k}_d\left[\left(k^i+k^j\right)\left(k^i+k^j-k_d^f-k_d\right)+k_d^f{k}_d\right]}{\left(k^j-k_d^f\right)\left(k^j-k_d\right)\left(k^i+k^j-k_d^f\right)\left(k^i+k^j-k_d\right)}}{e}^{-k^{j}t}\end{array}$$

(7)

$$[C]=1-{\displaystyle \frac{k^j{k}_d}{\left(k_d-k_d^f\right)\left(k^j-k_d^f\right)}}{e}^{-k_d^{f}t}+{\displaystyle \frac{k_d^f{k}^j}{\left(k_d-k_d^f\right)\left(k^j-k_d\right)}}{e}^{-k_{d}t}-{\displaystyle \frac{k_d^f{k}_d\left[\left(k^i+k^j\right)\left(k^i+k^j-k_d^f-k_d\right)+k_d^f{k}_d\right]}{\left(k^j-k_d^f\right)\left(k^j-k_d\right)\left(k^i+k^j-k_d^f\right)\left(k^i+k^j-k_d\right)}}{e}^{-k^{j}t}$$

(8)

In the calculations, $k_d^f$were 0.46 min⁻¹ [9] for β-LG and 3.5 min⁻¹ for β-CN, the other rate constants (k_d and kⁱ ) were taken from Table 2.

4.2. Estimation of the rate constants

The parameter$k_d^f$ was the largest value of the rate constant for hydrolysis of various peptide bonds corresponding to the following simple equation for dependence of the peptide bond concentration on the hydrolysis time:

$$N^j(t)=N_0^j\left(1-e^{-k^{j}t}\right)$$

(9)

By summing the concentrations of components ABC, AB and A from Equations (3), (4) and (6), we obtained the following equation for the sum [ABC]+[AB}+[A]:

$$[ABC]+[AB]+[A]=A_0\left[1-{\displaystyle \frac{k_d{e}^{-k_d^{f}t}}{(k_d-k_d^f)}}+{\displaystyle \frac{k_d^f{e}^{-k_{d}t}}{(k_d-k_d^f)}}\right]$$

(10)

A similar equation was obtained for the sum [ABC]+[BC}+[C]. The dependences of the sums of the components [ABC]+[AB}+[A] and [ABC]+[BC}+[C] on the hydrolysis time were calculated from published data [7] and used to estimate the parameter k_d with Equation (10).

The parameters $k_d^f$and k_d were then used to calculate the rate constants for the hydrolysis of different peptide bonds using the following equation expressing the dependence of the peptide bond concentration on the hydrolysis time [1,14]:

$$N^j(t)=N_0^j\left[1-{\displaystyle \frac{k_d{k}^j{e}^{-k_d^{f}t}}{(k_d-k_d^f)(k^j-k_d^f)}}-{\displaystyle \frac{k_d^f{k}^j{e}^{-k_{d}t}}{(k_d-k_d^f)(k_d-k^j)}}-{\displaystyle \frac{k_d^f{k}_d{e}^{-k^{j}t}}{(k^j-k_d^f)(k^j-k_d)}}\right]$$

(11)

The experimental concentrations of peptide fragments taken from Table S-4 [7] were divided by a factor of 50 to obtain the relative concentrations of the peptide fragments that were used in the calculations. The values of model parameters d_r and n were independent of this factor.