Dynamic Solvent Accessible Surface Areas of Proteogenic a-Amino Acids as a Function of Explored f,y Dihedral Angle Space Modelled using the Ac-Ala-Xaa-Ala-NH2 system and CHARMM36m Force Field

1. Introduction

The study of proteins and their function is a cornerstone in molecular biology, with the implications for the understanding of normal physiology, disease mechanisms and custom pharmaceutical design. Techniques such as X-ray crystallography, NMR, and cryo-EM have proven effective for resolving protein structures yet at the same time are of much higher initial upfront cost associated with the instrumentation and the expertise needed to interpret and build a model of the experimental results.[1,2,3,4,5,6] Understanding the mechanisms and thermodynamics of protein folding and unfolding are key to the refinement of both deterministic methods such as Molecular Dynamics (MD) as well as artificial intelligence driven algorithms such as AlphaFold2 (AF2) and AlphaFold3 (AF3).[7,8,9,10,11] During protein folding, sequestration of hydrophobic surface areas within the folded core as well as the exposure of polar and charged surface areas at the protein/solvent (water) interface take place, which are many times referred to as hydrophobic interactions. These processes have also been shown to drive the interactions at protein-protein interfaces and likely play a role in protein-ligand interactions. The thermodynamics of these “interactions” can be computationally approximated since the changes in solvent accessible surface area (SASA) of amino acids are proportional with the transfer free energies from the bulk solvent (water) to a hydrophobic solvent such as ethanol, octanol, or cyclohexane.[12]

Initial studies by Lee and Richards used the atomic coordinates of extended conformation of the Ala-Xaa-Ala or Gly-Xaa-Gly model peptides, where Xaa represents the amino acids of interest.[13] The validity of this approach, however, has been questioned due to the unrealistic conformation of the tri-peptides and their relationship to the folded or even unfolded state of a protein. Furthermore the Gly-Xaa-Gly model peptide does not provide a realistic model of SASA since in a protein, the Xaa residue is affected by its nearest neighbors which on average are bulkier than Gly which has a volume of 0.0638 nm³ compared to the mean volume of 0.1455±0.0418 nm³ (Max: 0.2317 nm³, Trp) for the 20 proteogenic α-amino acids.[14,15] If the frequency of amino acids within expressed proteins is taken into consideration, the weighted mean is 0.1383±0.0389 nm³.[16] Extended static models as proposed by Lee and Richards also cannot account for the dynamic nature of protein/peptide structures and the sampling of the ϕ,ψ space that defines protein secondary structure, Figure 1 and Figure 2.

The first study to account for dynamic sampling of a model peptide to calculate SASAs was completed by Zielenkiewicz and Saenger.[17] The authors used 1 ns, 368 K, in-vacuo MD simulations to sample SASAs of the MeNH-Ala-Xaa-Ala-Me model peptides. They used several additional solvent models to confirm the validity of the results by comparison to the in-vacuo calculations. Elevated temperature was used to account for the unfolded state. The study was limited by its use of the consistent valence force field (CVFF), which is not well parameterized for amino acids, non-physiological temperature, use of neutral amino acids only without counterions, and lack of explicit solvents.[18,19] These factors have been shown to effect simulation results and play an integral role in peptide/protein conformation, folding and dynamic solvent interactions. Despite the use of MD sampling, only mean SASAs were reported without standard deviations. The lack of reported ϕ,ψ dihedral angle space sampling makes assessment of simulation quality as well as conformational analysis of the results not possible.

SASAs of amino acids have also been calculated using high quality X-ray crystal structures from the protein Databank (PDB) combined with computationally constructed Gly-Xaa-Gly model peptides.[1,20] A set of 3197 high quality protein crystal structures was used to collect both SASA as well as ϕ,ψ dihedral angle data. The ϕ,ψ dihedral angles were then binned into a three-dimensional histogram with 5° of resolution for each angle and assigned to regions based on population size. The core region accounted for 80% of sampled data, the allowed region accounted for 97% of sampled data, and the generous region extended the allowed region by 20° in all directions. The maximum observed SASA was defined for each region. Then, the Gly-Xaa-Gly model peptides were used to systematically sample the ϕ,ψ dihedral angles with a rotational algorithm with 1° of resolution while relaxing the side chain χ dihedral angles to prevent unfavorable contacts/interactions with the side chain. Using a combination of observed and empirical calculations, a theoretical maximum for the core, allowed, generous, and all regions was determined. The study by Tien et al. accounts for some of the favorable sampling of regions of ϕ,ψ space, Figure 2, but does not clearly delineate regions into the descriptors assigned to secondary structure motifs: β, α, α_L, ε, and contiguous. A visual inspection of these regions as show in Figure 3 clearly demonstrates the effect of the ϕ,ψ dihedral angles on both the backbone and side chain components of the central Xaa residue of the Ac-Ala-Xaa-Ala-NH₂ system. The work also does not account for the fact that in solution, an ensemble of conformations exists within the ϕ,ψ dihedral angle space (ϕ,ψ space), sampling a range of energetically favorable dihedral angles and a resulting range of sampled SASA values.

Topham et al. used quantum mechanical/density functional theory (QM/DFT) to do gas phase energy minimization and asses SASA for the Ac-Ala-Xaa-Ala-NH₂ model peptides as a function of β-sheet and α-helical regions of ϕ,ψ space for the central Xaa residue.[21] QM/DFT provides a better estimate of molecular geometry then previous MD simulations or empirical studies. While the SASA calculations are accurate they do not account for the dynamic behavior of the peptide in solution and cannot be used to determine a comparison of the separate populations due to the lack of sufficient sampling numbers within ϕ,ψ space as well as the use of gas phase states of extended tripeptides do not represent physiological conditions.

Herein, we report the backbone, side chain, and whole residue SASAs of the Xaa residues within the Ac-Ala-Xaa-Ala-NH₂ model peptides as a function of the Xaa residues MD-sampled ϕ,ψ space. The sampled ϕ,ψ space of each Xaa residue was grouped in the β, α, α_L, ε, and contiguous regions of ϕ,ψ space using a density clustering algorithm. The peptide systems were simulated using the CHARMM36m force field with associated TIP3Pm explicit solvation model of water as implemented within GROMACS 2022.[22,23] Non-bonded interactions were calculated using the current updated electrostatic and Leonard Jones particle mesh Ewald summation parameters within the constant temperature and pressure ensemble at a physiologically relevant temperature, 310 K (37° C).[24,25] A robust statistical analysis was performed for each sampled population associated ϕ,ψ space of the central Xaa residue, and resulting backbone, side chain, and whole residue SASAs.

2. Results

Comprehensive figures and tables with associated statical analysis are provided in the Supplementary Materials: Figures S1 and S2; Tables S1 through S40.

2.1. Systems Equilibration.

Calculating systems densities and -N-Cα-C- backbone entropies as functions of time showed that all simulations of Ac-Ala-Xaa-Ala-NH₂ and (Ac-Ala-Cys-Ala-NH₂)₂ model systems plateaued within the first 0.2 µs. Trajectories were sampled between 0.2 µs and 1 μs in 2 ps steps resulting in 400,000 sampled conformations for each peptide system. The one exception to this was the (Ac-Ala-Cys-Ala-NH₂)₂ system because of its two Ac-Ala-Cys-Ala-NH₂ peptide chains, 800,000 sampled conformations were analyzed utilizing the central Cys residue ϕ,ψ dihedral angle and SASA information from both chains.

2.2. ϕ, ψ Space Analysis (β, α, α_L, ε, and Contiguous Regions).

Two dimensional histograms of the log₁₀-scaled probabilities (ρ) of the sampled ϕ,ψ dihedral angles of the central Xaa residues for the Ac-Ala-Xaa-Ala-NH₂ and (Ac-Ala-Cys-Ala-NH₂)₂ systems are shown in Figure S1. The histograms are overlayed onto Xaa residue specific probability contours representing 98% and 99.8% of the sampled populations.[26] The smallest possible sampled population size for density clustering was 1000 conformations within 10° which represents 0.25% of the total sample size (400,000) except for the (Ac-Ala-Cys-Ala-NH₂)₂ system where it represents 0.125%. Selection of these regions of ϕ,ψ space (β, α, α_L, ε, and contiguous) from the cluster analysis was based on the quality of the derived clusters adhering to the above listed regions of ϕ,ψ space as determined by visual inspection of the clusters and the goal of having the defined regions (β, α, α_L, and,ε) account for 98% of the sampled population with the contiguous region accounting for the remaining 2% (Figures S1 and S2 and Tables 1 and S1).

The contours representing the sampled ϕ,ψ space as well as the histogram probabilities vary between Xaa residues as a function of side chain size (i.e., Ala, Val, and Leu), side chain steric hindrance, (i.e., Leu compared to Ile) as shown in Figure 4, chemical properties across the amino acid sequences (acidic, basic, hydrophobic, aromatic, and polar), and charge (i.e., His:ND1 and ND2 compared to His⁺, Asp^- compared to neutral Asp:H, and neutral LysN compared to Lys⁺), as shown in Figure 5 and Figure 6.

The β region of ϕ,ψ space is preferentially sampled for all model systems, Table 1. The Gly residue has the lowest probability of β ϕ,ψ space, ρ = 0.4917 while Val and Pro:trans are the highest, ρ = 0.9086 and ρ = 0.9204, respectively. The second most sampled region of ϕ,ψ space is α with Leu, Asn, Gln, Asp:H, His⁺, Cys^- and Pro:cis residues having the highest probability of sampling this conformation; ρ = 0.2020 to ρ = 0.2422. The α_L region is the third most populated with most residues having sampling probabilities ranging from ρ = 0.01 to ρ = 0.05. There are outliers: His:NE2 (ρ = 0.1219), His⁺ (ρ = 0.1073), Asn (ρ = 0.1002), and Asp:H (ρ = 0.0806) being slightly more favorable that other residue types. The Cys^-, Pro:cis and Pro:trans residues do not sample the α_L region. There also appear to be significant differences in the probability of sampling certain regions of ϕ,ψ space especially the ε region which is more favorable to Gly but can also be sampled by the aromatic residues His:ND1, His:NE2, Phe, and Tyr but not Trp; His⁺ and Tyr^-; the hydrophobic residue Ala, the polar residues Ser, Asn, Asp:H, and Lys^N; and the negatively charged residue Asp. Therefore, it appears, that side chain size may play some role in the propensity to sample the ε region since bulky hydrophobic residues: Val, Leu, Ile, and Met do not favor it. Aromatic residues sample small populations of the ε region except for the bulkier Trp. Sampling of the ε region represents a complex relationship between side chain size, steric hindrance, and chemical properties.

Gly reveals that, despite its lack of side chain and increased flexibility, equal sampling of the β, α, α_L, and ε regions does not occur. The most favored region is β, ρ = 0.4917 while the second is ε, ρ = 0.3596. The α and α_L regions are sampled with probabilities of ρ = 0.0430 and ρ = 0.0586, respectively. The results differ somewhat from the observed ϕ,ψ plots generated from structural databases where the α_L conformation is more heavily favored but also from the ϕ,ψ plots generated from QM/MM simulations of Ace-Gly-NH₂ where sampling within the regions is more equal.[27,28] These difference suggest that sampling of ϕ,ψ space in Gly at least within the current Ac-Ala-Gly-Ala-NH₂ system is significantly influenced by the adjacent i-1 and i+1 amino acid residues.

Pro:cis is mostly β, ρ = 0.7959, with the next most common conformation α, ρ = 0.2104, the sampling of contiguous space is extremely small, ρ = 0.0037. Pro:trans is mostly β, ρ = 0.9204, with a small amount of α, ρ = 0.0328, and contiguous space, ρ = 0.0469, as shown in Figure 7. Neither Pro:cis nor Pro:trans sample the α_L or ε regions sufficiently to be assigned a population by cluster analysis. A ${\chi }^2$ analysis of Pro:cis and Pro:trans demonstrates statistically significant differences comparing β, α, and contiguous conformations sampled by each system: β (${\chi }^2$ = 28841, p-value < 0.0001), α (${\chi }^2$ = 59064, p-value < 0.0001) and contiguous (${\chi }^2$ = 15144, p-value < 0.0001). This demonstrates the importance of the ω dihedral angle that controls the cis/trans relationship of the Pro residue and its effects on ϕ,ψ space sampling.[29] Umbrella sampling QM/MM studies on a model Ac-Pro-Nme peptide in explicit solvent demonstrated that the trans state is more stable by approximately 4 kcal/mol with an energy barrier of approximately 20 kcal/mol separating the cis and trans states.[30]

The (Ac-Ala-Cys-Ala-NH₂)₂ has added complexity of non-bonded interactions between the two peptide chains affecting the sampled ϕ,ψ space of the central Cys residues. Like Cys:H and Cys^-, β conformation is favored: Cys-Cys, ρ = 0.7406; Cys:H, ρ = 0.7795; and Cys-, ρ = 0.7757. Despite similarities in ρ, statistically significant differences exist (${\chi }^2$ = 3006 p-value < 0.0001). These differences are also present for the α conformation: Cys-Cys, ρ = 0.1690; Cys:H, ρ = 0.1359; and Cys^-, ρ = 0.2134 (${\chi }^2$ = 8518 p-value < 0.0001). For the α_L conformation: Cys-Cys, ρ = 0.0449; Cys:H, ρ = 0.0445; and Cys^-, ρ = N/S (${\chi }^2$ = 1.248 p-value = 0.2639) these differences are not statistically significant. Cys-Cys, Cys:H and Cys^- do not significantly sample the ε region of ϕ,ψ space.

Statistically significant differences exist across all populations of ϕ,ψ space regions sampled by the Ac-Ala-Xaa-Ala-NH₂ peptides: β, α, α_L, ε, and contiguous with p-values < 0.0001, Table S1. A Marascuillo procedure analysis, with exclusion of Gly, Pro:cis, Pro:trans, and Cys-Cys, ϕ,ψ space for a total of 27 Ac-Ala-Xaa-Ala-NH₂ model systems were performed. This resulted in 351 pairwise comparisons for the β, α, and contiguous regions, 325 for the α_L region, and 45 for the ε region (Tables S2 − S11). For the β region, 330 of 351 pairwise comparisons had statistically significant differences. The α region, 322 of 351 pairwise had statistically significant differences. The contiguous region, 275 of 351 pairwise comparisons had statistically significant differences. The α_L region, 301 of 325 pairwise comparison had statistically significant differences; and the ε region, 44 of 45 pairwise comparisons had statistically significant differences. The pairwise comparisons demonstrate interesting relationships within the Xaa amino acids. Statistically significant differences occur in the β region population probabilities of His:ND1, His:ND2, and His⁺, Tyr to Tyr^-, Arg:NE to Arg:NH, Arg:NH to Arg, Asp:H to Asp, Glu:H to Glu, and Lys to Lys^N. The differences between Arg:NE to Arg and Cys:H to Cys^- are not statistically significant. For the α region, statistically significant differences are noted for: His:ND1 to His:NE2, His:ND1 to His⁺, His:NE2 to His⁺, Tyr to Tyr^-, Cys:H to Cys^-, Arg:NE to Arg:NH, Arg:NE to Arg, Arg:NH to Arg, Asp:H to Asp, Glu:H to Glu, and LysN to Lys. The same relationships hold true for the α_L conformation except for Cys:H to Cys^- since Cys^- does not sample this region of ϕ,ψ space sufficiently to be assigned a population density. The ε region is more sparsely populated with few amino acids having sufficient population density to be analyzed. Significant differences exist between His:ND1 to His:NE2 and Asp:H to Asp. Previous relationships as noted for β, α, and α_L are not noted secondary to regional sampling. These findings suggest that subtle differences in amino acid charge and tautomer state (His and Arg residues) can affect the probability of sampling a particular region of ϕ,ψ space.

2.3. ϕ, ψ Dihedral Angle Analysis.

With assignment of sampled ϕ,ψ space of the Xaa residues into their respective β, α, α_L, ε, or contiguous region populations; mean and standard deviations for each associated ϕ and ψ dihedral angle of the Xaa residue were calculated. Results for each region are shown in Table S12. Gly has a mean ϕ dihedral angle of the β region that is slightly more positive with a decreased standard deviation compared to other residues while the mean ψ dihedral angle that is slightly more negative with a slightly larger stand deviation compared to other residues. The α region of Gly is more like other residues with respect to the mean ϕ dihedral and its standard deviation while the mean ψ dihedral and its standard deviation are more variable compared to all residues. This relationship of variability is true for the ϕ,ψ dihedral angles for both the α_L and ε regions. The ϕ,ψ dihedral angles for the β, α, and contiguous regions of Pro:cis and Pro:trans were compared using a Welch’s t-test.[31] There are statistically significant differences for both ϕ and ψ dihedral angles for all regions: β, ϕ (t = -108.5433, p-value < 0.0001) ψ (t = 95.2015, p-value < 0.0001); α, ϕ (t=248.7765, p-value < 0.0001) ψ (t = 284.1296, p-value < 0.0001); and contiguous, ϕ (t = -7.0356, p-value < 0.0001) ψ (t = -13.3615, p-value < 0.0001). Neither the α_L nor ε regions are sampled by Pro:cis or Pro:trans. Cys:H, Cys⁻ and Cys-Cys were compared using the Welch’s ANOVA. Statistically significant differences for both the ϕ and ψ dihedral angles for the following conformations: β, ϕ (F = 135590.27, p < 0.0001) ψ (F = 7271.32, p-value < 0.0001); α, ϕ (F= 58329.27, p-value < 0.0001) ψ (F = 45566.36, p-value < 0.0001); and contiguous, ϕ (F = 1482.19, p-value < 0.0001) ψ (F = - 503.03, p-value < 0.0001). The differences for the ϕ dihedral angle of the α_L conformation are statistically significant (F = 43.75, p-value < 0.0001) while the ψ dihedral angles are not (F = 503.03, p-value = 1.0).

Comprehensive results of the Scheffe test are provided in Supplemental Table S13h S32 (Gly, Pro:cis, Pro:trans, and are Cys-Cys excluded). For the ϕ dihedral angles of the β region, 322 of 351 pairwise comparisons are statistically different. For the ψ dihedral angles of the β region, 330 of 351 pairwise comparisons are statistically different. For the ϕ dihedral angles of the α region, 326 of 351 pairwise comparisons are statistically different. For the ψ dihedral angles of the α region, 324 of 351 pairwise comparisons are statistically different. For the ϕ dihedral angles of the α_L region, 222 of 325 pairwise comparisons are statistically different. For the ψ dihedral angles of the α_L region, 245 of 325 pairwise comparisons are statistically different. . For the ϕ dihedral angles of the ε region, 36 of 45 pairwise comparisons are statistically different. For the ψ dihedral angles of the ε region, 43 of 45 pairwise comparisons are statistically different.

Like the pairwise population analysis discussed above, statistically significant differences in the ϕ,ψ dihedral angles for the β, α, α_L and ε regions are also demonstrated. For the β region, differences exist in the ϕ dihedral angle for: His:ND1 to His:NE2, His:ND1 to His⁺, His NE2 to His⁺, Tyr to Tyr^-, Cys:H to Cys^-, Arg:NE to Arg:NH, Arg:NE to Arg, Arg:NH to Arg, Asp:H to Asp, Glu:H to Glu and Lys^N to Lys. These differences also occur for the ψ dihedral angle with the exceptions of: His:ND1 compared to His⁺ and His:NE2 compared to His⁺. For the α region, the above listed differences exist in the ϕ dihedral angle as they do for the β region. For the corresponding ψ dihedral angle, the exception is His:ND1 to His:NE2. For the α_L region, differences exist for the ϕ dihedral angle with the exceptions of His:ND1 to His:NE2 and Arg:NE to Arg:NH. For the corresponding ψ dihedral angle, the exceptions are: Tyr to Tyr-, Arg:NE to Arg:NH, Arg:NE to Arg, and Arg:NH to Arg. The Cys^- system does not adequately sample the α_L region for a Cys:H to Cys^- interaction to be analyzed. Cys:H, Cys^-, Arg:NE, Arg:NH, Arg, Glu:H, Glu, and Lys do not adequately sample the ε region to perform a pairwise analysis. All other pairwise interactions demonstrate significantly differences except for His:ND1 to His:NE2 for the ψ dihedral angle. Like the population analysis, these results suggest that subtle differences in amino acid charge and tautomer state (His and Arg residues) can affect the sampling regions of ϕ,ψ space particularly the β and α regions.

2.4. Hydrogen Bond Analysis.

The possible role of hydrogen bond stabilization of secondary structure conformation was evaluated by measuring the i-1 (Ala 1) carbonyl oxygen to i+1 (Ala 3) amide nitrogen and the i-1 (Ala 1) amide nitrogen to i+1 (Ala 3) carbonyl oxygen distances of the Ac-Ala-Xaa-Ala-NH₂ peptides. Distances between 0.27 and 0.33 nm were considered consistent with hydrogen bond formation.[32] The global mean ± standard deviation donor-acceptor distances, mean ± standard deviation donor-acceptor distance of identified hydrogen bonds, probability of a hydrogen bond being present (ρ), and number (n) of hydrogen bonds identified for each Ac-Ala-Xaa-Ala-NH₂ peptide are summarized in Tables S33 and S34. The probability of hydrogen bonds between either the i-1 (Ala 1) carbonyl oxygen to i+1 (Ala 3) amide nitrogen or the i-1 (Ala 1) amide nitrogen to i+1 (Ala 3) carbonyl oxygen are low. Values range from ρ = 0.0000 (n = 2) for the Pro:cis residue to ρ = 0.0868 for the Pro:trans residue for the i-1 (Ala 1) carbonyl oxygen to i+1 (Ala 3) amide nitrogen and from ρ = 0.0000 (n = not sampled) for the Gly, Val, Leu, Phe, Ser, Thr, Cys:H, Arg, Lys, Glu, and Pro:trans residues to ρ = 0.0011 (n = 456) for the Pro:cis residue for the i-1 (Ala 1) amide nitrogen to i+1 (Ala 3) carbonyl oxygen. Formation of hydrogen bonds between the i-1 (Ala 1) amide nitrogen to i+1 (Ala 3) carbonyl oxygen are rarer with only 20 of the 31 peptides simulated forming this type of hydrogen bond. The data suggests the hydrogen bonds do not play a significant role in the stabilization of conformations sampled within the ϕ,ψ space regions.

2.5. Solvent Accessible Surface Area Analysis

The SASAs for the whole, backbone, and side chain components of the Xaa residues as a function of β, α, α_L, ε and contiguous regions of φ,ψ space are given in Tables 2 and Tables S35, S36, and S37, respectively. Statistically significant differences exist across the φ,ψ space regions: β, α, α_L, ε, and contiguous for each residue with p-values < 0.0001. Pairwise comparisons using the Scheffe test are provided in Tables S38 through S40. Considering the Xaa residues of the Ac-Ala-Xaa-Ala-NH₂ systems as whole residues (backbone and side chain combined), Table 2, the most solvent shield (lowest SASA) region of φ,ψ space is the β region except for Ile which is α_L, Phe which is ε, and Trp which is α_L. The most solvent accessible (highest SASA) φ,ψ space region is α for: Gly, Ala, Val, Leu, Ile, Trp, Thr, Cys:H, Arg, Cys^-, Pro:cis ,and Pro:trans; while it is α_L for: Met, His:ND1,His:NE2, Phe, Tyr, Ser, Asn, Gln, Arg:NE, Arg:NH, Asp:H, Glu:H, Lys^N, His⁺, Lys, Asp, Glu, Tyr^-, Cys-Cys. The contiguous region of φ,ψ space tends to have SASA values that are between the low values of the β region and the high levels of α and α_L. Pairwise comparisons across all φ,ψ space regions indicated statistically significant differences between all regions with respect to SASA with the follow exceptions: Val, β - α; His:ND1, ε - contiguous; Phe, β - ε; Trp, β - α_L; Ser, α - α_L; Thr, β - α_L; Asn, β - ε and α to α_L; Glu:H, α - α_L, Arg, α to α_L; and Asp, α to α_L; Table S38.

The results can be deconvoluted into separate backbone and side chain values, Tables S36 and S37. For the backbone, the φ,ψ space region corresponding to the lowest SASA is the β region for all residues. The φ,ψ space region corresponding to the highest SASA values is α_L with the exceptions of: Val, Trp, and Thr, which are contiguous and Cys-, Cys-Cys, Pro:cis, and Pro:trans which are a. Pairwise comparisons across all φ,ψ space regions indicated statistically significant differences between all conformations with respect to SASA with the follow exceptions: Val, α_L – contiguous_; Leu, α – contiguous; Ile, α – α_L; Tyr, α – ε; Trp, α – α_L; Thr, α – α_L; Cys, α – contiguous; Asn, ε – contiguous; Arg:NH, α - contiguous; Glu:H, α - contiguous; Lys^N, α – ε and ε – contiguous; Table S39. The pattern of least and greatest SASA does not hold however for the side chains where the φ,ψ space region with the highest SASA values is a with the exceptions of: His⁺, and Cys-Cys which are α_L and Pro:cis, and Pro:trans which are β. The SASA values in descending order tend to be much more variable with respect which φ,ψ space region is the lowest. Pairwise comparisons across all φ,ψ space regions indicated statistically significant differences between all conformations with respect to SASA with the follow exceptions: Ala, α_L – contiguous; Val, β – contiguous; Lys^N, β – α_L; Arg, β – α_L; His⁺, α – α_L; Lys, β – α_L; Cys-Cys, α – contiguous; Table S40.

4. Discussion

Solvent accessible surface area (SASA) of amino acid residues plays a significant role in inter-residue, protein-ligand, and protein-protein interactions. Therefore, we studied backbone, side chain, and whole residue SASAs of the Xaa residues within the Ac-Ala-Xaa-Ala-NH₂ model peptides using MD simulations. Our hypothesis was that rather than existing in extended conformations, the Ac-Ala-Xaa-Ala-NH₂ model peptides central Xaa residues’ ϕ,ψ dihedral angles would sample the β, α, α_L, ε, and contiguous regions of ϕ,ψ space in a manner specific to the Xaa residue type, effecting the resulting SASAs. We have demonstrated this using a comprehensive statistical analysis of the simulation data.

The Ac-Ala-Xaa-Ala-NH₂ model is a compromise for the determination of SASA values of the central Xaa residue in the unfolded state. The two most common systems used for calculation of these values have been variations of -Gly-Xaa-Gly- and -Ala-Xaa-Ala- sequences with different N- and C-terminal groups being used to maintain charge neutrality of the systems.[13,14,17,20,21] The main criticism with regards to use of -Gly-Xaa-Gly- sequence is that the Gly residues do not have side chains and as a consequence, do not realistically restrict sampling of the central Xaa residues ϕ,ψ dihedral angles.[17,20,21] It is well known that both the i-1 and i+1 residues surrounding the central Xaa residue can have an effect on it’s available sampling of ϕ,ψ space and if the peptide is long enough, effect the resulting secondary structure, potential tertiary structure, and resulting SASA values. An ideal situation would be a comprehensive understanding of the effects of every combination of amino acid surrounding the central residue of interest. For a three-residue system, considering all the potential combinations of charge and tautomers, this would represent a total of 29,791 simulations. The -Ala-Xaa-Ala- model system is therefore seen as a reasonable compromise for the present study.

An important consideration when interpreting the ϕ,ψ space of the Ac-Ala-Xaa-Ala-NH₂ data presented herein is that the central Xaa residue occupying a particular region of ϕ,ψ space (β, α, α_L, ε, or contiguous) does not mean that the Ac-Ala-Xaa-Ala-NH₂ is conforming to a particular secondary structure. Secondary structure elements (parallel and antiparallel β-sheet, α-helix, 3₁₀-helix, π-helix, turns, and bends) as observed in NMR, X-ray crystal, and cryo-EM published structures must adhere to sets of defined rules with respect to the ϕ,ψ dihedral angles of the involved residues as well as hydrogen bonding patterns and strand orientation.[33] The minimal secondary structure, with respect to number of amino acid residues involved (not counting the terminal Ac- or -NH₂), that the Ac-Ala-Xaa-Ala-NH₂ peptides could form would be either the γ-turn or γ’-turn.[34,35] The γ-turn or γ’-turn are defined by a hydrogen bond between carbonyl oxygen of residue i-1 and amide hydrogen of residue i+1. The ϕ,ψ dihedral angles of central residue i (Xaa in this case) are (65° to 75°, -55° to -65°) and (-65° to -75°, 55° to 65°) for the γ-turn or γ’-turn, respectively. Significant populations within the ϕ,ψ space of the Xaa residues representing either γ-turn or γ’-turn are not observed in this study, Figures S1 and S2. We confirmed the absence of stable hydrogen bonds between the i-1 and i+1 residues by measuring the distance between the carbonyl oxygen of the Ala 1 and the amide nitrogen of Ala 3, Table S33. We also examined the possibility of atypical hydrogen bonding between the i-1 and i+1 residues measuring the distance between the amide nitrogen of Ala 1 and the carbonyl oxygen of Ala 3, Table S34). The data indicate that the Ac-Ala-Xaa-Ala-NH₂ peptides exist in solution as an ensemble of conformations, not stabilized by hydrogen bonds or conforming to a defined class of secondary structures, with the central Xaa residue favoring the β region of ϕ,ψ space. These findings are consistent with the current parameterization of the CHARMM36m force field.[22]

Our SASA results can also be compared to those previously published by Zielenkiewicz and Saenger, and Topham and Smith.[17,21] Zielenkiewicz and Saenger used the atomic radii of Shrake and Rupley for their main results (whole, backbone, and side chain) but also published comparison values obtained by using the atomic radii of Rose et al. and Miller et al. for the whole residue.[14,36,37] Their work shows differences between the mean SASA of the whole residues and the previously published results of the other above listed authors. Much of these differences may be attributed to the sampling technique (MD) used by Zielenkiewicz and Saenger compared to the static conformation techniques of the other authors since the differences between the calculated SASA values from the MD simulations using the different atomic radii are smaller than the differences between those obtained by MD simulations and the previously published values.[17] Unfortunately, Zielenkiewicz and Saenger did not publish associated standard deviations for each amino acid so that a more robust statistical comparison could be performed.

Comparing the global mean SASA data of the whole residues in this study to the data reported in Zielenkiewicz and Saenger demonstrates that for the atomic radii of Shrake and Rupley, the absolute differences between values range from 0.009 nm² for Cys:H to 0.135 nm² for Ile. These represent percentage differences range from 0.473% for Phe to 8.21% for Val.[17,36] Comparing to the values of Rose et al. shows absolute differences ranging from 0.005 nm² for Ala to 0.173 nm² for Trp. Percentage differences range from 0.329% for Lys to 6.88% for Trp.[17,37] Comparing the values of Miller et al. shows absolute differences ranging from 0.003 nm² for Glu to 0.135 nm² for Trp. Percentage difference range from 0.169% for Glu to 6.78% for Pro.[14,17] Similar comparisons are provided for the Topham and Smith. data which utilized the atomic radii of Chothia et al. [21,28] The absolute differences in SASA range from 0.0028 nm² for Glu to 0.112 nm² for Pro:cis. Percentage differences range from 0.158% for Glu to 7.59% for Pro:cis. Care should be taken with respect to interpreting these values since they represent the SASA of a geometry-optimized extended conformation of the Ac-Ala-Xaa-Ala-NH₂ peptides using a B3LYP level of quantum chemical calculation, not dynamic sampling of conformations from MD.

The whole residues all conformations SASA values in Table 2 for each Xaa residue can be compared statistically to those published in Zielenkiewicz and Topham by making some simple assumptions. These assumptions are based on those used within the Welch’-Aspen t-test.[31,39] Since standard deviations of the mean were not published in these two works, an estimated standard deviation is calculated from the global mean and the standard deviation is the square root of the sum of the variances between each published value and the global mean. Assuming equal population sizes across all samples are equal, the values can be compared using the t-test.[40] Absolute differences range from 0.003 nm² for Ala to 0.205 nm² for Trp. Percentage differences range from 0.274% for Ala to 5.25% for Ile. The differences between all values for the whole residues reported in Table 2 are statistically different compared to the means calculated from the previously published results with p-values <0.0001.

There are multiple reasons these differences exist within the SASA data. First, as demonstrated by reviewing of the manuscripts of both Zielenkiewicz and Topham; there are multiple atomic radii that can be used to calculate SASA, including those of Tsai et al. which were utilized herein.[35,36,37,38,41] Each set of radii have subtle differences in value for each atom type and will contribute to different calculated SASA values. Second, although the surface area calculation algorithm of Lee and Richards used in Zielenkiewicz is the same as the one used herein; we utilized a sampling density of 162000 points per sphere compared to the 987 points per sphere which should improve that accuracy of the measurement. Third, there are differences in methodology between the study of Zielenkiewicz and the methodology used herein. As stated previously, the study by Zielenkiewicz utilized the CVFF force field for their simulations which were done in vacuo with a dielectric constant of 80 for charge screening/solvent simulation using a charge neutral amino acid model without counterions, an elevated temperature of 368 K and pairwise calculation of long range interactions (partial charges and Leonard-Jones potentials) with a short cut off distance of 1 nm and a switching function starting at 0.85 nm [17,18] We utilized the CHARMM36m force field with explicit TIP3Pm solvation and counterions, a more physiological temperature of 310 K, and particle mesh Ewald summations of long-distance interactions with a cut off distance of 1.2 nm which should result in a more accurate sampling of relevant conformational ensembles. Forth, is data reporting and analysis. Zielenkiewicz reports a mean SASA without a standard deviation. They also do not report the sampled conformational ensembles. Both these factors make a direct comparison of the data using robust statistical tools impossible.

5. Materials and Methods

This study of the central Xaa residue within the Ac-Ala-Xaa-Ala-NH₂ system was limited to the 20 encoded proteogenic α-amino acids, their tautomers (His:ND1, His:NE2, Arg:NE, and Arg:NH) and their acidic and basic forms. The following model peptides were constructed in their extended conformations using Pymol v. 2.30.[42]

Neutral Termini Tripeptide

Ac-Ala-Xaa-Ala-NH₂

and

Neutral internal Cys-Cys bond (Ac-Ala-Cys-Ala-NH₂)₂

All MD simulations were performed with the GROMACS 2022 software package using the CHARMM36m force field parameters with the CHARMM36m consistent version of TIP3Pm water. [22,23,43,44,45,46,47,48,49,50,51,52,53] Peptides were solvated in dodecahedral boxes with TIP3Pm water with 150 mM NaCl.[19] Additional Cl^– and Na⁺ ions were used to neutralize the charges of the systems. The minimum distance of the peptide to the edge of the dodecahedron was 1.4 nm. The solvated conformations were prepared for simulation using a modified protocol of Roe et al. [54] The particle mesh Ewald method was used to calculate long-range electrostatic and van der Waals interactions with cutoff distances of 1.2 nm, and a Fourier spacing was 0.15 nm.[24,25] Initial velocities were assigned to the solvent and peptide separately using a temperature dependent Maxwell-Boltzmann distribution unless otherwise specified.[55,56] The simulation protocol is given below:

• Energy-minimized with 5000 steps of steepest descent with strong positional restraints to the heavy atoms of the peptide, 209.29 kJ/mol/nm with the initial coordinates as a reference.
• 1000 ps of constant volume and temperature (NVT) MD with a weakly coupled Berendsen thermostat with τ=0.5 ps, a 1 fs time-step and strong positional restraints to the heavy atoms of the peptide, 209.29 kJ/mol/nm with initial coordinates as a reference.[57]
• Energy-minimized with 5000 steps of steepest descent with medium positional restraints to the heavy atoms of the peptide, 83.72 kJ/mol/nm with the initial coordinates as a reference.
• Energy-minimized with 5000 steps of steepest descent with weak positional restraints to the heavy atoms of the peptide, 4.19 kJ/mol/nm with the initial coordinates as a reference.
• Energy-minimized with 5000 steps of steepest descent without any positional restraints to the peptide.
• 100 ps of constant pressure and temperature (NPT) MD with a weakly coupled Berendsen thermostat and barostat with τ=1.0 ps, 1 fs time-step and medium positional restraints to the heavy atoms of the peptide, 83.72 kJ/mol/nm with the final energy minimized conformation as a reference. Initial velocities will be assigned using a Maxwell-Boltzmann distribution. The hydrogen atoms are restrained by the SHAKE algorithm.[58]
• 100 ps of constant pressure and temperature (NPT) MD with a weakly coupled Berendsen thermostat and barostat with τ=1.0 ps, 1 fs time-step and medium positional restraints to the heavy atoms of the peptide, 20.93 kJ/mol/nm with the final energy minimized conformation as a reference. Initial velocities should be the final velocities from step 6. The hydrogen atoms are restrained by the SHAKE algorithm.
• 100 ps of constant pressure and temperature (NPT) MD with a weakly coupled Berendsen thermostat and barostat with τ=1.0 ps, 1 fs time-step and medium positional restraints to the heavy atoms of the peptide, 4.19 kJ/mol/nm with the final energy minimized conformation as a reference. Initial velocities should be the final velocities from step 7. The hydrogen atoms are restrained by the SHAKE algorithm.
• 100 ps of constant pressure and temperature (NPT) MD with a weakly coupled Berendsen thermostat and barostat with τ=1.0 ps, 2 fs time-step and without restraints to the heavy atoms of the peptide. Initial velocities should be the final velocities from step 8. The hydrogen atoms are restrained by the SHAKE algorithm.

Production run of 1.0 μs NPT simulations performed at 310 K and 101.325 kPa pressure. The peptide and solvent with ions were separately coupled to a Parrinello-Rahman barostats and the temperatures were maintained by separate coupling to stochastic thermostats using the velocity-rescaling method of Bussi-Parrinello.[59,60] The LINCS algorithm will be used to constrain all bonds to their correct length, with a warning angle of 30°.[61,62] For analysis, a sampling frequency of 0.1 ns was utilized.

Trajectories were sampled for analysis after the system density reached a plateau.[54] An essential dynamics analysis of the trajectory was also be performed.[63,64,65] The covariance matrix for the backbone peptide atoms (-N-Cα-C-) was calculated using the covar module of GROMACS, the eigenvectors corresponding to the 10 highest eigenvalues were used to calculate the backbone configurational entropy as a function of time. This value also plateaued prior to trajectory sampling.

The time dependent ϕ,ψ dihedral angles of the central Xaa residue were extracted from the trajectories using the rama utility of GROMACS, and an in-house Python script.[66] The three-dimensional binned histograms of the normalized density were created for the ϕ,ψ dihedrals of Xaa with one degree of resolution using an R script and plotted using the MatPlotLib package in Python.[60,61,62,63,64,65,66,67] Contours representing the 98% and 99.8% of the sampled populations were generated using the smoothing algorithm of Lovell et al. and calculated using KIN2DCONT.[26,70] The ϕ,ψ dihedrals of Xaa for each peptide were subsequently grouped into their respective β, α, α_L, ε, and contiguous regions of ϕ,ψ space using density clustering methodology as implemented in DBSCAN within the Scikit Learn package for Python.[66,70,71,72,73] The minimum cluster size for DBSCAN was 1000 conformations within a 10° radius. To account for the periodicity of the ϕ,ψ dihedral angles, clustering was performed in both the -180° to +180° and 0° to 360° representations of ϕ,ψ space then deconvoluted using an inhouse PERL script to remove redundant data points from each cluster.[74] Sample proportions (ρ) and standard deviations (σ) for each population of β, α, α_L, ε, and contiguous regions of ϕ,ψ space of the Xaa residues of each peptide were calculated.[75] The sum of the probabilities for each Ac-Ala-Xaa-Ala-NH₂ system are such that;

$$1={\sum }_{i=1}^k{\rho }_i$$

(1)

there are k groups representing the β, α, α_L, ε, and contiguous regions and the total population size (n_Tot) is calculated as

$$n_{Tot}={\sum }_{i=1}^k{n}_i$$

(2)

where k and $n_i$ are defined above. The standard deviation of each ρ_ι is,

$${\sigma }_i=\sqrt{{\displaystyle \frac{{\rho }_i\left(1-{\rho }_i\right)}{n_{Tot}}}}$$

(3)

Statistical comparisons within each region of ϕ,ψ space (β, α, α_L, ε, and contiguous) were performed using two-tailed χ²-test.

$${\chi }^2={\sum }_{i=1}^k{\sum }_{j=1}^k{\displaystyle \frac{{\left(O_{ij}-E_{ij}\right)}^2}{E_{ij}}}$$

(4)

k is the number of groups being compared (the number of filled entries in a column of Table 1 minus any excluded residues), $O_{ij}$ is the observed value and $E_{ij}$ is the expected value.[75,76] An example data table is given below.

Category	Xaa	Xaai+1	…	Xaak	Totals
β, α, αL, or contig.	Oi	Oi+1	…	Oi=k	Oi,Total
Other	Oj	Oj+1	…	Oj=k	Oj,Total
					Oi,j,Total

Category	Xaa	Xaai+1	…	Xaak	Totals
β, α, αL, or contig.	Oi	Oi+1	…	Oi=k	Oi,Total
Other	Oj	Oj+1	…	Oj=k	Oj,Total
					Oi,j,Total

The expected value represents the observed value corrected to the weighted arithmetic mean of the observed values:

$$E_i=O_i{\displaystyle \frac{O_{i,Total}}{O_{i,j,Totsl}}}$$

(5)

and

$$E_j=O_j{\displaystyle \frac{O_{j,Total}}{O_{i,j,Totsl}}}$$

(6)

Results are considered statistically significant (the null hypothesis is rejected) if,

$${\chi }^2>{\chi }_{1-\alpha ,k-1}^2$$

(7)

where ${\chi }_{1-\alpha ,k-1}^2$, is the critical value with a significance level of α (α = 0.0001 for this study), and k-1 degrees of freedom.

Individual differences were deconvoluted using the Marascuillo procedure as implemented in R.[67,76,77,78,79] The methodology allows for the simultaneous comparison of all pairs of proportions when multiple populations are being tested. For a system of k populations, the number of unique pairs can be calculated from the binomial theorem.[75]

$$possiblepairs={\displaystyle \frac{k!}{2!\left(k-2\right)!}}$$

(8)

which can be shown to simplify to:

$$possiblepairs={\displaystyle \frac{k\left(k-1\right)}{2}}$$

(9)

The critical value $r_{ij}$ is calculated after choosing a significance level, ${\chi }_{1-\alpha ,k-1}^2$

$$r_{ij}=\sqrt{{\chi }_{1-\alpha ,k-1}^2}\sqrt{\left({\displaystyle \frac{p_i\left(1-p_i\right)}{n_i}}+{\displaystyle \frac{p_j\left(1-p_j\right)}{n_j}}\right)}$$

(10)

and using the proportions, $p_i$ and $p_j$, which are equal to

$$p_i={\displaystyle \frac{n_i}{{\sum }_{i=1}^k{n}_i}}$$

(11)

The absolute difference between proportions, $p_i$ and $p_j$, is compared to the critical value, $r_{ij}$, if the relationship between the two is such that,

$$\left|p_i-p_j\right|>r_{ij}$$

(12)

the null hypothesis is rejected and there is a statistically significant difference between the two pairs with a significance level of α. Results were considered statistically significant if the p-value < 0.0001.

The solvent accessible surface areas (SASA) of each β, α, α_L, ε, and contiguous region of ϕ,ψ space for the Xaa residue were calculated using the FreeSASA program with the atomic radii of Tsai et al. a water probe radius of 0.14 nm, 162000 points per sphere and the surface area algorithm of Lee and Richards.[13,41,80] Results were partitioned into whole residue, backbone, and side chain values and the mean and standard deviation for each SASA value were then calculated using the NumPy package for Python.[81] Statistical comparisons between and across Xaa residues were then performed using a Welch’s ANOVA as implemented in the Pingouin package for Python.[31,39,82] The Welch’s ANOVA is utilized due to the variation in β, α, α_L, ε and contiguous population sizes and lack of similar variances within the sampled ϕ,ψ dihedral angles for each residue caused a lack of homoscedasticity within the data sets. The f-statistic ($f_{{\nu }_1,{\nu }_2}$) for the Welch’s ANOVA is expressed as,

$$f_{{\nu }_1{,\nu }_2}={\displaystyle \frac{\frac{1}{k-1}{\sum }_{i=1}^k{w_i\left({\overline{X}}_{i\bullet }-{\overline{Y}}_{\bullet \bullet }\right)}^2}{1+\frac{2\left(k-2\right)}{k^2-1}\bullet {\sum }_{i=1}^k\left(\frac{1}{n_i-1}\right){\left(1-\frac{w_i}{w}\right)}^2}}$$

(13)

k is the number groups being analyzed, $n_i$ is the size of the population within each group, and the variances are weighted as follows,

$$w_i={\displaystyle \frac{n_i}{s_i^2}}$$

(14)

$s_i^2$ is the variance of each population i. The sum of all variances is given as:

$$w={\sum }_{i=1}^k{w}_i$$

(15)

The critical value of F is calculated with a value of significance (1 - α; α = 0.0001 in this study) and modified (${\nu }_1$) and (${\nu }_2$) degrees of freedom.

$${\nu }_1=k-1$$

(16)

$${\nu }_2={\displaystyle \frac{\left(k^2-1\right)}{3\bullet {\sum }_{i=1}^k\left(\frac{1}{n_i-1}\right){\left(1-\frac{w_i}{w}\right)}^2}}$$

(17)

The null hypothesis is rejected if the following relationship holds true.

$$f_{{\nu }_1{,\nu }_2}>F_{\left(1-\alpha \right),{\nu }_1,{\nu }_2}$$

(18)

values were considered statistically significant if the p-value <0.0001.[31,39,83,84] A pairwise difference between residues was compared using Scheffe’s method.[85] The number of possible pairwise relationships are the same as the Marascuillo procedure discussed above. The test statistic $f_s$ is calculated as,

$$f_s=\sqrt{\left(1-k\right)f_{{\nu }_1,{\nu }_2}MSE\left({\displaystyle \frac{1}{n_i}}-{\displaystyle \frac{1}{n_j}}\right)}$$

(19)

k is the number of groups within the ANOVA analysis, $f_{{\nu }_1,{\nu }_2}$ is the f-statistic calculated from the ANOVA, $n_i$ and $n_j$ are the respective sizes of the populations being compared and MSE is the mean squared error.

$$MSE={\displaystyle \frac{{\sum }_{i=1}^k\left(n_i-1\right)s_i^2}{n_{Tot}-k}}$$

(20)

The null hypothesis is rejected if the absolute difference between the means of the two pairs being compared (i and j) is greater than $f_s$.

$$\left|{\overline{X}}_i-{\overline{X}}_j\right|>f_s$$

(21)

All utilized Python packages required the use of the Pandas package for importation and data organization.[86]

6. Conclusions

The ϕ,ψ space sampling of the Xaa residues within the Ac-Ala-Xaa-Ala-NH₂ system simulated at physiologic temperatures (310 K (37° C)), and salt concentrations (150 mM NaCl) presented herein, most likely represents a realistic model of short peptides within solution. The Ac-Ala-Xaa-Ala-NH₂ peptides exist as an ensemble of conformations with the ϕ,ψ dihedrals of the central Xaa residues favoring the β region of ϕ,ψ space without the presence of secondary structure or significant hydrogen bond formations. Major conclusions of the present work are:

• Subtle differences in amino acid charge and tautomer state can affect the probability of sampling a particular region of ϕ,ψ space.

• Subtle differences in amino acid charge and tautomer state also effect the mean and standard deviation of the ϕ,ψ dihedral angles for each β, α, α_L and ε conformational region.

• Pronounced differences occur in the α_L and ε regions as a function of side chain sizes and chemical properties.

• The ε region of ϕ,ψ space is not significantly sampled for multiple amino acids while the α_L region is.

• The ε region is accessible to small hydrophobic, and hydrophilic residues (Ala, Ser, Asn, and Asp) but also bulky aromatics (His, Phe, Tyr, and Trp).

• The population density map of each amino acid’s ϕ,ψ space as shown in Figures S1 and S2 is unique and may be affected by the i-1 and i+1 adjacent residues.

• The lack of uniform sampling in the β, α, α_L, and ε regions of ϕ,ψ space in the Ace-Ala-Gly-Ala-NH₂ system reported here and in contrast to previously published QM/MM results for Ace-Gly-NH₂ indicate that the ϕ,ψ dihedral sampling of Gly is affected by the i-1 and i+1 adjacent residues.[27,28]

• Statistically significant differences exist in backbone and side chain SASA comparing β, α, α_L, ε, and contiguous regions of ϕ,ψ space.