A twist in the road less traveled: The AMBER ff15ipq-m force field for protein mimetics

The extraordinary characteristics of proteins originate in the diversity of amino acid sequences found in nature specifying a wide array of precise folded shapes. While proteins in nature are constructed primarily from a palette of 20 α-amino acids, chemists have long sought to demonstrate what is possible by expanding beyond this set, with respect to both side chain identity and backbone composition. Many oligomers with artificial amide-based backbones are capable of discrete folding behavior,¹ and these entities, termed foldamers,² have found widespread applications.³ After the inception of the field, work on folding in artificial backbones focused primarily on isolated secondary structures (e.g., helices) and, later, their assemblies (e.g., helix bundles).^4,5 More recently, attention has turned to developing artificial backbones able to recreate more sophisticated tertiary folding patterns that are typical of proteins.⁶ Advancing from secondary structure mimicry (peptidomimetics) to tertiary structure mimicry (proteomimetics) is important, as a tertiary structure is essential to function for most natural proteins. As the complexity of folds targeted for mimicry by artificial backbones has increased, so has the challenge of design and elucidating folding behavior in the resulting scaffolds. Computational tools have proven invaluable in shedding light on such issues in the context of natural proteins; however, the application of such methods to artificial backbones has been limited.

Some prior efforts have adapted computational tools to treat various non-canonical building blocks for structure prediction, enabling the design of oligooxopiperazines, oligo-peptoids, β-peptides, hydrogen bond surrogate helices, and oligosaccharides.⁷ Other work has sought to broaden the scope of backbone compositions amenable to the atomistic MD simulation in the form of CHARMM force fields for peptoids⁸ and for N- and α-methylated amino acids.⁹ Among experimental strategies for creation of artificial backbone proteomimetics with predictable and complex folding patterns, a particularly useful approach that has emerged is the blending of many artificial building block types alongside natural α residues in a single chain.¹⁰ When such “heterogeneous backbone” oligomers display a side chain sequence derived from a natural protein, the resulting agents can manifest a range of interesting folds and functions.^11–18 Biophysical analysis of folding behavior in some of these tertiary structure mimetics has revealed interesting thermodynamic effects accompanying backbone alteration.^14,19–21 MD simulation would be an invaluable tool to shed light on issues related to structure and dynamics in these and related systems; however, this requires a force field that is able to treat a range of artificial residue classes in a general manner that is internally consistent among all the residues involved, both natural and artificial.

Here, we report the development and validation of a general force field, AMBER ff15ipq-m, for the simulation of heterogeneous backbone protein mimetics. The force field is an expansion of the AMBER ff15ipq force field for canonical proteins,²² which features implicitly polarized atomic charges in the presence of explicit solvent, i.e., the best that fixed-charge force fields can offer toward modeling condensed-phase electrostatics in the absence of explicit polarization (e.g., Drude oscillators²³ and inducible multipoles²⁴). In ff15ipq-m, the original force field is expanded to include four new classes of artificial backbone units commonly used in protein mimetics (see Fig. 1): D-α-, C_α-methylated α- or aminoisobutyric acid (Aib), β³, and the cyclic β- (β^cyc) residues, trans-aminopyrrolidine carboxylic acid (APC) and trans-2-aminocyclopentane-1-carboxylic acid (ACPC). Given that the D-α-residues required no new parameterization, the expansion includes parameters for a total of 25 backbone units, consisting of 22 β³ residues (side chains corresponding to the 20 canonical amino acids including, three protonation states for His), Aib, APC, and ACPC. As shown in Fig. 1, β residues exhibit greater flexibility relative to α-residue counterparts due to the additional carbon in the backbone, which introduces a third backbone torsion θ in addition to the usual φ and ψ backbone torsions. With the aid of the GPU-accelerated AMBER molecular dynamics (MD) engine,^25–28 we have extensively validated the force field by simulating both peptides and proteins that include each backbone modification type—all on the μs timescale, yielding 30 µs of aggregate simulation time. Collectively, these results provide a useful new computational model for enabling the simulation of protein-inspired artificial backbones at the atomic level and describe a general workflow for expanding the force field further to encompass additional monomer classes in future work.

The AMBER ff15ipq force field is the latest in the lineage of Implicitly Polarized Charge (ipq) protein force fields that are based almost entirely on quantum mechanical data.^22,29 A hallmark of ipq force fields is that the atomic charges are implicitly polarized, meaning that the charges derived are intended to reproduce the mean field electron density of a molecule in explicit solvent.³⁰ For ff15ipq, the explicit SPC/E_b (extended simple point charge) water solvent model³¹ was chosen to take advantage of (i) the model’s accurate parameterization of rotational protein diffusion in solution and (ii) the relatively low cost of a three-point model compared to other four-point alternatives.

The ff15ipq force field exhibits a number of features which, when taken together, distinguish it from other contemporary force fields.²² First, the force field yields reasonable propensities for salt-bridge formation, as it was originally developed to address the overstabilization of salt bridges, a common limitation of many force fields.^22,32 In particular, newly derived atomic charges at the MP2/cc-pVTZ level of theory and atomic radii for polar hydrogens reproduce experimental probabilities of salt-bridge formation. Second, ff15ipq is far from a limited adjustment of its predecessor, ff14ipq,²⁹ but rather a complete rederivation that yields the expected balance of secondary structures for both the globular and non-globular (“disordered”) peptides/proteins. The force field not only includes new atomic charges and radii but also a greatly expanded torsion parameter set and new angle parameters. Finally, as a benefit of its parameterization with SPC/E_b explicit water, ff15ipq yields reasonable NMR observables, including J-coupling constants for the Ala₅ peptide that are comparable to those of force fields that were specifically parameterized to reproduce such experimental values (e.g., ff03w³³ and ff14SB³⁴). Moreover, a recent modification of ff15ipq methyl-side chain torsions has yielded accurate NMR relaxation involving methyl groups.³⁵ The ff15ipq force field has been implemented in dynamics engines that take advantage of GPU acceleration, e.g., AMBER²⁸ and OpenMM.³⁶ 

New parameters for each of the 25 artificial residues in the AMBER ff15ipq-m force field were derived using a workflow that was as consistent as possible with that of the parent ff15ipq protein force field.²² The same fitting classes and convergence criteria were used for deriving atomic charges and torsion parameters. The fitting classes consisted of the following: neutral, negatively charged, positively charged, glycine, and proline. The neutral class was divided into four subclasses based on the type of the first carbon of the side chain (the Cγ atom in β³ residues): alanine, two-branch residues whose Cγ atoms are attached to two heavy atoms, three-branch residues whose Cγ atoms are attached to three heavy atoms, and aromatic residues. Atomic charges for each residue were derived in the same context, i.e., within a blocked dipeptide containing acetyl (Ace) and N-methyl (NMe) caps, in the presence of explicit SPC/E_b water molecules. One minor difference from the original ff15ipq workflow is in the generation of dipeptide conformations for the derivation of both new charges and torsion parameters. In particular, the diversity of the conformations was improved by progressively restraining the backbone torsions at evenly spaced intervals and energy-minimizing the conformation rather than using high-temperature simulations.

A total of 12 new atom types were introduced for the ff15ipq-m force field (see Fig. S1). For the β³ residues, new atom types were assigned to the additional backbone methylene carbon between the carbonyl carbon and the carbon bearing the side chain. For the C_α-methylated α-residue Aib, a new atom type was introduced for the C_α atom. For the β^cyc residue APC, a new atom type was assigned to the ring nitrogen of its side chain since this form of sp³-hybridized nitrogen is not present in the side chains of other residues. For the other β^cyc residue, ACPC, it was not necessary to introduce any new atom types for its side chain since the ring carbons are sufficiently similar to those present in proline.

For each artificial residue, IPolQ atomic charges were derived by (i) generating a set of conformations for a blocked dipeptide containing the artificial residue with Ace and NMe caps; (ii) calculating the electrostatic potential for each conformation at the MP2/cc-pVTZ level of theory in both vacuum and explicit solvent; (iii) using the average electrostatic potential for all conformations to fit two sets of atomic point charges, one in vacuum and one in solvent; (iv) averaging the vacuum-phase and solvent-phase point charges to obtain an “implicitly polarized” set of charges; and (v) repeating steps (i) through (iv) until a desired convergence of the calculated charges was achieved.

To generate a set of conformations for each blocked dipeptide (20 for Aib; 30 for all other residues), φ, ψ, and θ backbone torsions were progressively restrained in vacuum between −180° and 180° using harmonic restraints with a force constant of 32 kcal/mol and subjected to energy minimization. Each conformation was then solvated in a truncated octahedral box of SPC/E_b water molecules with a 12 Å buffer and subjected to energy minimization while applying harmonic restraints to the solute (the dipeptide). To equilibrate the solvent, the solvated system was subjected to 100 ps of dynamics at constant temperature (25 °C) and pressure (1 atm) while applying harmonic restraints to the solute. A further 500 ps equilibration of the solvent was carried out during which the coordinates of surrounding solvent molecules were collected to generate a collection of point charges representing the solvent reaction field potential. This collection consists of an inner cloud of point charges based on the coordinates of the solvent molecules within 5 Å of the solute and three outer shells of point charges fitted to reproduce longer-range contributions to the solvent reaction field potentials.

Next, electrostatic potentials were calculated for each peptide conformation in two environments: a vacuum-phase environment and a solvent-phase environment, which includes the solvent reaction field potential as modeled by the collection of point charges representing the surrounding explicit solvent molecules. The calculations were carried out at the MP2/cc-pVTZ level of theory using the ORCA 4.1.2 software package.^37,38

The resulting electrostatic potentials in the vacuum and solvent phases were used to fit two corresponding sets of nuclear-centered atomic partial charges. For each artificial residue, the fitting procedure utilized the MP2/cc-pVTZ electron densities for all 30 dipeptide conformations and was carried out using the FitQ fitting procedure in the mdgx program.²⁵ Consistent with the parent ff15ipq force field, negatively charged, positively charged, and neutral residues were fit as separate groups. To account for differences in molecular charge densities due to the addition of a backbone ring, β³Pro and ACPC were fit separately from the other neutral residues, and APC was fit separately from the other positively charged residues. In addition, the C_α-methylated α-residue Aib was fit separately. Atomic charges of the Ace and NMe capping groups for each blocked dipeptide were fixed to their net neutral values in the ff15ipq force field.

In the final step of the IPolQ method of charge derivation, the vacuum-phase and solvent-phase charges were averaged to obtain an “implicitly polarized” set of point charges. The resulting set of charges were used as a starting point for further iterations of the above procedure, whereby another 30 conformations were generated using charges from the previous iteration, and an improved set of implicitly polarized atomic point charges were generated. This procedure was repeated until the charges remained within a 10% deviation of those from the previous iteration.

For each artificial residue, backbone torsion parameters were derived by (i) obtaining initial trial parameters from an existing force field for each angle and torsion involving a new atom type; (ii) generating a set of 1000 conformations for the artificial residue in the context of a blocked dipeptide with Ace and NMe caps in vacuum; (iii) calculating the QM energy of each conformation at the MP2/cc-pVTZ level of theory alongside the MM energy using the current force field’s torsion parameters; (iv) adjusting the backbone torsion parameters of the residue in order to minimize the error between the QM and MM energies; and (v) repeating steps (i)–(iv) until the error between the QM and MM energies does not change.

Initial angle and backbone torsion parameters were preferably taken from the parent ff15ipq force field if a direct analog existed. In the event that an angle or torsion term involving a new atom type did not have an exact analog in the ff15ipq force field, starting parameters were taken from the most similar parameters in the GAFF2 force field.³⁹ Angle terms were not fit in our force field and in most cases retained their fitted values from either ff15ipq or GAFF2.

An array of conformations of each blocked dipeptide were generated in vacuum by progressively restraining the φ, ψ, and θ angles at 36 evenly spaced intervals between −180° and 180° with a force constant of 32 kcal/mol and then energy minimized in the absence of restraints. In this way, the trial parameters allowed the dipeptide to extensively explore its conformational space. All conformations were inspected using φ/ψ, φ/θ, and θ/ψ plots to ensure that the conformations did not sample regions of configurational space that were too unfavorable and that the conformations did not all fall into only a few minima. Across all iterations of torsion fitting >125 000 total conformations were generated.

Next, the QM energy of each conformation was calculated at the MP2/cc-pVTZ level of theory to match the original ff15ipq protocol and subsequently fit to a function that describes the molecules MM energy as a function of torsions around the three backbone angles. The calculated QM energies were screened to exclude conformations whose energies were >10% different from the average in order to ensure a better fit to MM energies. The QM energies of all conformations were then used in a fitting procedure to generate torsion terms that gave MM energies as close to the QM energies as possible. The addition of new atom types for the residues, as previously discussed, has led to the torsion classes and subclasses by decoupling the torsion angles. Backbone torsion φ, θ, and ψ classes were separated for the following groupings of residues: neutrals, negatively charged, positively charged, glycine, proline, ACPC, APC, and Aib. The φ′ and θ′ torsion angles of the neutral residues were further separated into groups to account for differences in the first carbon in the side chain (see Table S3). The new torsion terms were then used to begin a new iteration of conformation generation. The conformation generation and fitting procedure were repeated until the root-mean-square error (RMSE) between the QM and MM energies deviated by <1% from the previous iteration.

Starting coordinates for simulations of peptides 1–4 (see Fig. 2) were built from scratch using Avogadro.^40,41 Starting coordinates for protein systems were obtained from x-ray crystal structures previously deposited in the RCSB Protein Data Bank (PDB; accession codes in Table I).^{19,20,42–45} All systems were solvated in truncated octahedral boxes of explicit SPC/E_b water molecules with a clearance of at least 16 Å for the peptides and 12 Å for the proteins. Heavy-atom coordinates for ubiquitin, GB1, and their variants were taken from the corresponding experimental structures (see Table I). Any missing atoms were built in using the SCAP program,^46–48 and net system charges were neutralized by adding Na⁺ or Cl⁻ ions, which were treated with Joung/Cheatham ion parameters intended for the SPC/E water model.⁴⁹ Hydrogens were added according to ionizable states present at pH 5.6 for GB1 and its two variants; pH 2 for peptides 1–4; and pH 7 for ubiquitin and its d-Gln35 mutant.

All simulations were carried out using the GPU-accelerated pmemd module of the AMBER 18 software package.^25–28 After preparing the system as described above and subjecting the system to energy minimization, the solvent was equilibrated in two stages, with harmonic restraints implemented on the proteins. In the first stage, a 20 ps equilibration was carried out at 25 °C at constant volume. In the second stage, a 1 ns equilibration was carried out at a constant temperature of 25 °C and a constant pressure of 1 atm. Finally, unrestrained 1 ns simulations were carried out at 25 °C and 1 atm for the peptides and proteins in Table I before data were collected for 0.5 µs and 2 µs, respectively. Consistent with NMR experiments performed with peptides 1–4 (see Sec. III F), each peptide was capped with Ace and NH₂. For each of the peptides, five 500 ns simulations were started from an extended conformation and five 500 ns simulations were started from an α-helical conformation, enforcing the total number of atoms in each simulation to be the same and yielding an aggregate simulation time of 5 µs per peptide. Both sets of five simulations result in similar distributions of φ backbone torsion angles (see Fig. S2), demonstrating convergence of the simulations. Temperatures were maintained using the Langevin thermostat with a frictional constant of 1 ps⁻¹, while pressure was maintained using the Monte Carlo barostat⁵⁰ with 100 fs between attempts to adjust the system volume. van der Waals and short-range electrostatic interactions were truncated at 10 Å, while long-range electrostatic interactions were calculated using the particle mesh Ewald method.⁵¹ A 2 fs time step was enabled by constraining bonds to hydrogen to their equilibrium values using the SHAKE algorithm.⁵² Coordinates were saved every 100 ps for analysis.

To generate Ramachandran plots of backbone conformations favored for the various artificial residues (φ, θ, and ψ torsions), umbrella sampling was carried out for a series of tetrapeptides Ace-Ala-Xaa-Ala-NMe (Xaa = each of 25 new monomers in the force field). Next, the weighted histogram analysis method (WHAM)⁵³ was applied to the umbrella sampling simulations to yield unbiased free energy profiles as a function of the torsion angles of interest for the central residue. To avoid the ring puckers of the APC and ACPC β^cyc residues becoming kinetically trapped in minor conformations, the θ torsion angles of these residues were set to 70° in the starting models of the corresponding tetrapeptides. Prior to umbrella sampling, each tetrapeptide was solvated and equilibrated as described above for unrestrained simulations with the only difference being that the duration of the second stage of the equilibration was 100 ps instead of 1 ns. For each tetrapeptide, umbrella sampling involved 1296 windows for each pair of torsions, and each window corresponded to the application of a harmonic penalty function with a force constant of 8 kcal/mol rad² to restrain the pair of torsions to angles at 10° intervals. For C_α-methylated tetrapeptides, the φ and ψ torsions of the Aib residue were restrained; for β³ residues and β^cyc residues (APC and ACPC), the φ, θ, and ψ torsions were restrained. Each window was subjected to a 0.2 ns incrementally restrained equilibration prior to a 1 ns restrained simulation at a constant temperature (25 °C) and pressure (1 atm). Coordinates of the tetrapeptide were saved every 0.05 ps. For the collective set of 25 tetrapeptides, umbrella sampling involved an aggregate simulation time of 208 µs. All umbrella sampling simulations were carried out using the GPU-accelerated pmemd module of the AMBER 18 software package.

Backbone amide J-coupling constants for peptides 1–4 were calculated from our simulations using the following Karplus relation:⁵⁴ 

³J_HN-Hα values were calculated for α residues and ³J_HN-Hβ for β³ residues to probe the φ backbone torsion angles and compare to NMR measurements (see Sec. III A).^55–57 We calibrated separate sets of Karplus coefficients for β³Ala, ACPC, and Ala based on quantum mechanical calculations using the H–N–C–H torsion that is directly measured and setting ∆ to 0°. For each residue, these calculations were based on a set of 39 dipeptide conformations containing the residue of interest. Consistent with the material used in the NMR experiments (see Sec. III A), the dipeptides were blocked with Ace and NH₂ caps (Ace-Xaa-NH₂). Conformations were generated by first restraining the φ backbone torsion to angles at 30° intervals between −180° and 180° and then at 3° intervals between −180° and −90° while subjecting the dipeptide to energy minimization with the AMBER ff15ipq-m force field. Quantum mechanical (DFT) calculations were carried out in Gaussian 09revD1⁵⁸ using the OLYP functional⁵⁹ and the pcJ-1 (triple-ζ) basis set optimized for spin–spin coupling calculations.⁶⁰ Test calculations showed that going to the larger pcJ-2 basis set changed the results by <0.2 Hz, as did changing to the hybrid functions such as PBE0 or B3LYP. The Karplus coefficients were determined by a least-squares fit to the J-coupling values computed using DFT (see Fig. S3). As noted below (see Table II and Fig. S4), results for the Ala residue in this new parameterization (which we call DFT-3) are very close to earlier fits (DFT-1 and DFT-2)⁶¹ despite differences in the density functional, basis set, and method of calculation. We have too few examples here to draw any conclusions from these small differences, especially given the experimental uncertainties in couplings of ∼0.4 Hz.

For each residue of interest, a ligand structure search was performed to determine if the monomer was present in a PDB deposited structure and, if so, the corresponding three-letter identifier. All monomers, except the β³ analogs of His and Ser, were present in at least one structure. The resulting list of identifiers (AIB, B3A, BCX, B3D, B3E, 3FB, BAL, BIL, B3K, B3L, B3M, B3X, EOE, B3Q, HMR, B3T, 1VR, HT7, B3Y, XCP, and XPC) was used to construct a query for structures containing any of these monomers as part of a polymer chain (i.e., not a free ligand). The dataset was filtered for structures that were determined by x ray and refined to a resolution of ≤2.5 Å. The resulting dataset consisted of 107 structures of backbone modified proteins. Backbone φ, θ, and ψ torsion angles were then calculated for each artificial residue, omitting the φ and ψ torsion angles for C-terminal and N-terminal residues, respectively.

Peptides 1–4 (see Table I and Fig. 2) were synthesized via the microwave-assisted Fmoc solid-phase peptide synthesis (CEM MARS microwave reactor) on Sieber amide resin (0.05 mmol scale). Protected amino acids were activated in situ by combining N-α-Fmoc-protected amino acid (4 equiv), HCTU (4 equiv), and DIEA (6 equiv) in N-methyl-2-pyrrolidone (final concentration 0.1M amino acid). Couplings were performed at 70 °C for 4 min. Deprotections of the Fmoc group were performed at 80 °C for 2 min using 20% v/v 4-methylpiperidine in DMF. Each peptide was capped at the N-terminus with an acetyl group. DMF/DIEA/acetic anhydride (8:2:1 by volume; 2 ml) was added to the resin and allowed to stir at room temperature for 20 min. Post synthesis, the resin was successively washed with DMF, DCM, and methanol and allowed to dry in a vacuum desiccator for 30 min. Cleavage of the peptide from the resin was performed by treatment of the resin with 1% v/v TFA in DCM for 20 min–30 min with gentle stirring. The cleavage flowthrough was collected in a scintillation vial, and excess DCM evaporated under nitrogen. For peptide 4, the resin was washed an additional five times with cleavage solution after incubation to maximize the cleavage yield. Cleaved peptides were purified by preparative HPLC on a Phenomenex Luna Prep C18 column (10 µM particle size and 100 Å pore size). The identity of each peptide was confirmed by electrospray-ionization mass spectrometry (ESI-MS; see Table S1), and the purity of each peptide was determined by analytical HPLC (see Fig. S5) and NMR (see Figs. S7–S10). Purified peptides were lyophilized prior to NMR analysis.

Lyophilized peptides were weighed and dissolved in 9:1 H₂O:D₂O with 0.2 mM DSS and the pH adjusted to 2. NMR experiments were performed on a Bruker Avance 700 MHz spectrometer at 300 K. Spectra acquired for each peptide consisted of 1D ¹H (128 scans), 2D NOESY (200 ms mixing time and 16 scans), and 2D COSY (8 scans). All spectra were processed using TOPSPIN, employing DSS as an internal standard for chemical shift and concentration. Final concentrations were 10 mM, 9 mM, 17 mM, and 20 mM for peptides 1–4, respectively. J-coupling values were measured on TOPSPIN through deconvolution and peak fitting of amide H_N doublets (see Fig. S6) using a mixture of Lorentzian and Gaussian curves; reported experimental uncertainties (0.4 Hz) were estimated as half the spectral resolution. Variable-concentration 1D ¹H experiments for peptide 1 in the range of 1 mM–10 mM showed that spectra were concentration-independent in this range, demonstrating that the peptide is monomeric under the conditions of the measurement (see Fig. S11).

To assess the accuracy of the newly derived torsion and angle parameters for each of the 25 artificial amino acid residues, we monitored the RMSE of the AMBER ff15ipq-m molecular mechanical (MM) potential energy surface from the target quantum mechanical (QM) potential energy surface. As shown in Fig. 3, the RMSE values range from 1.3 kcal/mol to 2.3 kcal/mol. These values are slightly higher, by 0.1 kcal/mol–1 kcal/mol, relative to those of their canonical counterparts in the parent AMBER ff15ipq force field²² due, in part, to the increased backbone flexibility that results from the additional methylene group in the β³ residues and a greater likelihood for forming backbone hydrogen bonds.

To characterize the backbone conformational preferences of each artificial residue modeled by the ff15ipq-m force field, we carried out umbrella sampling of Ace-Ala-Xaa-Ala-NMe tetrapeptides with the artificial residue at the central position, progressively restraining the backbone torsion angles. Given that the simultaneous variation of all three backbone torsion angles in the β³ and β^cyc residues is computationally prohibitive within the same umbrella sampling simulation, only two of the backbone torsion angles were restrained at a time, yielding a total of three free energy surfaces (i.e., φ/ψ, φ/θ, and ψ/θ Ramachandran plots) after the application of the weighted histogram analysis method (WHAM; see Sec. III). To validate the calculated Ramachandran plots from the simulations, we compared predicted torsional preferences to experimentally observed conformations from 107 x-ray structures of proteins containing one or more of the newly parameterized residues from the PDB (see Data S1).

Overall, the backbone torsional preferences observed for each artificial residue in our tetrapeptide simulations are in qualitative agreement with the behavior of the corresponding residue seen in the PDB dataset (see Fig. 4 and Fig. S12). Of note, the experimental dataset is dominated by structures where the artificial residue is in an α-helical context as α-helix mimicry has historically been a major focus in development of protein mimetics. As a result, the corresponding regions in the β-residue Ramachandran plots (at around φ/ψ −110°/−110°, ψ/θ −110°/90°, and ψ/θ −110°/90°) are heavily populated by experimental data points.

Consistent with the principles that informed earlier application to control β-residue conformational preferences in homogeneous β-peptide backbones⁶² and heterogeneous backbones displaying biological side chain sequences,⁶³ the Ramachandran plot for the β^cyc residue ACPC is highly restricted and favors a right-handed helical conformation. Ramachandran plots for β³ residues bearing proteinogenic side chains show more accessible minima, many outside the helical fold seen in experimental structures. As expected, βGly exhibits the greatest diversity of conformations and falls within the same free energy well that is populated by the PDB dataset. Conformations of representative γ-branched (β³Ile) and charged (β³Lys) residues also fall within the free energy wells populated by the PDB dataset. Taken together, the above findings validate the conformational preferences of the newly parameterized artificial residues modeled by ff15ipq-m. Furthermore, given that the PDB dataset consists of structures for entire proteins containing the artificial residues, the qualitative agreement of our simulation-based Ramachandran plots with those from the PDB dataset indicates that the conformational preferences observed by simulations in the context of tetrapeptides are predictive of those in the context of proteins.

To complement our qualitative comparisons for the individual artificial residues, we quantitatively assessed the accuracy of backbone conformational predictions by ff15ipq-m for selected artificial residues by calculating backbone amide J-coupling constants for tetrapeptides 1–4 (see Fig. 2 and Table I) and comparing to the observations in NMR experiments on the same compounds. J-coupling constants corresponding to the φ torsion angle (³J_HN-Hα for α residues and ³J_HN-Hβ for β residues) were calculated using the Karplus relation and three sets of DFT-based Karplus coefficients: DFT-1, DFT-2, and DFT-3 (see Table S2). The DFT-1 and DFT-2 coefficients were previously derived for canonical residues,⁶¹ whereas the DFT-3 coefficients were derived in this study for the β³Ala and ACPC residues and re-calibrated for the Ala residue.

As shown in Table II, when using the newly derived DFT-3 coefficients, the ff15ipq-m force field yields J-coupling constants that are generally within error of experiment. The only exception is that the J-coupling constant for β³Ala is just outside of the error for the experimental value, indicating a conformational ensemble with a slightly less extended character in the simulation compared to the experiment. The β³Ala and ACPC β^cyc residues appear to have significantly different J-coupling constant distributions than those of the Ala and d-Ala residues (see Fig. S4), likely, due to the preference toward an extended conformation and the increased flexibility of the peptide backbone. Relative to the DFT-1 and DFT-2 coefficients, the DFT-3 coefficients yield improved agreement with experimental J-coupling values for the β³Ala and ACPC residues, while preserving the expected trends in the ensemble-averaged values and distributions (see Fig. S4) of J-coupling constants for the Ala, d-Ala, β³Ala, and ACPC residues.

To assess the kinetic stability of protein mimetics with the ff15ipq-m force field, we carried out μs-timescale simulations for a set of three backbone modified proteins that have high-resolution crystal structures available and together contain the various types of artificial backbone units that have been parameterized for the force field (i.e., d-α-, Aib, β³, and β^cyc residues). The set of protein mimetics consists of (i) a ubiquitin variant that contains d-Gln35 in place of l-Gln35; (ii) a GB1 variant that contains Aib24, β³Lys28, β³Lys31, and Aib35; and (iii) a GB1 variant that contains ACPC24, β³Lys28, β³Lys31, and ACPC35. For comparison, simulations were also performed for the corresponding wild-type proteins. All simulations were carried out at 25 °C. For further simulation details, see Table I and Sec. III.

As shown in Fig. 5, the overall structures of both the wild-type and variant proteins remained stable over their entire simulations. Only small deviations from the corresponding crystal structures were observed, and these were similar in magnitude to deviations observed in simulations of the corresponding wild-type proteins. The average backbone root-mean-square deviation (RMSD) from the corresponding crystal structure is <2.3 Å for all systems examined. Notably, the secondary structures of both wild-type ubiquitin and GB1 are maintained throughout the simulations in their corresponding variants (see Fig. S13).

In summary, we have reported here the development and validation of a general force field, AMBER ff15ipq-m, for the molecular modeling of protein mimetics with backbones bearing combinations of canonical α residues and four classes of artificial backbone units: d-α-, Aib, β³, and β^cyc residues (APC and ACPC). The ff15ipq-m force field is an expansion of the AMBER ff15ipq force field for natural proteins and includes 472 unique atomic charges and 148 unique backbone torsion terms. As was the case for creating the ff15ipq force field, the derivation of new IPolQ parameters was greatly facilitated by the mdgx program of AMBERTools via a “sweeping optimization” approach.

The AMBER ff15ipq-m force field features all the following for the modeling of the newly parameterized artificial backbone units: (i) a charge set that accounts directly for water-induced polarization, as consistent with the AMBER IPolQ lineage of force field; (ii) conformational preferences of individual residues that are consistent with those found in crystal structures of protein mimetics; (iii) strong agreement with NMR J-coupling constants for short peptides using newly derived Karplus coefficients; and (iv) maintenance of globular protein folds with altered backbones on the μs timescale.

The new force field reported here and the path to its development have a number of potential implications with respect to understanding folding behavior in protein-inspired artificial backbones. The ff15ipq-m force field can be used to conduct atomistic simulations on an array of experimentally characterized systems, shedding light on fundamental issues related to the impact of altered backbone composition on folding and function. Furthermore, simulations on simple de novo systems may provide insights into regions of backbone configurational space that are accessible to heterogeneous backbones but not previously observed experimentally due to the historical focus on mimicry of α-helical structures. Finally, the present results demonstrate the applicability of the IPolQ methodology to backbones beyond α peptide. The workflow described should prove readily applicable to additional residue classes, raising the possibility of future expansion of ff15ipq-m to encompass a broad array of macromolecules in chemical space that are adjacent to but distinct from proteins.

See the supplementary material for new atom types for artificial backbone units; probability distributions of torsion angles and J-coupling constants for central residues of peptides 1–4; DFT-3 Karplus coefficients and least-squares fits for extracting these coefficients; HPLC traces, ESI-MS data, and NMR spectra for peptides 1–4; Ramachandran plots for all 25 artificial amino acid residues in the context of tetrapeptides; secondary structure content of ubiquitin and GB1 variants in simulations on the μs timescale; backbone torsion classes of the AMBER ff15ipq-m force field; PDB codes of the 107 crystal structures used for qualitative validation of computed Ramachandran plots; and ff15ipq-m force field files for use with the AMBER software package.

A.T.B., H.E.P., and J.M.G.L. contributed equally to this work .

The data that support the findings of this study are available within the article and its supplementary material. All AMBER ff15ipq-m force field parameters will be included as an update to the AmberTools20 distribution and are provided in the supplementary material.

This work was supported by the University of Pittsburgh Arts and Sciences Graduate Fellowship to A.T.B., the NSF (Grant No. CHE-1807301) to L.T.C. and W.S.H., the NIH (Grant No. 1R01GM115805-05) to L.T.C., and the NIH (Grant No. U54AI50470) private clusters in the laboratory of D.A.C. Computational resources were provided by the University of Pittsburgh’s Center for Research Computing and private clusters in the CASE Laboratory.