Towards Ab Initio Protein Folding: Inherent Secondary Structure

Jan 15, 2019 - Towards Ab Initio Protein Folding: Inherent Secondary Structure Propensity of Short Peptides from the Bioinformatics and Quantum-Chemic...
0 downloads 0 Views 2MB Size
Subscriber access provided by UNIV OF LOUISIANA

B: Biophysics; Physical Chemistry of Biological Systems and Biomolecules

Towards Ab Initio Protein Folding: Inherent Secondary Structure Propensity of Short Peptides from the Bioinformatics and Quantum-Chemical Perspective Martin Culka, Jakub Galgonek, Ji#í Vym#tal, Jiri Vondrasek, and Lubomír Rulíšek J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.8b09245 • Publication Date (Web): 15 Jan 2019 Downloaded from http://pubs.acs.org on January 20, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Towards Ab Initio Protein Folding: Inherent Secondary Structure Propensity of Short Peptides from the Bioinformatics and Quantum-Chemical Perspective Martin Culka, Jakub Galgonek, Jiří Vymětal, Jiří Vondrášek* and Lubomír Rulíšek* 1

Institute of Organic Chemistry and Biochemistry of the Czech Academy of Sciences,

Flemingovo náměstí 2, 166 10, Praha 6, Czech Republic * [email protected]; [email protected]

Abstract By combining bioinformatics with quantum-chemical calculations we attempt to address quantitatively some of the physical principles underlying protein folding. The former allowed us to identify tripeptide sequences in existing protein three-dimensional structures with a strong preference for either helical or extended structure. The selected representatives of pro-helical and pro-extended sequences were converted into ‘isolated’ tripeptides – capped at N- and C-termini and these were subjected to an extensive conformational sampling and geometry optimization (typically thousands to tens of thousands conformers for each tripeptide). For each conformer, QM(DFT-D3)/COSMO-RS free energy value was then calculated, Gconf(solv). The Gconf(solv) is expected to provide an objective, unbiased and quantitatively accurate measure of the 1 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

conformational preference of the particular tripeptide sequence. It has been shown that irrespective of the helical vs. extended preferences of the selected tripeptide sequences in context of the protein, most of the low-energy conformers of isolated tripeptides prefer the R-helical structure. Nevertheless, pro-helical tripeptides show slightly stronger helix preference than their proextended counterparts. Furthermore, when the sampling is repeated in the presence of a partner tripeptide to mimic the situation in a β-sheet, pro-extended tripeptides (exemplified by the VIV) show bigger free energy benefit than pro-helical tripeptides (exemplified by the EAM). This effect is even more pronounced in a hydrophobic solvent, which mimics less polar parts of a protein. This is in line with our bioinformatic results showing that majority of pro-extended tripeptides are hydrophobic. The preference for a specific secondary structure by the studied tripeptides is thus governed by the plasticity to adopt to its environment. In addition, we show that most of the ‘naturally occurring’ conformations of tripeptide sequences, i.e. those found in existing threedimensional protein structures, are within ~10 kcal.mol-1 from their global minima. In summary, our ‘ab initio’ data suggest that complex protein structure may start to emerge already at the level of their small oligopeptidic units which is in line with a hierarchical nature of protein folding.

2 ACS Paragon Plus Environment

Page 2 of 38

Page 3 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

1. Introduction Full understanding of protein folding – a key mechanism in translation of genetic information into the functional machinery - is one of the major scientific challenges for this millennium.1 The protein folding process in the cell can be steered by helper complexes called chaperones2 or even by the ribosome3,4 upon the gradual translation of the mRNA into a protein chain. However, more than 50 years ago it has been shown that an isolated protein can be unfolded and refolded also in solution,5 without any helper complexes. Therefore, the three-dimensional structural information is somehow encoded in the primary sequence of amino acids, which is in turn derived from the genetic information. Since the average protein is composed of several hundreds of amino acid residues6, it is difficult to envisage the process of unassisted protein folding as a completely random search for the global energy minimum. The scientific question yet to be solved is how the protein avoids this combinatorial explosion (so-called Levinthal’s paradox).7 It is clear that not only the final 3D structure, but also the folding mechanism has to be somehow encoded in the primary amino acid sequence. Furthermore, one protein sequence can fold into multiple stable 3D structures. This alternative folding, also termed misfolding, is, for example, the mechanism behind neurodegenerative diseases.8 Thus, it is not surprising that the protein folding has been subject of intense experimental and theoretical research over the past several decades.9 Given the fact that the folding process is relatively fast (usually < 1 s),10 it is not easy to study kinetic details of the process experimentally. Nevertheless, studies based on a single-molecule fluorescence resonance energy transfer (smFRET)11–14 and hydrogen exchange (HX) methods15,16 with nuclear magnetic resonance (NMR) or mass spectrometry (MS) analysis were able to provide many useful insights into protein folding mechanisms.

3 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

On the computational side, molecular dynamics (MD) is a straightforward choice to study a folding process in atomistic resolution. MD simulates temporal evolution of the protein structure within a physical energy model and if the simulation time is long enough, the folding process should be observable. Although the state-of-art supercomputers are able to run MD of the folding process of small stable proteins,17 the required millisecond timescales are still inaccessible for the majority of regular-sized proteins. Other computational approaches to protein folding either use sampling enhancement within physical energy model (e.g. biased MD or Monte Carlo techniques),18 or go further and replace the pure physical model with a set of constrains based on known experimental structures (homology modelling19) or their fragments (e.g. Rosetta20). Although those techniques have success in reconstructing previously unknown protein structures (e.g. CASP experiments21,22), they provide little clue about the folding process itself. Based on these experimental and computational studies, several theories of the protein folding mechanisms were formulated. All of them agree on the fact that to fight the Levinthal’s paradox, the native (folded) protein structure has to be assembled in a gradual fashion. The proposed theories differ, though, in the order of events, i.e. if the secondary structure elements are formed first and then aggregate into tertiary structure, or if the packed globular structure formation precedes the definite formation of secondary structures.23 Recent studies speak about foldons – independent protein fragments with defined supersecondary structure that are repeatedly folding and unfolding and eventually aggregate into the final tertiary structure.10 Specific foldons have been identified in cytochrome c,24,25 ribonuclease H,26 or ribosomal protein S6.27 In the light of the gradual nature of protein folding, it is interesting to look at the very initial steps of the process that might be represented by the conformations of the protein building blocks – short peptide fragments. In this study, we present the results of the protein structural database survey to 4 ACS Paragon Plus Environment

Page 4 of 38

Page 5 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

identify which of the 8000 (20x20x20) possible tripeptide fragments composed of naturally occurring amino acids have the highest tendency to prefer particular secondary structure type (helical or extended) in protein three-dimensional structures. The bioinformatic data are complemented by a large-scale quantum chemical study, employing calibrated methods presumably approaching the ~1-2 kcal.mol-1 accuracy in the computed energetics.28 Our primary aim is to correlate the occurrence of secondary structure elements in three-dimensional protein structures with the conformational spectrum of the corresponding tripeptide in solution to answer the question to what extent is the structural propensity encoded in the very small protein elements. In contrast to the vast majority of computational protein folding studies we used quantum mechanical (QM) free energies in model solvents described by presumably accurate COSMO-RS solvation theory. Although the QM methods have been used for evaluation of small peptides in past,29–31 this study provides a unique connection between database analysis and accurate QM freeenergy evaluation and may represent the initial step in understanding the protein structure from the first principles.

2. Methods Bioinformatic analysis of the protein structural database. In order to avoid redundancies and low-resolution structures, we used the Top8000 selection from PDB database as a large and representative sample of protein structures.32 Each of all 8000 possible variations of tripeptides composed of 20 canonical amino acids was assigned secondary structure employing the STRIDE algorithm.33 The same algorithm (STRIDE) was also used to identify discontinuities in the protein chain (i.e., omitted residues in experimental X-ray structures) and such tripeptides were ignored. Similarly, histidine anchors at both termini of the proteins, defined as NAA ≤ 8 amino acid sequences 5 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

followed by at least three consecutive His residues, were not included in our analysis. Alternatively, the secondary structure was assigned based on main chain dihedral angles (Figure 1) as described previously.34

Figure 1: Symbolic representation of a Ramachandran plot of main chain dihedral angles (φ and ψ). The plot is divided into three sectors, which represent R-helical, extended and other structure.

The propensity of a tripeptide X1X2X3 for helical or extended structure was calculated as a ratio between a number of tripeptides classified as HHH or EEE (all three amino acids in helical or all three in extended structure) and the total number of occurrences of X1X2X3 in the Top8000 database:

ℎ(𝑋1𝑋2𝑋3) =

|{𝑠 ∈ 𝑆(𝐴):𝑠𝑒𝑞(𝑠) = 𝑋1𝑋2𝑋3 ∧ 𝑠𝑡𝑟𝑡(𝑠) = 𝐻𝐻𝐻}| |{𝑠 ∈ 𝑆:𝑠𝑒𝑞(𝑠) = 𝑋1𝑋2𝑋3}|

(1) 6

ACS Paragon Plus Environment

Page 6 of 38

Page 7 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

𝑒(𝑋1𝑋2𝑋3) =

|{𝑠 ∈ 𝑆(𝐴):𝑠𝑒𝑞(𝑠) = 𝑋1𝑋2𝑋3 ∧ 𝑠𝑡𝑟𝑡(𝑠) = 𝐸𝐸𝐸}| |{𝑠 ∈ 𝑆:𝑠𝑒𝑞(𝑠) = 𝑋1𝑋2𝑋3}|

(2)

Initial Conformer Sets. Two conformer sets for each of the studied tripeptides were considered for the QM free energy evaluation (vide infra). First (reference) set was extracted from the Top8000 database taking all continuous occurrences of the tripeptides with the resolved structure of the side chains. All computationally studied tripeptides were capped at both termini: N-acetylated at the Nterminus and N-methylated at the C-terminus ( Figure 2).

Figure 2: Model of a capped tripeptide used for atomistic modeling throughout this study Second set was acquired using well-tempered35 bias-exchange36 metadynamics simulation of a capped tripeptide model Figure 2) in periodic box of explicit water using GROMACS37 and PLUMED38. The capped peptides in extended conformation were placed into a cubic simulation box with 35 Å edges and solvated by approx. 1400 TIP3P water molecules. Systems containing charged peptides were neutralized by addition of equivalent number of counterions (sodium cations and chlorine anions). All simulations utilized Amber ff99SB-ILDN39 force field (as provided in Gromacs force field library). Initially, the energy of the system was minimized by several hundreds of steps of steepest descent. Energy minimization was followed by 100 ps of solvent equilibration with restrained 7 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

peptide geometries (harmonic position restraints at non-hydrogen atoms with force constant 24 kcal.mol-1.Å-2). All molecular dynamics simulations were performed using 2 fs time steps, v-rescale thermostat (reference temperature 298.15 K, coupling constant 1.0 ps-1) and weak coupling isotropic (Berendsen) barostat (reference pressure of 1 bar and coupling constant 1.0 ps-1). All bond lengths were constrained by LINCS40 algorithm. Electrostatic interactions were evaluated by Particle-Mesh-Ewald (PME)41 algorithm. Lennard-Jones interactions were truncated at 14 Å. For the production run, we employed well-tempered bias-exchange metadynamics protocol for efficient sampling of free energy landscape. 8-11 replicas were utilized for systematic sampling of all φ, ψ and χ1 torsion angles in a peptide. One replica (neutral replica) were always kept without any bias, the others were biased by one-dimensional potential acting on a single torsion. The bias potential was updated according to the well-tempered algorithm (biasfactor 10) each 2 ps by a Gaussian with initial height of 0.06 kcal.mol-1 and a width of 0.2 rad. Exchanges between randomly selected replicas were attempted every 2 ps. All replicas were run for 100 ns of simulation time. We performed 2 independent metadynamics runs for each studied peptide. The peptide coordinates were saved every 20 ps. The structures from neutral replicas of both independent runs were collected (100,000 structures for each peptide in total) and taken for the quantum-chemical free energy evaluation (vide infra). Conformer sets of the X1X2X3 tripeptides with their β-sheet partners (also modeled as tripeptides) were obtained using 10 ns MD simulation of the two tripeptides in which the starting position was set to an antiparallel β-sheet orientation. The positions of main-chain heavy atoms of the β-sheet partner were kept frozen to emulate a rigid protein scaffold whereas no constraints were applied for the probe tripeptide (X1X2X3). We used NAMD42 with AMBER ff99SB43 force field for this sampling. The simulation has been carried out in the periodic box of explicit water molecules 8 ACS Paragon Plus Environment

Page 8 of 38

Page 9 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

employing the standard protocol: (1) energy minimization comprising 1000 conjugate-gradient steps; (2) heating from 0 to 298 K with velocity rescaling and harmonic potential restraints on all heavy atoms in the protein chain (k = 25 kcal.mol-1.Å-2); (3) 5ps equilibration; (4) 10 ns of production run. Conformers (snapshots) were saved every 25 ps. To acquire a more independent conformer sets, we performed above sampling procedure for both tripeptide pairs tested and then used PyMol44 to mutate the conformers into the other tripeptide pair (sometimes, this procedure is denoted as ‘morphing’). In total, 800 conformers of each β-sheet-like pair were taken for the quantum-chemical free energy evaluation (vide infra). The same MM setup was also used to obtain MM energies of the conformers in the set. In this case, each conformer was minimized using conjugate gradient method in Generalized-Born solvent model45,46 as implemented in NAMD. All conformer sets were filtered for redundancy using a Biopython47,48 script based on main chain and side chain dihedral angles using 20° tolerance. QM(DFT-D3)/COSMO-RS free energy calculations. The conformer’s free energies were evaluated employing previously calibrated DFT-D3/COSMO-RS methodology.49 This composite protocol may presumably yield the GDFT-D3/COSMO-RS values of individual conformers with ~1-2 kcal.mol-1 accuracy.28 All quantum-chemical calculations were performed using TURBOMOLE 7.2.50 Each conformer was subjected to geometry optimization in COSMO implicit solvent model51 with εr = 80, corresponding to water as a solvent. We employed BP8652,53 DFT functional with a moderate DZVP-DFT basis set54, which was recently shown to have the best price:performance ratio.55 The energies were corrected using D3BJ dispersion correction.55–57 The resulting optimized conformers were filtered for redundancy based on main chain and side chain dihedral angles using a Biopython47,48 script with a 20° tolerance. The non-redundant conformers were subjected to a 9 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

separate geometry optimization in the gas phase and subsequent frequency calculation to obtain the gas-phase thermodynamic quantities (Eq. 3 below). We then used the thermo program58 to calculate thermal contribution to the free energy using free-rotor rigid-rotor harmonic-oscillator (FRRRHO, sometimes denoted as quasi-RRHO) model with smooth transition between rigid and free rotors at 100 cm-1. Furthermore, we calculated two single-point energies on each nonredundant conformer using BP86 functional with def2-TZVPD59 basis set and D3BJ dispersion correction - one in gas phase and the second one using COSMO with εr = ∞. We then used the COSMOtherm17 program,60 employing the “BP_TZVPD_FINE_C30_1701.ctd” parametrization file and FINE cavities ($cosmo_isorad keyword), to calculate the accurate COSMO-RS61,62 solvation energy in three representative solvents of increasing hydrophobicity - water, 1-octanol and hexane. The resulting free energy for each conformer was calculated as a sum of the BP86D3/def2-TZVPD gas-phase energy, COSMO-RS solvation energy and FRRRHO thermal contribution (Eq. (3)).

𝐺𝑐𝑜𝑛𝑓(𝑠𝑜𝑙𝑣) = 𝐸𝑒𝑙(𝐷𝐹𝑇 ― 𝐷3𝐵𝐽) + 𝐺𝑠𝑜𝑙𝑣(𝐶𝑂𝑆𝑀𝑂 ― 𝑅𝑆) + 𝐸𝑍𝑃𝑉𝐸 ―𝑅𝑇𝑙𝑛(𝑞𝑡𝑟𝑎𝑛𝑠𝑞𝑟𝑜𝑡𝑞𝑣𝑖𝑏) +𝑝𝑉 (3) For calculation of β-sheet interaction energies, we performed geometry optimization for the complex and all single point and frequency calculations for both the complex and isolated partners. The interaction energy Gint was calculated by subtracting the free energies of individual partners from the complex energy. Statistical analysis. We used previously defined34 regions of the Ramachandran plot (Figure 1) to assign R-helical, extended or other structure character for each conformer of each tripeptide considering the main-chain dihedral angles of the central amino acid (in contrast to the situation in 10 ACS Paragon Plus Environment

Page 10 of 38

Page 11 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

the protein, the conformation of the peripheral amino acids cannot be considered as well defined). The propensity of each conformer set to each secondary structure type was calculated as a Boltzmann-weighted ratio using the QM-derived free energy and temperature of 298 K according to the formula:

∑𝑠𝑠𝑒 ―𝐺 𝑅𝑇

𝑃𝑠𝑠 = ∑

𝑒

(4)

―𝐺 𝑅𝑇

𝑎𝑙𝑙

where Pss is propensity for particular secondary structure (R-helix, extended or other), G is the free energy of the conformer, T is the absolute temperature (here we use 298 K) and R is the universal gas constant. Given the fact that our free energies suffer from an inherent error (vide supra), we added a random error from Gaussian distribution with σ = 1 kcal.mol-1 to each energy as a model of nonsystematic error and calculated the secondary structure propensity using the Eq. (4) thousand times. We then used the average secondary structure propensity and its standard deviation as a representative value and an error bar, respectively. Hydrogen bond analysis. In this study, we recognize a hydrogen bond based on hydrogen – acceptor distance lower than 2.5 Å and donor – hydrogen – acceptor angle greater than 130°. Within our tripeptide model, we consider hydrogen bonds between main chain CO and NH groups and between polar side chains and respective counterparts on the main chain. The script for the hydrogen bond analysis was written using Biopython.47,48

11 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 38

3. Results 3.1. Bioinformatics: Identification of tripeptides with strong secondary structure preference in Top8000 database. Using the Top8000 subset of PDB database and the STRIDE algorithm, each of 8000 tripeptides was statistically analyzed with respect to the preference for the secondary structure. The full set of data is available in the SI whereas the tripeptides with the strongest preference for the either the right-handed (R) helices or extended structures are listed below in

Table 1 and Table 2. For comparison, we also calculated secondary structure propensity based purely on regions of the Ramachandran plot34 (c.f. Tables 1 and 2). It can be seen that in general, the pro-R-helical tripeptides contain charged residues and alanines, while the pro-extended are rather hydrophobic or have aromatic side chains. For further analysis (vide infra), we have chosen two tripeptides with the highest propensity for each secondary structure type: EAM and KAM (-1 and +1 charge, respectively) from the R-helical set and hydrophobic VIV and VVV from the extended set. In addition, we have chosen untypical representative of each group for further analysis – hydrophobic pro-R-helical ALA tripeptide and polar pro-extended IYI tripeptide. The selected representative tripeptides are highlighted in bold in

Table 1 and Table 2.

Table 1. Triplets having count of samples >= 100 and high R-helix propensities triplet

count

STRIDE(h)

STRIDE(e)

Ramachandran (h)

Ramachandran (e)

EAM KAM

234 171

0.803 0.789

0.03 0.029

0.838 0.830

0.060 0.058

12 ACS Paragon Plus Environment

Page 13 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

EAL EEA AEA QEA EAA QAL MEQ EAI ALA AAR QAA EMA AMQ XLA AAL EEL MEA EEM

1470 895 1063 350 1086 787 102 782 1724 960 658 204 119 104 1657 1162 184 208

0.776 0.759 0.746 0.743 0.742 0.740 0.735 0.733 0.731 0.730 0.728 0.725 0.723 0.721 0.720 0.719 0.717 0.716

0.033 0.023 0.048 0.046 0.038 0.047 0.088 0.083 0.073 0.048 0.038 0.074 0.034 0.087 0.078 0.025 0.071 0.024

0.829 0.808 0.791 0.791 0.795 0.778 0.784 0.753 0.776 0.778 0.766 0.789 0.748 0.74 0.774 0.815 0.761 0.798

0.061 0.051 0.094 0.080 0.067 0.097 0.118 0.121 0.108 0.071 0.058 0.098 0.076 0.144 0.109 0.065 0.092 0.077

Table 2. Triplets having count of samples ≥ 100 and high extended conformation propensities triplet

count

STRIDE(h)

STRIDE(e)

Ramachandran (h)

Ramachandran (e)

VIV VVV VII IIV III VLV VVI IVV IVI VFV VIF VYV VYI FVV LVV IVC YVI IYI FVI ILV

516 840 358 315 245 859 579 495 353 329 232 276 198 335 759 103 206 151 220 525

0.081 0.093 0.12 0.089 0.102 0.135 0.09 0.093 0.102 0.103 0.116 0.083 0.141 0.122 0.154 0.155 0.194 0.086 0.159 0.185

0.802 0.794 0.779 0.778 0.763 0.753 0.746 0.737 0.737 0.729 0.724 0.721 0.712 0.704 0.701 0.699 0.689 0.689 0.686 0.686

0.081 0.096 0.12 0.092 0.106 0.14 0.09 0.095 0.105 0.106 0.142 0.087 0.146 0.128 0.154 0.175 0.194 0.099 0.164 0.192

0.868 0.848 0.835 0.851 0.829 0.767 0.864 0.848 0.83 0.781 0.75 0.804 0.763 0.776 0.775 0.728 0.714 0.755 0.75 0.73

13 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 38

3.2. Secondary structure propensity at the tripeptide level obtained by the DFTD3/COSMO-RS conformational sampling. To address the central question of this study; i.e. whether the propensity for secondary structure is somehow encoded in the smallest building blocks of proteins, such as tripeptides, a thorough conformational sampling of the six tripeptides discussed above (pro-R-helical EAM, KAM, ALA, pro-extended VIV, VVV, IYI) was carried out. The models were capped on both termini ( Figure 2). Two sets of conformers for each of the six tripeptides were considered: (i) ConfSetEXP; the reference (experimental) set comprised all conformers of the particular tripeptide extracted from the Top8000 database (i.e., the same set that was used for the bioinformatic analysis); (ii) ConfSetMETDYN; in silico-generated probe set comprising the conformers obtained by an extensive molecular dynamics (MD) sampling in explicit solvent using the bias-exchange metadynamics method. In this simulation, biasing potential was gradually added to discourage the MD simulation from sampling the part of the dihedral angle space that has already been visited. Furthermore, each simulation was run in several replicas, where different dihedral angles were biased and the conformers were occasionally exchanged among the replicas. Thus, full coverage of dihedral angle space of each tripeptide was ensured. The both sets of conformers – ConfSetEXP and ConfSetMETDYN - were subjected to QM(DFT-D3)/COSMO-RS post-processing, which resulted in identifying the unique local minima on presumably accurate potential (free) energy surface (PES), defined by Eq. (3). It can be seen from Table 3 that the ConfSetEXP set is up to two orders of magnitude smaller than the ConfSetMETDYN set.

Table 3. Number of unique conformers of every optimized tripeptide in both sets EAM

KAM

ALA

VIV

14 ACS Paragon Plus Environment

VVV

IYI

Page 15 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

ConfSetEXP (from Top8000) ConfSetMETDYN (metadynamics)

101

103

83

71

60

50

23258

31216

839

2522

1513

7617

The computed conformational free energies exhibit roughly Gaussian distribution ( Figure 3 and Figure 4). We divide the conformers into three groups based on secondary structure propensity; R-helical and extended, both are common in protein structures, and other structures (which are, in fact, not that much populated in proteins and can be found in unstructured parts of proteins or rather rare L-helices34). It can be seen in both conformer sets ( Figure 3 and Figure 4) that the energy histograms of the pro-R-helical tripeptides are shifted towards lower energies compared to their pro-extended and other counterparts (1–5 kcal.mol-1). The extended conformers are slightly more disfavored relative to R-helical conformers in the case of charged EAM and KAM, which are strongly pro-R-helical according to our database survey. This observation is the most pronounced in hexane. The pro-R-helical but hydrophobic ALA behaves similarly as the pro-extended VIV, VVV and IYI.

15 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 3. Free energy histograms for the two conformer sets in two solvents of the pro-R-helical tripeptides (EAM, KAM, ALA) divided according to the secondary structure types. The energies are reported as relative to the lowest conformer present in the data set. For each histogram, the solid colored line represents the position of the lowest bar (G0), while the dashed line represents the mean value of the Gaussian fit (μ). The histograms were truncated at 20 kcal.mol-1 relative to the lowest conformer. 16 ACS Paragon Plus Environment

Page 16 of 38

Page 17 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 4. Free energy histograms for the two conformer sets in two solvents of the pro-extended tripeptides (VIV, VVV, IYI) divided according the secondary structure types. The energies are reported as relative to the lowest conformer present in the data set. For each histogram, the solid colored line represents the position of the lowest bar (G0), while the dashed line represents the mean value of the Gaussian fit (μ). The histograms were truncated at 20 kcal.mol-1 relative to the lowest conformer.

17 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Next noticeable thing is that the conformers in the ConfSetEXP (derived from Top8000 database) lie in the lower energy part of the much larger set ConfSetMETDYN. This indicates that at the tripeptide level, the nature has chosen to “employ” a low energy conformer subset for building the proteins. However, the free energy window spanned by the conformers in the ConfSetEXP is rather broad, up to 16 kcal.mol-1 in water (Figure 3 and Figure 4). Admittedly, this is a bit higher that one would expect and it implies that the forces, which stabilize the particular tripeptide conformation in tertiary structure of the native protein, can overcome such values. We calculated the weighted secondary structure propensity of the two conformer sets using the free energies as Boltzmann weights (Eq. (4)) and performed thousand statistical trials to account for inherent error in computed free energies (see Methods). The resulting average propensities and standard deviations are depicted in Figure 5. It can be seen that in water, all studied tripeptides show the significant preference for the helical structure. This helical preference decreases with increasing hydrophobicity of the solvent. For the pro-extended tripeptides (VIV, VVV and IYI), it can be noticed that the conformers from the ConfSetEXP set reproduce the propensity for the extended structure in hexane quite well

18 ACS Paragon Plus Environment

Page 18 of 38

Page 19 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 5. Computed propensities (Pss, see Eq. (4)) for R-helical structure – Phelix - of the central amino acid of the tripeptide for the ConfSetEXP and ConfSetMETDYN sets in different solvents. Error bars represent standard deviations when Gaussian-distributed random error is added thousand times. The propensities for extended structure are complementary (and approximately equal to Pext = 1 - Phelix, since those for other structures are close to zero).

3.3. Role of intramolecular hydrogen bonds. Since internal hydrogen bonds dominate the helical structure, we analyzed how the presence of internal hydrogen bonds correlate with the energy of the conformer. We found out that the low energy conformers indeed contain one or more intramolecular hydrogen bonds. On the other hand, when we divide the energy histograms into conformers with at least one hydrogen bond and the 19 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ones without any hydrogen bond we can see that both energy histograms span over a wide energy range (Figure 6 and Figures S1-S4 in the SI). Concerning the ConfSetEXP set, we can clearly see that the pro-helical tripeptides (EAM, KAM, ALA) contain more conformers with hydrogen bonds than the pro-extended ones (VIV, VVV and IYI). The calculations also predict that, in general, the conformers without hydrogen bonds have higher energies. Interestingly, these are consistently shifted towards lower energies with increasing solvent hydrophobicity (decreasing polarity). The same is true for the ConfSetMETDYN conformers. Moreover, the relatively large number of conformers in the ConfSetMETDYN set enabled us to perform Gaussian curve fits and calculate the respective mean values, μ. We can see that for the pro-R-helical tripeptides the conformers containing hydrogen bonds are systematically shifted towards lower energies (1-2 kcal.mol-1), which increases with the solvent hydrophobicity, up to 8 kcal.mol-1. In the case of pro-extended conformers the shift is almost zero in water (c.f. EAM and VIV in Figure 6) and increases to only 2 kcal.mol-1 in hexane.

20 ACS Paragon Plus Environment

Page 20 of 38

Page 21 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 6. Energy histograms of the conformers of the pro-R-helical tripeptide EAM and proextended tripeptide VIV from the ConfSetMETDYN set divided into two groups according the presence/absence of intramolecular hydrogen bonds. Dashed lines represent fits with a Gaussian function with the respective mean values (). 3.4. β-sheet partners increase the relative stability of extended conformers. It can be expected that extended conformers can be stabilized by the interaction with their partners in the β-sheet structures; this term is missing in our simple tripeptide model. To evaluate this energetic contribution, we chose a particular VIV conformer from the ConfSetEXP set that has 21 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

two β-sheet partners in the parent protein structure (PDB ID 1DL5, chain B). Energetically, the conformer is 8.5 kcal.mol-1 higher than the lowest VIV conformer. To be consistent, we truncated the β-sheet partners to capped tripeptides (antiparallel FLF and parallel IFV). Each partner tripeptide has 3-4 hydrogen bonds to the central VIV tripeptide (Figure 7).

Figure 7. VIV tripeptide fragment and its β-sheet partners (antiparallel FLF and parallel IFV) in the 1DL5 PDB structure.

For comparison, we identified a rare case where the otherwise strongly pro-R-helical tripeptide EAM assumes extended conformation in the middle of an antiparallel β-sheet (PDB ID 1H2C). This extended EAM conformer is energetically 10.0 kcal.mol-1 above the global EAM minimum and has PIW tripeptide as its interaction partner. Interaction free energy of each pair of tripeptides (VIV-FLF, VIV-IFV and EAM-PIW) was calculated as difference between free energy of the complex and the sum of free energies of the two isolated tripeptide partners. As can be seen in

22 ACS Paragon Plus Environment

Page 22 of 38

Page 23 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Table 4, the interaction free energy is -2 to -6 kcal.mol-1 in water and it increases with the hydrophobicity of the solvent, as can be expected. In the gas phase, the interaction free energies (Gint ) of the two tripeptides in the pair reach -24 to -44 kcal.mol-1, which is in agreement with a previous study carried out for polyalanine.63

Table 4. Interaction free energies (kcal.mol-1) of a high-energetic VIV conformer with two β-sheet partner tripeptides Water 1-octanol

Hexane Gas-phase

VIV-FLF (antiparallel)

-5.1

-9.6

-15.9

-24.2

VIV-IFV (parallel)

-5.9

-9.1

-17.4

-26.5

EAM-PIW (antiparallel) -2.4

-9.4

-27.4

-44.1

To further corroborate the above findings, we carried out an independent (and unbiased) computational experiment by performing MD simulations in explicit water for the VIV-FLF and EAM-PIW pairs and their swapped counterparts (EAM-FLF and VIV-PIW) with the fixed β-sheet conformations of the partners (FLF, PIW). The conformers obtained by this procedure were subjected to the same DFT-D3/COSMO-RS post-processing resulting into conformer free energies (Gconf(solv)) and interaction free energies (Gint). It has been found that majority of resulting conformers show propensity for extended structure. For both VIV and EAM, the free energies of the tripeptide conformers obtained by the above sampling (with a fixed β-sheet partner) are energetically higher (> 3 kcal.mol-1) with respect to their counterparts in the ConfSetMETDYN set. (vide supra). Adding the computed Gint values to Gconf(solv) enabled us to estimate the relative stabilization of the tripeptide in the corresponding pair. The free energy histograms for the VIV23 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

FLF, EAM-FLF, EAM-PIW and VIV-PIW pairs are depicted in Figure 8. It can be seen that for both VIV and EAM, the energies are lowered towards (or even below) zero, which means that the interaction with FLF or PIW in a β-sheet-like structure makes this extended conformer preferred over the pro-helical conformer, the latter being favored in the isolated tripeptide. Clearly, the VIV conformers benefit from the interaction to a greater extent than the EAM conformers, since the mean of the interaction-corrected free energy histograms for VIV is negatively shifted with respect to EAM in all the cases (Figure 8), even though the differences are slightly smaller in the case of PIW (natural partner of EAM in the 1H2C structure). Interestingly, the histograms shift away from each other with increasing hydrophobicity of the solvent. This means that while increasing fraction of VIV conformers prefer to be in β-sheet when the hydrophobicity of the surroundings increases, for EAM the number of conformers that benefit from interaction with FLF or PIW decreases with the increasing hydrophobicity. When we analyze these structures visually, we typically see an interaction of the charged glutamate with the main chain amino hydrogens, which competes the βsheet-like interaction.

24 ACS Paragon Plus Environment

Page 24 of 38

Page 25 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 8. Histograms of free energies of VIV and EAM tripeptide conformers relative to the independent VIV and EAM conformer sets corrected by interaction energies with a β-sheet-like partner tripeptide FLF or PIW. The partner tripeptides were chosen based on β-sheets found in Xray structures – VIV-FLF from PDB ID IDL5 and EAM-PIW from PDB ID 1H2C. Dashed lines represent fits with a Gaussian function, where μ is the mean value. 25 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

3.5. Comparison of QM(DFT-D3)/COSMO-RS energies with force fields (MM). Admittedly, the QM methodology used here to obtain free energies is rather time-consuming, compared to the common force-field-based molecular mechanics (MM) approaches. However, by comparing the MM energies within the Generalized-Born solvation model with the DFTD3/COSMO-RS energies, we may immediately see that the results for the secondary structure propensity are dramatically different. In general, the propensities are closer to 0.5, or even sometimes, the extended structure is preferred over the helical, which is in contrast to the QMbased results. In addition, the errors derived from the thousand statistical trials are much bigger for the propensities, where MM energies were used as Boltzmann weights in Eq. (4) (Figure 9). On the other hand, the energies obtained by the lower-level QM method used for geometry optimization (BP86-D3/DZVP-DFT//COSMO) yield similar results as the more sophisticated free energies used throughout this study.

26 ACS Paragon Plus Environment

Page 26 of 38

Page 27 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 9. Propensity for R-helical structure for the selected tripeptides – comparison of the expected propensity from the Top8000 database (black) with energy-weighted propensity derived from metadynamics sampling followed by geometry optimization and energy calculation in water using different levels of theory.

4. Discussion In the present study, we identified tripeptide sequences that occur predominantly in only one type of secondary structure – R-helical or extended (α-helix or β-sheet in broader protein context). This preference can be governed by many factors – among others hydrophobicity, number of hydrogen bond donors and acceptors, or bulkiness of the sidechains. We noticed that the pro-R-helical tripeptides usually contain alanines and charged amino acids, while the pro-extended are mostly 27 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

hydrophobic. This reflects the fact that the extended structure is stabilized by long-range mainchain interactions in β-sheets, which might be further stabilized by formation of hydrophobic patches from the apolar and aromatic side chains. We next set to find out how much of this secondary structure preference can be reproduced on the level of capped tripeptides, i.e. to find out if the tripeptides themselves possess a strong driving force towards the preferred secondary structure type and can thus serve as folding initiation sites. We compared isolated tripeptide fragments extracted from the experimental structures (Top8000 database (denoted ConfSetEXP), with a large independent set of conformers derived from enhanced sampling molecular dynamics (bias-exchange metadynamics, denoted ConfSetMETDYN) and used calibrated quantum-chemical (QM) methodology to obtain accurate conformational free energies (Gconf(solv)). It has been shown that the selected representatives of pro-R-helical (EAM, KAM, ALA) and proextended (VIV, VVV, IYI) tripeptides always favor the R-helical conformation in implicit water solution. However, if hexane is used as the implicit solvent model instead, the expected secondary structure preference is restored in the case of the conformers of pro-extended (VIV, VVV, IYI) tripeptides derived from the Top8000 database. Thus, it seems that the hydrophobicity of the environment accounts for a substantial amount of the stabilization energy, which is in line with the fact that the majority of pro-extended tripeptides are hydrophobic. However, in a much larger and unbiased set of conformers generated by the metadynamics (ConfSetMETDYN), it seems that the proextended tripeptides (VIV, VVV, IYI) have little reason to prefer the extended conformation, since the formation of the inner main-chain hydrogen bonds is obviously not hindered enough by the bulky hydrophobic sidechains during the sampling procedure. Thus, there is no other driving force to favor extended conformers as it is in the ConfSetEXP set and the secondary structure propensity 28 ACS Paragon Plus Environment

Page 28 of 38

Page 29 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

of the pro-extended tripeptides is also R-helical in the independent ConfSetMETDYN set. Still, the charged pro-R-helical EAM and KAM tripeptides show slightly stronger preference for R-helix than the rest of the studied tripeptides, which are uncharged (1 kcal.mol-1 average free energy difference of the low-lying conformers). Although the data from the ConfSetMETDYN set indicate that the capped tripeptides prefer Rhelical structure in all the solvents considered, the standard deviations derived from the random error addition provide further interesting information. It can be seen that the strongly pro-R-helical tripeptides EAM and KAM show strong R-helix propensity with a small deviation in water, while the R-helix preference is less pronounced for the non-polar tripeptides. This may indicate that the neutral tripeptides with non-polar side chains (ALA, VIV, VVV and IYI) will be more susceptible to give up the R-helical structure whenever they have an opportunity to do so. This general helicity preference can be related to the fact that pro-R-helical structures form intramolecular hydrogen bonds, which are perfectly reproduced within our tripeptide model. Indeed, we observed that lowest conformers always contain at least one intramolecular hydrogen bond, in agreement with previously published QM studies29,64 and also with the gas-phase QM calculations coupled with gas-phase UV/IR experiments showing that capped tripeptides form helical structure.65 When various secondary structures are compared using high-level QM methods and implicit water solvent model, the helical structures come out as most stable for trialanine and tetraalanine.31 Admittedly, the situation might differ in explicit water solution, where the water molecules can form bridges between adjacent carbonyl and amide groups, thus stabilizing the extended conformations, albeit in non-capped (zwitterionic) tripeptides.66 In several cases, the conformers denoted as other are the second-best within the ConfSetMETDYN set. These conformers often resemble L-helix, which contains the inner hydrogen bonds as well. 29 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The presence of an intramolecular hydrogen bond seems to be enough to prefer the pro-L-helical conformers over the extended conformers in some cases which might be an inauspicious aspect of our rather small tripeptide model. Note that although the L-helical conformations are not very common in experimental structures, short L-helices were also identified in the protein structural database.67 In addition, few conformers denoted as other were found in our ConfSetEXP set for the pro-helical tripeptides EAM, KAM and ALA. We inspected the respective protein structures and found out that the corresponding tripeptides come exclusively from the loop or turn regions and they contain inner hydrogen bonds in some of the cases. These presumably rare conformers might be necessary to connect the secondary structure elements and we have indeed found them in ConfSetMETDYN set. It is clear from above that the difference between pro-R-helical and pro-extended tripeptides cannot lie only in the difference between the low-energetic conformers, but rather in the behavior of the whole ensembles of the conformers, also given the wide range (up to 16 kcal.mol-1) of energies represented by the conformers derived from the experimental protein structures. On one hand we observe that the pro-helical tripeptides (EAM, KAM, ALA) benefit more from internal hydrogen bonding, since bigger fraction of low-energy conformers of these tripeptides contain at least one hydrogen bond. On the other hand, the extended main chain conformations are usually stabilized by the long-range hydrogen bonds in the β-sheet structures in proteins, which are missing in the separated tripeptide model and these conformers are artificially penalized in the separate tripeptide model by at least 2 kcal.mol-1. On the example of a VIV and EAM conformers from βsheet we found out that interaction free energies within a β-sheet reach 2-5 kcal.mol-1 in water and increase with the solvent hydrophobicity. When we compare whole conformer ensembles for proextended VIV and pro-helical EAM obtained from sampling in a β-sheet-like pair of tripeptides, 30 ACS Paragon Plus Environment

Page 30 of 38

Page 31 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

we see that VIV benefits more from the interaction with the partner than EAM. This difference even increases with solvent hydrophobicity, which is in line with the observation on separated conformers from ConfSetEXP presented above. VIV thus turns out to be more susceptible to form a β-sheet when given an opportunity and this preference is further enhanced if this interaction happens in a hydrophobic patch. Indeed, it is well known that the environment polarity is crucial for the stabilization of the β-sheet68 and that side chain interactions also play an important role.69 Given the fact that majority of pro-extended tripeptides found in our database survey are hydrophobic, these observations from ab initio tripeptide treatment are in line with trends in experimental protein structures. Furthermore, our observations that pro-extended tripeptides might prefer helical structure when separated and only switch to extended conformers when a suitable partner is present is in line with outcomes of NMR and CD experiments, which have shown that separated β-sheet elements of the native structure transiently fold to helical structure and only upon tertiary contacts during the folding process form the final β-sheets.70,71

5. Conclusions In summary, we identified tripeptide fragments with strong R-helical or extended structure preference in the protein structural database (PDB). Using quantum-chemical (QM) calculations, we have shown that all tested tripeptide fragments when taken as isolated species, have preference to R-helical structure. Nevertheless, the pro-R-helical tripeptides benefit more from internal hydrogen bonds than the pro-extended ones. When given an extended β-sheet-like tripeptide partner, the extended structure preference is observed for VIV tripeptide, which indeed prefers the extended structure in proteins. In contrast, the pro-R-helical EAM benefits less from the interaction with its partner in the β-sheet. The reason why we find tripeptide fragments with strong secondary structure preference in the protein database (PDB) thus lies in their plasticity to interact with the 31 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

environment rather than in different conformational space of isolated tripeptides. Although the process of protein folding is complex, our results indicate that the preference to form different secondary structure types might be tracked back to small tripeptide fragments using ab initio conformer sets and accurate QM methods. It is tentative to speculate that the tripeptide fragments with a strong preference for one secondary structure type might serve as folding initiation sites and prime the formation of secondary structure elements.

Acknowledgements. The financial support of the Grant Agency of the Czech Republic is gratefully acknowledged (grant 17-24155S). This work was supported by The Ministry of Education, Youth and Sports from the Large Infrastructures for Research, Experimental Development and Innovations project “IT4Innovations National Supercomputing Center – LM2015070“. Access to computing and storage facilities owned by parties and projects contributing to the National Grid Infrastructure MetaCentrum provided under the programme "Projects of Large Research, Development, and Innovations Infrastructures" (CESNET LM2015042), is greatly appreciated.

Supporting Information. PDF file with figures related to the hydrogen bond analysis (Figures S1-S4). Excel spreadsheet with a table of secondary structure propensities for all 8000 possible (20x20x20) tripeptides identified in Top8000 protein structure database (an extension of Table 1 and Table 2). ZIP file with all structures of capped tripeptides together with DFT-D3/COSMO-RS free energies. References (1)

So Much More to Know. Science. 2005, 309, 78–102.

(2)

Mashaghi, A.; Kramer, G.; Lamb, D. C.; Mayer, M. P.; Tans, S. J. Chaperone Action at the SingleMolecule Level. Chem. Rev. 2014, 114, 660–676.

(3)

Nilsson, O. B.; Hedman, R.; Marino, J.; Wickles, S.; Bischoff, L.; Johansson, M.; Müller-Lucks, A.; Trovato, F.; Puglisi, J. D.; O’Brien, E. P.; et al. Cotranslational Protein Folding inside the Ribosome Exit Tunnel. Cell Rep. 2015, 12, 1533-1540. 32 ACS Paragon Plus Environment

Page 32 of 38

Page 33 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(4)

Javed, A.; Christodoulou, J.; Cabrita, L. D.; Orlova, E. V. The Ribosome and Its Role in Protein Folding: Looking through a Magnifying Glass. Acta Crystallogr. Sect. D Struct. Biol. 2017, 73, 509–521.

(5)

Anfinsen, C. B.; Haber, E.; Sela, M.; White, F. H. The Kinetics of Formation of Native Ribonuclease during Oxidation of the Reduced Polypeptide Chain. Proc. Natl. Acad. Sci. U. S. A. 1961, 47, 1309–1314.

(6)

Skovgaard, M.; Jensen, L. J.; Brunak, S.; Ussery, D.; Krogh, A. On the Total Number of Genes and Their Length Distribution in Complete Microbial Genomes. Trends Genet. 2001, 17, 425–428.

(7)

Levinthal, C. How to Fold Graciously. In Mössbaun Spectroscopy in Biological Systems Proceedings; 1969; pp 22–24.

(8)

Bolshette, N. B.; Thakur, K. K.; Bidkar, A. P.; Trandafir, C.; Kumar, P.; Gogoi, R. Protein Folding and Misfolding in the Neurodegenerative Disorders: A Review. Rev. Neurol. (Paris). 2014, 170, 151–161.

(9)

Abaskharon, R. M.; Gai, F. Meandering Down the Energy Landscape of Protein Folding: Are We There Yet? Biophys. J. 2016, 110, 1924–1932.

(10)

Englander, S. W.; Mayne, L. The Nature of Protein Folding Pathways. Proc. Natl. Acad. Sci. 2014, 111, 15873–15880.

(11)

Soranno, A.; Buchli, B.; Nettels, D.; Cheng, R. R.; Muller-Spath, S.; Pfeil, S. H.; Hoffmann, A.; Lipman, E. A.; Makarov, D. E.; Schuler, B. Quantifying Internal Friction in Unfolded and Intrinsically Disordered Proteins with Single-Molecule Spectroscopy. Proc. Natl. Acad. Sci. 2012, 109, 17800–17806.

(12)

Voelz, V. A.; Jäger, M.; Yao, S.; Chen, Y.; Zhu, L.; Waldauer, S. A.; Bowman, G. R.; Friedrichs, M.; Bakajin, O.; Lapidus, L. J.; et al. Slow Unfolded-State Structuring in Acyl-CoA Binding Protein Folding Revealed by Simulation and Experiment. J. Am. Chem. Soc. 2012, 134, 12565– 12577.

(13)

Pirchi, M.; Ziv, G.; Riven, I.; Cohen, S. S.; Zohar, N.; Barak, Y.; Haran, G. Single-Molecule Fluorescence Spectroscopy Maps the Folding Landscape of a Large Protein. Nat. Commun. 2011, 2, 493.

(14)

Liu, J.; Campos, L. A.; Cerminara, M.; Wang, X.; Ramanathan, R.; English, D. S.; Munoz, V. Exploring One-State Downhill Protein Folding in Single Molecules. Proc. Natl. Acad. Sci. 2012, 109, 179–184.

(15)

Krishna, M. Hydrogen Exchange Methods to Study Protein Folding. Methods 2004, 34, 51–64.

(16)

Krishna, M. M. G.; Lin, Y.; Mayne, L.; Walter Englander, S. Intimate View of a Kinetic Protein Folding Intermediate: Residue-Resolved Structure, Interactions, Stability, Folding and Unfolding Rates, Homogeneity. J. Mol. Biol. 2003, 334, 501–513.

(17)

Lindorff-Larsen, K.; Piana, S.; Dror, R. O.; Shaw, D. E. How Fast-Folding Proteins Fold. Science. 2011, 334, 517–520.

(18)

Li, B.; Fooksa, M.; Heinze, S.; Meiler, J. Finding the Needle in the Haystack: Towards Solving the Protein-Folding Problem Computationally. Crit. Rev. Biochem. Mol. Biol. 2018, 53, 1–28.

(19)

Saxena, A.; Sangwan, R. S.; Mishra, S. Fundamentals of Homology Modeling Steps and Comparison among Important Bioinformatics Tools: An Overview. Sci. Int. 2013, 1, 237–252. 33 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(20)

Rohl, C. A.; Strauss, C. E. M.; Misura, K. M. S.; Baker, D. Protein Structure Prediction Using Rosetta. In Methods in enzymology; 2004; Vol. 383, pp 66–93.

(21)

Moult, J. A Decade of CASP: Progress, Bottlenecks and Prognosis in Protein Structure Prediction. Curr. Opin. Struct. Biol. 2005, 15, 285–289.

(22)

Moult, J.; Fidelis, K.; Kryshtafovych, A.; Schwede, T.; Tramontano, A. Critical Assessment of Methods of Protein Structure Prediction: Progress and New Directions in Round XI. Proteins Struct. Funct. Bioinforma. 2016, 84, 4–14.

(23)

Compiani, M.; Capriotti, E. Computational and Theoretical Methods for Protein Folding. Biochemistry 2013, 52, 8601–8624.

(24)

Maity, H.; Maity, M.; Krishna, M. M. G.; Mayne, L.; Englander, S. W. Protein Folding: The Stepwise Assembly of Foldon Units. Proc. Natl. Acad. Sci. 2005, 102, 4741–4746.

(25)

Hu, W.; Kan, Z.-Y.; Mayne, L.; Englander, S. W. Cytochrome c Folds through Foldon-Dependent Native-like Intermediates in an Ordered Pathway. Proc. Natl. Acad. Sci. 2016, 113, 3809–3814.

(26)

Hu, W.; Walters, B. T.; Kan, Z.-Y.; Mayne, L.; Rosen, L. E.; Marqusee, S.; Englander, S. W. Stepwise Protein Folding at near Amino Acid Resolution by Hydrogen Exchange and Mass Spectrometry. Proc. Natl. Acad. Sci. 2013, 110, 7684–7689.

(27)

Haglund, E.; Danielsson, J.; Kadhirvel, S.; Lindberg, M. O.; Logan, D. T.; Oliveberg, M. Trimming Down a Protein Structure to Its Bare Foldons. J. Biol. Chem. 2012, 287, 2731–2738.

(28)

Řezáč, J.; Bím, D.; Gutten, O.; Rulíšek, L. Toward Accurate Conformational Energies of Smaller Peptides and Medium-Sized Macrocycles: MPCONF196 Benchmark Energy Data Set. J. Chem. Theory Comput. 2018, 14, 1254–1266.

(29)

Yu, W.; Wu, Z.; Chen, H.; Liu, X.; MacKerell, A. D.; Lin, Z. Comprehensive Conformational Studies of Five Tripeptides and a Deduced Method for Efficient Determinations of Peptide Structures. J. Phys. Chem. B 2012, 116, 2269–2283.

(30)

Li, H.; Lin, Z.; Luo, Y. A Fragment Based Step-by-Step Strategy for Determining the Most Stable Conformers of Biomolecules. Chem. Phys. Lett. 2014, 610–611, 303–309.

(31)

Lanza, G.; Chiacchio, M. A. Comprehensive and Accurate Ab Initio Energy Surface of Simple Alanine Peptides. ChemPhysChem 2013, 14, 3284–3293.

(32)

Hintze, B. J.; Lewis, S. M.; Richardson, J. S.; Richardson, D. C. Molprobity’s Ultimate RotamerLibrary Distributions for Model Validation. Proteins Struct. Funct. Bioinforma. 2016, 84, 1177– 1189.

(33)

Frishman, D.; Argos, P. Knowledge-Based Protein Secondary Structure Assignment. Proteins Struct. Funct. Genet. 1995, 23, 566–579.

(34)

Vymětal, J.; Vondrášek, J. Critical Assessment of Current Force Fields. Short Peptide Test Case. J. Chem. Theory Comput. 2013, 9, 441–451.

(35)

Barducci, A.; Bussi, G.; Parrinello, M. Well-Tempered Metadynamics: A Smoothly Converging and Tunable Free-Energy Method. Phys. Rev. Lett. 2008, 100, 020603.

(36)

Piana, S.; Laio, A. A Bias-Exchange Approach to Protein Folding. J. Phys. Chem. B 2007, 111, 4553–4559.

(37)

Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly 34 ACS Paragon Plus Environment

Page 34 of 38

Page 35 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435–447. (38)

Bonomi, M.; Branduardi, D.; Bussi, G.; Camilloni, C.; Provasi, D.; Raiteri, P.; Donadio, D.; Marinelli, F.; Pietrucci, F.; Broglia, R. A.; et al. PLUMED: A Portable Plugin for Free-Energy Calculations with Molecular Dynamics. Comput. Phys. Commun. 2009, 180, 1961–1972.

(39)

Lindorff-Larsen, K.; Piana, S.; Palmo, K.; Maragakis, P.; Klepeis, J. L.; Dror, R. O.; Shaw, D. E. Improved Side-Chain Torsion Potentials for the Amber Ff99SB Protein Force Field. Proteins Struct. Funct. Bioinforma. 2010, 78, 1950-1958.

(40)

Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS: A Linear Constraint Solver for Molecular Simulations. J. Comput. Chem. 1997, 18, 1463–1472.

(41)

Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577–8593.

(42)

Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kalé, L.; Schulten, K. Scalable Molecular Dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781–1802.

(43)

Hornak, V.; Abel, R.; Okur, A.; Strockbine, B.; Roitberg, A.; Simmerling, C. Comparison of Multiple Amber Force Fields and Development of Improved Protein Backbone Parameters. Proteins Struct. Funct. Bioinforma. 2006, 65, 712–725.

(44)

DeLano, W. L. PyMOL: An Open-Source Molecular Graphics Tool. CCP4 Newsl. Protein Crystallogr. 2002.

(45)

Onufriev, A.; Bashford, D.; Case, D. A. Modification of the Generalized Born Model Suitable for Macromolecules. J. Phys. Chem. B 2000, 104, 3712–3720.

(46)

Onufriev, A.; Bashford, D.; Case, D. A. Exploring Protein Native States and Large-Scale Conformational Changes with a Modified Generalized Born Model. Proteins Struct. Funct. Bioinforma. 2004, 55, 383–394.

(47)

Hamelryck, T.; Manderick, B. PDB File Parser and Structure Class Implemented in Python. Bioinformatics 2003, 19, 2308–2310.

(48)

Cock, P. J. A.; Antao, T.; Chang, J. T.; Chapman, B. A.; Cox, C. J.; Dalke, A.; Friedberg, I.; Hamelryck, T.; Kauff, F.; Wilczynski, B.; et al. Biopython: Freely Available Python Tools for Computational Molecular Biology and Bioinformatics. Bioinformatics 2009, 25, 1422–1423.

(49)

Gutten, O.; Bím, D.; Řezáč, J.; Rulíšek, L. Macrocycle Conformational Sampling by DFTD3/COSMO-RS Methodology. J. Chem. Inf. Model. 2018, 58, 48–60.

(50)

TURBOMOLE V7.2 2017, a Developement of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989-2007, TURBOMOLE GmbH, since 2007, Available at Http://Www.Turbomole.Com, (last accessed December 5, 2018).

(51)

Klamt, A.; Schüürmann, G. COSMO: A New Approach to Dielectric Screening in Solvents with Explicit Expressions for the Screening Energy and Its Gradient. J. Chem. Soc., Perkin Trans. 2 1993, 5, 799–805.

(52)

Becke, A. D. Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior. Phys. Rev. A 1988, 38, 3098–3100. 35 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(53)

Perdew, J. P. Density-Functional Approximation for the Correlation Energy of the Inhomogeneous Electron Gas. Phys. Rev. B 1986, 33, 8822–8824.

(54)

Godbout, N.; Salahub, D. R.; Andzelm, J.; Wimmer, E. Optimization of Gaussian-Type Basis Sets for Local Spin Density Functional Calculations. Part I. Boron through Neon, Optimization Technique and Validation. Can. J. Chem. 1992, 70, 560–571.

(55)

Hostaš, J.; Řezáč, J. Accurate DFT-D3 Calculations in a Small Basis Set. J. Chem. Theory Comput. 2017, 13, 3575–3585.

(56)

Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132, 154104.

(57)

Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the Damping Function in Dispersion Corrected Density Functional Theory. J. Comput. Chem. 2011, 32, 1456–1465.

(58)

Grimme, S. Supramolecular Binding Thermodynamics by Dispersion-Corrected Density Functional Theory. Chem. - A Eur. J. 2012, 18, 9955–9964.

(59)

Weigend, F.; Ahlrichs, R. Balanced Basis Sets of Split Valence, Triple Zeta Valence and Quadruple Zeta Valence Quality for H to Rn: Design and Assessment of Accuracy. Phys. Chem. Chem. Phys. 2005, 7, 3297–3305.

(60)

Eckert, F.; Klamt, A. Fast Solvent Screening via Quantum Chemistry: COSMO-RS Approach. AIChE J. 2002, 48, 369–385.

(61)

Klamt, A. Conductor-like Screening Model for Real Solvents: A New Approach to the Quantitative Calculation of Solvation Phenomena. J. Phys. Chem. 1995, 99, 2224–2235.

(62)

Klamt, A.; Jonas, V.; Bürger, T.; Lohrenz, J. C. W. Refinement and Parametrization of COSMORS. J. Phys. Chem. A 1998, 102, 5074–5085.

(63)

Takano, Y.; Kusaka, A.; Nakamura, H. Density Functional Study of Molecular Interactions in Secondary Structures of Proteins. Biophys. Physicobiology 2016, 13, 27–35.

(64)

Valdés, H.; Řeha, D.; Hobza, P. Structure of Isolated Tryptophyl-Glycine Dipeptide and Tryptophyl-Glycyl- Glycine Tripeptide: Ab Initio SCC-DFTB-D Molecular Dynamics Simulations and High-Level Correlated Ab Initio Quantum Chemical Calculations. J. Phys. Chem. B 2006, 110, 6385–6396.

(65)

Brenner, V.; Piuzzi, F.; Dimicoli, I.; Tardivel, B.; Mons, M. Spectroscopic Evidence for the Formation of Helical Structures in Gas-Phase Short Peptide Chains †. J. Phys. Chem. A 2007, 111, 7347–7354.

(66)

Lanza, G.; Chiacchio, M. A. Effects of Hydration on the Zwitterion Trialanine Conformation by Electronic Structure Theory. J. Phys. Chem. B 2016, 120, 11705–11719.

(67)

Novotny, M.; Kleywegt, G. J. A Survey of Left-Handed Helices in Protein Structures. J. Mol. Biol. 2005, 347, 231–241.

(68)

Minor, D. L.; Kim, P. S. Context Is a Major Determinant of β-Sheet Propensity. Nature 1994, 371, 264–267.

(69)

Merkel, J. S.; Sturtevant, J. M.; Regan, L. Sidechain Interactions in Parallel β Sheets: The Energetics of Cross-Strand Pairings. Structure 1999, 7, 1333–1343. 36 ACS Paragon Plus Environment

Page 36 of 38

Page 37 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(70)

Kuroda, Y.; Hamada, D.; Tanaka, T.; Goto, Y. High Helicity of Peptide Fragments Corresponding to β-Strand Regions of β-Lactoglobulin Observed by 2D-NMR Spectroscopy. Fold. Des. 1996, 1, 255–263.

(71)

Kuwata, K.; Shastry, R.; Cheng, H.; Hoshino, M.; Batt, C. A.; Goto, Y.; Roder, H. Structural and Kinetic Characterization of Early Folding Events in β-Lactoglobulin. Nat. Struct. Biol. 2001, 8, 151–155.

37 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

TOC Graphics

38 ACS Paragon Plus Environment

Page 38 of 38