Prediction of protofilament structures

CABS-dock standalone [25] version 0.9.18 was used for the test prediction of six distinct protofilament structures composed of in-register parallel peptide chains. The following experimental structures were extracted from the Protein Data Bank database, PDB IDs: 1) 2E8D: a beta2-microglobulin fragment probed by solid-state NMR [53], 2) 2NNT: a structural model of the CA150.WW2 protofilament constructed based on magic-angle-spinning (MAS) NMR spectroscopy, alanine scanning and EM experiments [54], 3) 6TI5: a structural model of Aβ40 fibrils derived from solid-state NMR (SSNMR) experiments [55,56], 4) 6ZRQ: a protofilament amyloid structure of the S20G variant of human amylin (IAPP, islet amyloid polypeptide) solved using cryo-electron microscopy (cryo-EM) [57], 5) 7YAT: a protofilament structure of the hamster prion peptide (sHaPrP, sequence 108–144) solved by cryogenic electron microscopy [58] and 6) 7Q66: a shorter protofilament fragment of the FG-repeat domain of human nucleoporin 98 (Nup98) obtained using a combination of NMR spectroscopy and cryo-electron microscopy [59]. For a detailed description of experimental structures see Table 1.

The predicted protofilament structures were generated during a docking simulation followed by a scoring procedure and selection of the best models (the CG docking and scoring scheme used for predicting protofilament structures is shown in Fig 1). During each docking simulation, five identical peptide monomers (copies of a single peptide chain) were simultaneously docked to each other. According to the CABS-dock input convention, the first peptide chain was treated as a receptor molecule and the four other monomer chains were treated as docked peptide ligands. Weak distance restraints were imposed between the side-chain (SC) atom pairs of the corresponding amino acid residues in adjacent peptide monomers in the predicted protofilament model. The restraint distance was set to 5 Å which corresponded to the average distance observed in the available protofilament experimental structures. The docking simulation started with a random position and random conformation of each peptide chain. All peptides were fully flexible during the docking calculations. Each simulation was performed using the same CABS-dock input parameters except for the amino acid sequence of each peptide (see supplementary materials, S1 Table). During a single docking simulation, 10,000 protofilament models were created and 1,000 lowest energy structures were selected. For each protofilament system, 40 independent docking simulations were completed, resulting in a set of 40,000 low-energy structures.

The selection of the best models was completed according to peptide interaction energy values estimated using CABS statistical potentials or the peptide-chain-average-RMSD (pcaRMSD) parameter or a combination of both. The pcaRMSD value was calculated for each predicted protofilament model using a set of RMSD values obtained during the mutual comparison of all five monomer chains (after structural fitting of two peptide chains for every pair). In this way, the conformations of the five peptides were compared with each other in a given protofilament model. In our earlier work, we showed that low pcaRMSD values correlate with a high level of translational symmetry of monomers within the protofilament structure [19]. The pcaRMSD values were calculated for all predicted models using the following formula:

Eq 1. The pcaRMSD formula used for scoring the predicted protofilament models.

where M is the number of peptide chains in the oligomer model, N is the number of Cα atoms, X_k is the coordinate vector for the k target Cα atom, Y_k is the coordinate vector for the k reference Cα atom, whereas i and j indicate indexes of a particular pair of compared peptide chains in the oligomer model.

Identification of enzymatic cleavage sites

CABS-dock docking simulations were used to identify cleavage sites in the peptide substrates of proteolytic enzymes. The procedure was similar for all the analyzed systems and included the following steps. First, a large number of docking simulations of a peptide substrate to an enzyme protein were performed without any prior knowledge of the localization of the enzyme’s active site. The peptide ligand was fully flexible and could freely sample the entire surface of the enzyme protein. During a single docking simulation, 10,000 structures of enzyme-substrate complexes were generated. Then, 1,000 models with the lowest protein-peptide interaction energy were selected. Finally, the top 10 models were identified using hierarchical clustering as the central structures of the ten largest clusters. The selected models were reconstructed to all-atom representation using the Modeller program [60]. A large number of docking simulations were completed for each analyzed system. The position of the cleavage site on each peptide was identified based on the analysis of a large number of resulting enzyme-substrate complexes generated using the above procedure (the main stages of the procedure used for the identification of cleavage sites in enzyme substrates are shown in Fig 2). For this task, we assumed that the peptide could be cleaved only in the substrate region that directly interacted with the enzyme’s catalytic site. More specifically, the cleavage site was assigned based on the histogram analysis of contacts between the carbonyl oxygen atoms of the peptide bonds in the substrate molecules and the center of the enzyme active site (the peptide bond for which the number of recorded contacts was the highest was considered as a potential cleavage site). Therefore, the center of the active site was approximated by the center of the line segment connecting the Cα atoms of the two catalytic Asp residues of the enzyme molecule. We considered a carbonyl oxygen atom to be in direct contact with the enzyme active site when the measured distance d_n between the oxygen atom of the substrate peptide and the center of the enzyme catalytic site was less than or equal to d_n ≤ 7 Å. The graphic illustrating the scheme used for distance d_n measurements is shown in Fig 3.

Fifteen different peptide substrates were analyzed. The first five were docked to a pepsin molecule (PDB ID: 4PEP [61]) and the procedure included 100 independent docking simulations, resulting in 1,000 all-atom models for each enzyme-substrate complex. The same procedure was applied for three substrates of renin protease (using crystal structure PDB ID: 3K1W [62]). The other seven peptides were docked to HIV-1 protease (structure extracted from PDB ID: 3EL1 [63]). In this case, 500 independent docking simulations were conducted resulting in 5,000 enzyme-substrate models for each system. A larger number of docking simulations were completed to account for the possible movement of the two flaps covering the HIV-1 protease active site. Therefore, the two protein fragments (residues 45 to 55 in chain A and chain B) were fully flexible during the docking simulation and no internal restraints were used. For a detailed description of the analyzed systems, see S2 Table.

Prediction of amyloid fibril structures

Structures of six different protofilaments with known experimental structures were predicted. The structures were assembled during docking simulations of five interacting peptide monomers. The final protofilament structures were then selected from a large set of protofilament models using structural and energetic criteria. To evaluate the accuracy of this procedure and the quality of the resulting models, we compared the obtained models with their experimental structures.

The standalone version of CABS-dock offers a wide range of adjustable parameters that allow modification of the docking procedure and system properties [20,25]. For instance, applying distance restraints can limit the conformational space of the simulated system, rejecting highly improbable states and allowing more efficient sampling of biologically relevant states. Due to the nature of the modeled system and the structural properties of the protofilaments created by in-register, parallel and identical peptide chains (which resulted in high translation symmetry of peptide chains forming amyloid fibrils) we applied weak distance restraints between the SC united atoms of corresponding amino acid residues in adjacent peptide chains (for a graphical representation of distance restraints, see Fig 1, panel b). The distance restraints did not determine the resulting shape of the predicted protofilament, but they defined the order of the monomer chains and favored in-register parallel stacking of peptides.

Projecting protein structures onto the lattice used in the CABS model results in a discrete representation of coordinates and allows very fast and efficient sampling of the conformational space. A large number of models can be generated in a relatively short time, which does not require large computing resources. Typically, this large set of models contains models that are very close to their native (or experimentally derived) structures and the most challenging step is to identify the best model. Our results show that among all the predicted models there were protofilament structures with low RMSD values, when compared to their experimental structures, in the range of 1.82 Å to 3.32 Å (see Table 2). Those structures were almost identical to the reference structures (see Fig 4).

Model scoring and selection of the most accurate structures in the proposed procedure were conducted based on two parameters: the interaction energy values estimated for the bonded monomers and the values of the pcaRMSD parameter, related to the translational symmetry of bonded peptide monomers in the protofilament structure. Scoring based solely on interaction energy values failed to identify high-quality models (see Table 2). Better results were obtained by scoring based on the pcaRMSD value, which was consistent with the fact that the structure of the protofilament, with its in-register parallel peptide chain architecture, should show high translational symmetry of interacting peptide monomers. The analysis of the interaction energy and pcaRMSD values calculated for all generated models showed a correlation with model accuracy, as evidenced by the low RMSD values when compared to the reference structure (refer to Figs 5 and S1). On average, the 1,000 models with the lowest energy/pcaRMSD values had an RMSD approximately 3.9 Å lower than the rest of the set, see S3 Table. However, it was impossible to identify the single most accurate structure among the entire set of models. To improve the scoring results, we tested a combination of these two parameters in two variants of the scoring procedure. In the first variant, the 1,000 models with the lowest interaction energy values were selected, and then the best model was chosen based on the lowest pcaRMSD value. The second variant of the scoring was similar. First, the 1,000 models with the lowest pcaRMSD values were selected, and then the model with the lowest interaction energy was chosen as the best-predicted structure. The first variant of the procedure proved to be more accurate and identified models with lower RMSD values (based on a comparison to the reference structures). The data for the best-identified structures are shown in Table 2, and the predicted structures for six protofilaments are shown in Fig 4.

The accuracy of the proposed protocol for predicting protofilament structures can be further increased by applying an MD-based refinement and model scoring procedure. In our recent work, we showed that even short MD simulations of CABS-generated models with reconstructed atomistic details can be efficiently used for model validation [26]. The receptor-peptide interaction energy values estimated during all-atom MD simulations with explicit solvent can be successfully used to identify the most accurate models, which are close to experimentally derived structures [26]. A similar procedure, followed by the analysis of the translational symmetry of monomers in the obtained protofilament models after MD optimization (e.g., using pcaRMSD parameter values), could probably increase prediction accuracy. Other scoring functions or geometry optimization protocols or consensus scoring functions based on the combination of different methods (e.g., Rosetta [64], GalaxyPepDock [65], HADDOCK [66], pepATTRACT [67]) could also improve the quality of the selected models.

Fast and accurate protocols for modeling protofilament structures can be the first step toward designing a method for predicting the structures of whole amyloid fibrils. It has been demonstrated that long fibril structures (formed by amyloidogenic insulin fragments) can be assembled using optimized protofilament structures as building blocks repeated along the long fibril axis [19]. We are currently working on a similar procedure to predict fibrils for other aggregating peptides, and the results will be published soon. The CABS-dock method enables the modification of distance restraints between the desired united atoms (or protein fragments) in a simulated system, expanding the possibilities for modeling different protofilament architectures with different symmetry patterns [37].

Identifying the positions of enzymatic cleavage sites in peptide substrates

The positions of enzymatic cleavage sites were identified using CABS-dock docking simulations for fifteen peptides: five pepsin substrates consisting of 8 to 21 amino acid residues, three renin substrates consisting of 8 to 18 amino acid residues and seven HIV-1 protease substrates containing of 8 to 16 amino acid residues. The predicted positions of the cleaved peptide bonds were then compared to experimental data from the MEROPS database [68] and published papers [69–74].

The three selected proteolytic enzymes were aspartate proteases, meaning they had two catalytic aspartate residues in their active sites. The exact location of the cleavage site on each peptide substrate was predicted by analyzing substrate-enzyme contact information from a large number of generated molecular models of enzyme-substrate complexes. Specifically, histogram analysis of all distances measured from the carbonyl oxygen atoms of all the substrate peptide bonds to the center of the enzyme’s active site was performed. According to the enzymatic catalysis mechanism of aspartate proteases, proteolytic cleavage begins with a nucleophilic attack of the catalytic water molecule on the carbonyl carbon atom of the scissile bond [38,39,75]. In this study, we assumed that cleavage occurred at the peptide bond whose carbonyl oxygen atom was most often located near the center of the protease’s active site in the set of predicted complex structures. The results of histogram analysis for the four selected systems are shown as bar plots in Fig 6.

The catalytic site of the pepsin molecule lies in the middle of a long cleft located between the two domains of the enzyme molecule (residues 1 to 175 and residues 176 to 327, respectively) [76]. Since the two domains of pepsin are not symmetrical, the local orientation of the substrate peptide chain at the center of the active site during enzymatic degradation must be strictly defined. This fact must be taken into account when predicting enzyme-substrate complexes. In the histogram analysis of pepsin-substrate systems, we included only those models in which the orientation of the substrate peptide chain along the cleft near the active site was similar to that observed in pepsin inhibitor structures (e.g., PDB ID: 1PSA [77]) and was also consistent with our previous studies [21]. In the case of renin-substrate complexes, we followed a similar approach as it has a structure comparable to that of pepsin, including the shape of the catalytic site (PDB ID: 3K1W [78]).

The positions of the cleavage sites identified using histogram analysis for two pepsin substrates KPAEF—IRL and HPHLSF—MAI show excellent agreement with experimental data [69]. For the first system, in 115 predicted models, the carbonyl oxygen atom in the 5^th peptide bond was located within the cutoff distance (d_n ≤ 7 Å) of the active site (see the left side of the panel (a) of Fig 6). The remaining oxygen atoms of other peptide bonds did not tend to stay near the binding site. In the second examined pepsin substrate, the cleavage site was also clearly identified on the 6^th peptide bond. The number of structures in which the oxygen atom was within the cutoff distance was 110. For the 5^th oxygen atom, the number of structures fulfilling the distance criteria was also high (67 models) but significantly smaller than for the 6^th oxygen (see the left side of the panel (b), Fig 6). The peptide chains of both ligands preferred to adopt an extended conformation along the cleft surrounding the active site. Two structures showing a ligand-bound conformation that allows cleavage of the accurately identified peptide bond are shown on the right side of panels (a) and (b), Fig 6.

For systems, which included HIV-1 protease, we conducted the histogram analysis by considering two different substrate peptide binding modes that were characterized by the opposite direction of the peptide chain extending along the enzyme active site (gray and black bars, Fig 6). This was necessary because HIV-1 protease is composed of two identical domains which make the enzyme molecule symmetrical. Therefore the same substrate can bind in the opposite direction of the peptide chain along the cleft between the two enzyme domains, which in turn shows two equivalent binding modes. The cleavage site was identified on the 5^th peptide bond of the two analyzed HIV-1 protease substrates GAETF—YVDGA, IRKIL—FLDGI which was positioned in the center of the two peptide chains (see Fig 6, panels c) and d)). To assess whether the length of the docked peptide could affect the prediction results of the cleavage site position (which was exactly in the center of the substrate peptide chain), we performed two additional docking simulations using longer variants of the HIV-1 substrate (seq.: GAETF—YVDGA). Since this peptide is part of GAG polyprotein [71], the main structural protein of HIV-1 and all other retroviruses [79], we constructed two longer substrate variants by adding additional amino acid residues to the C- or N-terminus of the substrate peptide chain according to the GAG sequence. The first substrate variant (seq.: TEPISGAETF—YVDGA) was created by adding five residues at the beginning of the peptide chain, while the second substrate variant (seq: GAETF—YVDGAANRET) was extended by adding five residues to the C-terminus. The cleavage sites in the two longer substrate variants were correctly identified using the proposed docking procedure, suggesting that length did not influence prediction accuracy. Bar plots presenting histogram analysis are shown in S2 Fig.

The active site of HIV-1 protease is covered by two identical flaps, residues from 45 to 55 of each domain [63]. The flap region is critical for binding substrates or inhibitors and also for catalytic activity of the enzyme [80]. Flap movement has been also reported in numerous works [81–83]. During all docking simulations, the flaps were fully flexible and transitions between open and closed conformations were observed. Interestingly, analysis of the MC trajectory (Cα trace representation only) generated during the docking simulation shows a correlation between the closed state of the two flaps and the ligand binding mode that is optimal for substrate cleavage (Fig 7). For bound ligands (with a distance of less than 9 Å between the center of the enzyme catalytic site and the scissile bond in substrate), we observed closed flap conformations in over 91% of cases (with a distance of less than 12 Å between Cα atoms of GLY:52 in chain A and GLY:52 in chain B of enzyme). In contrast, for unbound ligands, closed flap conformations were observed only 25% of the time. This was in agreement with previous findings showing that the flaps need to first open for the peptide to bind and that the protease interaction with the bound substrate influenced the flap opening frequency and interval [84].

The positions of cleavage sites identified for all fifteen enzyme-substrate systems are listed in Table 3. For comparison, we have included the prediction results obtained for pepsin and HIV-1 protease using PeptideCutter and HIVcleave servers, respectively. Our method accurately predicted cleavage sites for 12 systems. However, it should be noted that the predictions generated using the servers may overestimate the number of possible cleavage positions, depending on the chosen cutoff threshold value [85]. The servers identified between one and five cleavage sites. The HIVcleave method failed to identify cleaved bonds correctly for two HIV-1 protease substrates (seq.: SYFNLNPFEVL—QIDPE and seq.: NVVNSGGMVM—MVPGAG). To the best of our knowledge, no methods are currently available for the straightforward prediction of cleavage sites on renin peptide substrates.

Currently, computational methods used for predicting enzymatic cleavage sites in peptide or protein substrates can be divided into two groups based on their methodologies: sequence scoring methods and machine learning techniques [44]. Sequence scoring methods are generally less accurate, but they can produce results promptly. These methods use scoring functions based on large datasets of experimentally verified cleavage site positions in a large number of analyzed sequences. Predictions are made by comparing the query sequence to these known datasets. On the other hand, machine learning methods are more sophisticated and generally perform better, but they require well-assembled datasets for training, selection of an appropriate training model, and model evaluation and optimization [44]. In addition, machine learning methods can incorporate various heterogeneous features of cleaved substrates, including evolutionary information, physiochemical properties, and structural features during model training, which leads to better prediction accuracy. The proposed procedure for identifying cleavage sites based on docking simulations differs from the typical methods mentioned above. It relies entirely on enzyme-peptide interactions, meaning that the structural features of enzyme molecules and the dynamic character of the peptide ligand are essential for correctly positioning the substrate and allowing proteolytic degradation of a specific peptide bond. This approach is slower and requires more computational resources, but it takes into account all structural and dynamic features of the bound substrate, such as secondary structure preferences, multiple peptide chains, intramolecular disulfide bonds, exposed fragments of the peptide chain, and local substrate unfolding. In our recent work, we used a similar method to identify pepsin cleavage sites on the insulin molecule, and the simulation results were consistent with experimental findings [21]. We were able to correctly identify several cleavage sites on the substrate molecule. It is also worth noting that the molecular structure of insulin includes two peptide chains stabilized by three internal disulfide bonds, which determine its topology and conformation [21]. We are currently developing this method to enable the scanning of longer peptide/protein chains for potential cleavage sites and for the design of potent enzyme inhibitors.

In summary, we propose two new less obvious CABS-dock-based docking protocols. Since the CABS-dock algorithms are very flexible and allow easy-to-apply modifications of their default options, the range of integrative modeling schemes is very broad. The first of the two proposed modeling protocols addresses the early stages of amyloid formation: the assembly of small protofilaments. The test predictions were performed on several oligopeptide complexes whose three-dimensional structures have been determined experimentally. Five copies of peptide forming protofilaments were subject to simultaneous docking simulations, where one of the five peptides was treated as a flexible receptor, while the remaining four peptides were treated as fully flexible ligands. Known experimentally determined common structural regularities of expected protofilaments were imposed in a form of weak distance restraints, favoring parallel, in-register contacts of the peptide chains. The resulting structures after appropriate ranking were of good accuracy, providing structural data for further computational studies of amyloid formation. The second docking protocol was dedicated to predicting the positions of cleavage sites on peptide substrates of proteolytic enzymes. By analyzing many docking poses and specific geometries of the receptor-substrate complexes, those apparently associated with proteolytic degradation were identified in the analyzed systems. It has been shown that this strategy allows accurate identification of the cleaved peptide bonds, which is of great importance for structural biology and molecular medicine. Both protocols presented in this work prove that using a combination of multiscale modeling and known experimental facts, along with properly targeted search regions in efficient docking procedures, provides numerous opportunities for new computational studies of biomolecular complexes.