# Robust Conformational Space Exploration of Cyclic Peptides by Combining Different MD Protocols and Force Fields

Samuel Murail, Jaysen Sawmynaden, Akli Zemirli, Maud Jusot, Fabio Pietrucci, Jacques Chomilier, Pierre Tufféry, Dirk Stratmann

<figure><img src="image_1.png" /></figure>

## To cite this version:

Samuel Murail, Jaysen Sawmynaden, Akli Zemirli, Maud Jusot, Fabio Pietrucci, et al. Robust Conformational Space Exploration of Cyclic Peptides by Combining Different MD Protocols and Force Fields. Journal of Chemical Theory and Computation, 2025, 21 (19), pp.10018-10034. 10.1021/acs.jctc.5c01123. [hal-05415840](https://hal.science/hal-05415840v1)

HAL Id: hal-05415840

https://hal.science/hal-05415840v1

Submitted on 15 Dec 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
# Robust Conformational Space Exploration of Cyclic Peptides by Combining Different MD Protocols and Force Fields

Samuel Murail, Jaysen Sawmynaden, Akli Zemirli, Maud Jusot, Fabio Pietrucci, Jacques Chomilier, Pierre Tufféry, and Dirk Stratmann*

<figure><img src="image_3.png" /></figure>

Cite This: J. Chem. Theory Comput. 2025, 21, 10018–10034

<figure><img src="image_4.png" /></figure>

Read Online

## ACCESS

<figure><img src="image_5.png" /></figure>

Metrics & More

<figure><img src="image_6.png" /></figure>

Article Recommendations

<figure><img src="image_7.png" /></figure>

Supporting Information

**ABSTRACT:** Cyclic peptides are an important class of pharmaceutical drugs. We used replica-exchange molecular dynamics (REMD) and simulated tempering (ST) simulations to explore the conformational landscape of a set of nine cyclic peptides. The N-ter to C-ter backbone-cyclized peptides of 7-10 residues were previously designed for high conformational stability with a mixture of L- and D-amino acids. Their experimental NMR structures are available in the protein data bank (PDB). For each peptide, we tested several force fields, namely, Amber96, Amber14, RSFF2C, and Charmm36m in implicit and explicit solvents. We find that the variability of the free energy maps obtained from several protocols is larger than the variability obtained by just repeating the same protocol. Running multiple protocols is therefore important for the convergence assessment of REMD or ST simulations. The majority of the free energy maps showed clusters with a high RMSD compared to the native structures, revealing the residual flexibility of this set of cyclic peptides. The high RMSD clusters had in some cases the lowest free energy, rendering the prediction of the native structure more difficult with a single protocol. Fortunately, the combination of four implicit solvent REMD and ST simulations, mixing the Amber96 and Amber14 force fields, predicted robustly the native structure. As implicit solvent simulations in the REMD or ST setup are up to one hundred times faster than explicit solvent simulations, running four implicit solvent simulations is a good practical choice. We checked that the use of an explicit solvent REMD or ST simulation, taken alone or combined with implicit solvent simulations, did not significantly improve our results. It results in our combination of four implicit solvent simulations being tied in terms of success rate with much more expensive combinations that include explicit solvent simulations. This may be used as a guideline for further studies of cyclic peptide conformations.

<figure><img src="image_8.png" /></figure>

## INTRODUCTION

Peptides have attracted increasing attention over the past decade as a viable alternative to small molecules for developing drug molecules capable of interfering with protein-protein interactions particularly well.¹⁻⁵ While linear peptides composed of only natural amino acids suffer from degradation, cyclic peptides display a stronger resistance to proteases.⁶ In addition, their conformational space is usually more restricted than that of linear peptides,⁷ reducing the entropy loss upon binding and leading to a stronger binding affinity.⁸,⁹ Cyclization is usually achieved by a peptide bond between the N-ter and C-ter amino acids or by disulfide bonds between two cysteine residues. Combining both approaches allows the design of hyperstable cyclic peptides.¹⁰ Conformational stability can also be obtained by a combination of L- and D-amino acids,¹¹ without the need for disulfide bonds at least for small cyclic peptides up to about ten residues.¹² Combining the cyclization with chemical modifications, like noncanonical amino acids¹³,¹⁴ or N-methylation, can increase the specificity to the target¹⁵ and reduce the conformation flexibility.¹⁶ Cyclic

peptides that incorporate chemical modifications are named macrocycles, and they represent a class of drugs on their own right,¹⁷ with an important field of research. The chemical space available to macrocycles is huge, and as for small molecules, virtual libraries of billions of compounds have been developed.¹⁸ To interfere with intracellular targets, drugs must pass through the cell membrane. Fortunately, cyclic peptides can pass this barrier, as documented in the recent database CycPeptMPDB gathering experimental membrane permeability measurements of over 7000 cyclic peptides, of which more than the half present high membrane permeability,¹⁹ and in vitro workflows have been developed recently for generating orally available cyclic peptides.²⁰ Overall, there

Received: July 9, 2025

Revised: September 1, 2025

Accepted: September 8, 2025

Published: September 26, 2025

<figure><img src="image_9.png" /></figure>
is a bright future for cyclic peptides and macrocycles in therapeutic drug development²¹ with an average of about one cyclic peptide drug approved per year.⁵ The numerous possibilities for in vitro and in silico design approaches for peptide therapeutics in general are summarized in several recent reviews.²²⁻²⁸

While cyclic peptides are less flexible than linear peptides, some degree of flexibility remains depending on the length and composition of the amino acid sequence, as well as the use of chemical modifications.¹⁶,²⁹ In a drug design context, it is important to be able to access the entire conformational landscape accessible to a cyclic peptide for two main reasons: First, a poorly flexible cyclic peptide will likely lose less conformational entropy upon binding than a more flexible peptide.⁸,⁹ To quantify this, the accessible conformational space can be a possible indicator. Second, a flexible cyclic peptide usually has a different conformation in its free form than in its bound form inside the protein–peptide complex, due to additional protein–peptide interactions.³⁰ For an optimal binding affinity, it is important that the bound form is already explored inside the conformational ensemble of the free form (conformational selection) or that the free form can be converted into the bound conformation (induced fit) without having to cross a too high free energy barrier.³¹ For some applications, high flexibility is needed. For example, it is sometimes necessary to travel “a tunnel” to the binding site. So, cyclic peptides cannot always be designed as “rocks” but have to incorporate an appropriate level of flexibility. Experimental structure determination methods like NMR or X-ray crystallography will resolve the lowest free energy structure and perhaps a second minor conformation but not the whole conformational space and the free energy barriers among several metastable conformations. Therefore, conformational space explorations are usually done with in silico methods.

The conformational space explored by cyclic peptides can be characterized in silico with molecular dynamics (MD) (reviewed in ref 32), Monte Carlo simulations, and other sampling methods. The Rosetta software suite³³ uses Monte Carlo simulations and kinematic closure algorithms to sample cyclic peptides conformations.¹⁰ A similar kinematic closure algorithm has been implemented in EGSCyP (Exhaustive Grid Search for Cyclic Peptides)³⁴ that uses an exhaustive exploration for cyclic pentapeptides composed of standard, D-amino acid, and N-methylated amino acids. Conventional molecular dynamics simulations (cMD) are, in general, not sufficient to pass the high free energy barriers that separate different conformations of cyclic peptides. Therefore, enhanced sampling molecular dynamics simulations are used, like replica-exchange MD³⁴⁻⁴² and bias-exchanged metadynamics (BE-MetaD).⁴⁰,⁴¹,⁴³⁻⁵² Machine learning can even accelerate further the exploration of the conformational space of cyclic peptides.⁵³,⁵⁴ AlphaFold2⁵⁵ has also been applied to cyclic peptides,⁵⁶⁻⁵⁸ but it is limited by its design to the prediction of the 3D structure and not to the identification of a conformational ensemble. Whereas AlphaFold2 is limited to natural amino acids, the recent AlphaFold3⁵⁹ may be used to predict the structure of cyclic peptides that incorporate modified amino acids, as it can be applied to any small molecule in general, but this remains the subject for further investigations.

In the present study, we assess the performance and convergence of several force fields in molecular dynamics (MD) simulations for the exploration of the conformational

space of a set of nine cyclic peptides, ranging from seven to ten amino acids. We run in parallel two enhanced sampling methods, replica-exchange molecular dynamics (REMD) and simulated tempering (ST), to exclude bias from either sampling method.

A very recent study⁶⁰ of the Yu-Shan Lin group tested several force fields on a set of 12 smaller cyclic peptides, ranging from five to seven amino acids. The smaller size stabilized the eight sequences without D-amino acids. The four remaining pentapeptides had one or two D-amino acids. They performed BE-MetaD simulations in an explicit solvent for most of the force fields tested. In our study, we focused on implicit solvent simulations to assess if this kind of solvent model is sufficient for the elucidation of the conformational space of free-form cyclic peptides. We also focused on medium-sized cyclic peptides, as the conformational space can be more difficult to be resolved. Therefore, our study is complementary to the similar study⁶⁰ of the Yu-Shan Lin group.

## MATERIALS AND METHODS

**Cyclic Peptide Data Set.** Most cyclic peptides present in the Protein Data Bank (PDB) are bound with other molecules or have chemical modification (like N-methyl). In the present study, we choose a set of nine cyclic peptides available in free form in the PDB and having a medium size from seven to ten amino acids. All but one are peptides with a mixture of L- and D-amino acids, optimized for conformational stability by the David Baker group using Rosetta.¹² We employ here the same peptide names as in their paper, like “7.B” for a peptide that has seven amino acids and that is the second design of this size (Table 1). Peptide 8.C in Table 1 is not present in the study of

<table><caption>Table 1. Cyclic Peptide Data Set<sup>a</sup></caption><thead><tr><th>PDB code</th><th>peptide</th><th>sequence</th><th>atoms</th><th>NMR ensemble RMSD [Å]</th></tr></thead><tbody><tr><td>6BE9</td><td>7.A</td><td>TkNDTnp</td><td>104</td><td>0.9 ± 0.38</td></tr><tr><td>6BEW</td><td>7.B</td><td>hPdqseP</td><td>100</td><td>0.74 ± 0.22</td></tr><tr><td>6BF3</td><td>7.C</td><td>QDPpKtd</td><td>105</td><td>0.45 ± 0.087</td></tr><tr><td>6BE7</td><td>8.A</td><td>DDPTprQq</td><td>124</td><td>0.56 ± 0.15</td></tr><tr><td>6BEN</td><td>8.B</td><td>rQpqRePQ</td><td>142</td><td>0.44 ± 0.077</td></tr><tr><td>6WPV</td><td>8.C</td><td>GTVAVQFL</td><td>119</td><td>0.23 ± 0.056</td></tr><tr><td>6BEO</td><td>9.A</td><td>pPYhPKDLq</td><td>150</td><td>0.51 ± 0.2</td></tr><tr><td>6BEQ</td><td>10.A</td><td>AARvpRltPE</td><td>160</td><td>0.65 ± 0.18</td></tr><tr><td>6BER</td><td>10.B</td><td>EvDPehpNap</td><td>141</td><td>0.8 ± 0.18</td></tr></tbody></table>

<sup>a</sup>All structures of the peptides in free form (i.e., not bound to a protein) have been solved by NMR. Upper case letters are L-amino acids and lower case letters are D-amino acids. As a reference structure, we used the best representative conformer of the NMR ensemble which is here the first model for all nine PDB files. The conformational variability among the NMR ensemble is quantified with the average RMSD of the 19 models compared to the first model.

Hosseinzadeh et al. It is composed of standard amino acids only.⁶¹ We added it to see if a cyclic peptide of this size can be as stable as the peptides that are stabilized by D-amino acids and prolines.

**Implicit Solvent REMD Protocol.** Our REMD protocol is similar to the one we used in a former study on cyclic pentapeptides³⁴ and inspired from the protocol presented in the work of Wakefield et al.³⁹

The starting structures of the cyclic peptides were constructed manually with ChimeraX.⁶² Topology files were
generated with Ambertools⁶³ and converted in GROMACS format with acyppe.⁶⁴

Temperature REMD simulations were set up with eight replica with the OBC (Onufriev, Bashford, and Case)⁶⁵ GBSA implicit solvent and run with GROMACS 5.1.2;⁶⁶ more recent GROMACS versions do not include implicit solvents any more. A short minimization was first done with an alternation of one step of steepest descent every 500 conjugate gradient steps. The maximum number of steps was set to 50,000 with a step size of 0.01 nm and an energy convergence criterion of 10 kJ/(mol·nm). The thermalizations and simulations were done in a temperature range from 300 to 455 K. The temperature values of the eight replica were chosen to keep the probabilities of accepting exchanges the same (300, 318, 337.97, 358.81, 380.85, 404.27, 429.12, and 455.50 K). The acceptance ratio was high, as 55 to 65% of all exchange attempts were successful, thanks to the small number of atoms (100 to 150). Only a short phase of thermalization in NVT was realized during 100 ps (50,000 steps, 2 fs time step). A modified Berendsen thermostat with an additional stochastic term that ensures a correct kinetic energy distribution⁶⁷ was employed with the time constant $τ = 0.1$ ps.

A good sampling in REMD depends on the time step between exchange attempts between neighboring replicas.⁶⁸ No consensus choices exist for this parameter. The time step used in the literature is between 0.01 and 100 ps,³⁵,⁶⁹ and it depends on the system and the process studied. Here, we chose 1 ps as the time step between exchange attempts. For each force field and peptide, we produced trajectories of 1 µs with a 2 fs time step for each replica, i.e., 8 µs for each REMD simulation. The atomic coordinates were recorded every 10 ps. We repeated each REMD simulation five times (runs 1 to 5) to access convergence. In total, we produced in this study with the REMD protocol 8 µs × 5 runs × 9 peptides × 3 force fields = 1080 µs of cumulative trajectory time (see Table 2).

<table><caption>Table 2. MD Protocols<sup>a</sup></caption><thead><tr><th>protocol</th><th>force field</th><th>solvent</th><th>replica</th><th>runs</th><th>cumulated<br>trajectory time [µs]</th></tr></thead><tbody><tr><td>REMD</td><td>Amber96</td><td>implicit</td><td>8</td><td>5</td><td>360</td></tr><tr><td>REMD</td><td>Amber14</td><td>implicit</td><td>8</td><td>5</td><td>360</td></tr><tr><td>REMD</td><td>RSFF2C</td><td>implicit</td><td>8</td><td>5</td><td>360</td></tr><tr><td>REMD</td><td>RSFF2C</td><td>explicit</td><td>32</td><td>1</td><td>288</td></tr><tr><td>REMD</td><td>Charmm36m</td><td>implicit</td><td>8</td><td>2</td><td>144</td></tr><tr><td>ST</td><td>Amber96</td><td>implicit</td><td>1</td><td>2</td><td>180</td></tr><tr><td>ST</td><td>Amber14</td><td>implicit</td><td>1</td><td>1</td><td>90</td></tr><tr><td>ST</td><td>Amber14</td><td>explicit</td><td>1</td><td>1</td><td>90</td></tr><tr><td>ST</td><td>Charmm36m</td><td>implicit</td><td>1</td><td>1</td><td>90</td></tr></tbody></table>

<sup>a</sup>Each MD protocol was run on the nine peptides from Table 1. The trajectory length for REMD is 1 µs and for ST 10 µs. The cumulated trajectory time is the product of this trajectory length with the number of replica and with the number of repeated runs. The total cumulated trajectory times of the implicit and explicit solvent runs are 1584 µs and 378 µs, respectively.

To produce 1 µs per replica, the cost is between 1000 and 2000 CPU hours per REMD simulation. So, in total, we consumed about 5 runs × 9 peptides × 3 force fields × 1000–2000 CPU hours = 135–270 k CPU hours, which is quite reasonable thanks to the use of an implicit solvent.

Four force fields were tested with the REMD protocol: Charmm36m,⁷⁰ Amber96,⁷¹ Amber14SB,⁷² and the coil-

library-based Residue Specific Force Field RSFF2C⁷³ which is derived from the Amber14SB force field and has shown better performance for peptides,⁶⁰ proteins, and loops.⁷⁴ We included the older Amber96 force field, as we employed it in our former study³⁴ and as it still showed a good performance on cyclic peptides.

**Explicit Solvent REMD Protocol.** We also set up an explicit solvent REMD protocol with the RSFF2C force field, which has been tested only in explicit solvents in previous studies. Using the standard TIP3P explicit solvent model, a cubic water box with a minimum distance between the peptide and the box walls of 1.0 nm was equilibrated with GROMACS 2018.8.⁶⁶ In the case that the peptide does not have a zero net electric charge, the total charge of the peptide + water system has been neutralized by adding sodium or chloride ions to the solvent. The energy of the whole system was then minimized with a steepest descent algorithm, and each replica was equilibrated to its temperature in two steps, first with a simulation of 100 ps in the NVT ensemble followed by a simulation of 200 ps in the NPT ensemble. The time step of the equilibration and production simulations was set to 2 fs. A nonbonded cutoff of 1.0 nm was used with the Particle Mesh Ewald (PME) algorithm for long-range electrostatic interactions. The same Berendsen thermostat as for the implicit solvent REMD protocol was employed with the time constant $τ = 0.1$ ps and the Parrinello-Rahman pressure coupling⁷⁵,⁷⁶ in the NPT ensemble with the time constant $τ = 2.0$ ps.

We employed 32 replicas, as this value fitted well with our computer architecture. The temperature values were generated with the web-server⁷⁷ using the algorithm published in ref 78 using the same temperature range as for the implicit solvent protocol (see Supporting Information for the temperature list).

We obtained an acceptance ratio of 25 to 35% for replica exchanges, depending on the number of atoms, which varied from 4300 to 6900 atoms (solvent + peptide) as a function of the size of the peptide. Usually an acceptance ratio of around 30% is targeted, and we were pleased to see that the final acceptance ratios were about twice the values (12–15%) predicted by the web-server mentioned above.

We obtained 50 to 75 ns/day on two recent 32 core AMD EPYC 7452@2.35 GHz CPUs. Running all nine peptides on 18 CPUs with 32 cores in parallel, we obtained after about 13 to 20 days a single 1000 ns trajectory per peptide. In total, this single run consumed about 200,000 CPU hours (15 days × 64 cores × 9 peptides). Repeating the simulation with five runs like for the implicit solvent was out of reach for this study, neither doing the whole benchmark of the three force fields in explicit solvent with REMD, as taken together it would have consumed 15 times more CPU hours, i.e., about 3 million CPU hours.

**ST Simulations.** The peptide structures were prepared using the pdbfixer module of the OpenMM package, where a bond was specifically added in the peptide topology between the N backbone atom of the first residue and the C backbone atom of the last residue.⁷⁹ For the simulation run with an explicit solvent, peptides were solvated in truncated octahedron boxes with a padding of 1.5 nm and the TIP3P water model. Counter ions were added to counter the charge of the peptide, together with a 150 mM concentration of Na⁺Cl⁻ ions. Simulations were conducted with OpenMM 7.7 Molecular Dynamics software.⁷⁹

Energy minimization was run for up to 10,000 steps, followed by an equilibration step of 10 ns in the NPT ensemble.
<figure><img src="image_10.png" /><figcaption>**Figure 1.** (a) REMD is a random walk, and conformational space exploration is not the same with short (<100 ns) and long (>500 ns) simulations. The determination of the simulated time necessary to get a converged simulation is important to conclude which clusters are most populated. In this example, the population of the clusters converged after 500 ns. The RMSD values of the shown profiles are calculated by comparing the frames to the experimental reference structure. (b) Comparison of the RMSD profiles of five REMD runs of 1 µs (first 10% skipped) to the mean RMSD profile which combines all five profiles. (c) The Jensen-Shannon divergence (JSD) is used as a metric to quantify the difference between two profiles. If the JSD value is above 0.05, we consider that the RMSD profile is different from the mean RMSD profile. Convergence of all five REMD runs is reached if all five curves are below 0.05; here, only three of five runs converge, but the other two runs (runs 3 and 4 here) are not far from the chosen threshold.</figcaption></figure>

ensemble. During production and equilibration, a 4 fs integration time step was used with the help of constraining all bonds involving a hydrogen atom, as well as using heavy hydrogens assigned to 3 atomic mass units. A nonbonded cutoff of 1.0 nm was used with the PME algorithm for long-range electrostatic interactions. A pressure of 1.0 bar was applied during equilibration and production using the Monte Carlo algorithm<sup>80,81</sup> and a 25 time steps interval. Temperature was equilibrated using the Langevin Middle integrator<sup>82</sup> and 1 ps<sup>-1</sup> friction coefficient.

During the equilibration, the temperature was fixed to 300.0 K. ST simulations were conducted using a protocol inspired from a previous work<sup>83</sup> using a python script written by Peter Eastman and modified to implement the weight calculation of Park and Pande<sup>84</sup> and the on-the-fly weight calculation of Nguyen et al.<sup>85</sup> During ST simulations, temperature exchanges were attempted every 10 ps, and depending on the size of the peptide, 15 to 20 temperature ladders spaced exponentially between 300.0 and 500.0 K were used. The ST simulation production run had a length of 10 µs.

Simulations with an implicit solvent were run in the NVT ensemble with no cutoff for nonbonded interaction. We used only 8 temperature ladders, and the temperature exchange was attempted every 6 ps. Implicit solvent models were GBN2<sup>86</sup> for the simulation run with the Amber14 force field, as we used the OBC model<sup>87</sup> with Amber96 and Charmm36m force fields.

**Analysis of REMD and ST Simulations.** All analyses in this article were done on the REMD trajectory at 300 K and the parts of the ST trajectory that were at 300 K. We skipped the first 10% of the trajectory to reduce the bias of the starting conformation. We compared each frame to the experimental reference structure, which is here the best representative conformer of the NMR ensemble (first model here). Two measures were calculated for each frame: the backbone RMSD to the experimental reference structure and the radius of gyration. Both were calculated with the mdtraj<sup>88</sup> python package on the heavy atoms of the backbone (N, C<sub>α</sub>, C, and O). Free energy maps as a function of both metrics were calculated and plotted with python using a custom function adapted from the PyEMMA package.<sup>89</sup> Therefore, a 2D histogram is generated with numpy (100 bins for each axis) and translated into free energy values according to the formula

$F = -kT \times \log(\rho)$, with $\rho$ being the density at each point of the free energy map.

**Convergence Analysis.** To determine the convergence of an MD simulation is not a trivial problem. As MD is a random walk, short simulations cannot sample enough conformational landscape and explore all free energy minima. Moreover, the populations of the clusters can change during longer trajectories, as shown in Figure 1a.

To quantitatively assess the convergence of our implicit solvent REMD simulations, we performed five independent runs for each peptide and force field combination. If the RMSD density profiles of all five runs superpose within a certain error margin, then we consider that the REMD simulations converged. The RMSD values are calculated by comparing the frames to the experimental reference structure, but in principle, any structure could be used as a reference here. To quantify this superposition, we first build a mean density profile combining all five runs (see Figure 1b). To quantify how much each RMSD profile differs from the mean RMSD profile, the Jensen-Shannon divergence (JSD<sup>90</sup>) is used:

$$ \text{JSD}(P||Q) = \frac{1}{2} D_{\text{KL}}(P||M) + \frac{1}{2} D_{\text{KL}}(Q||M) \quad (1) $$

where P is the RMSD density profile, Q is the mean RMSD density profile, M = 1/2(P + Q), and D<sub>KL</sub> is the Kullback-Leibler (KL) divergence.

Unlike KL, JSD is a symmetric metric (JSD(A, B) = JSD(B, A)), and the range of its values is between 0 (similar profiles) and 1 (completely different profiles). If the JSD value is above 0.05, we consider that the RMSD profile is different from the mean RMSD profile (see Figure 1b as an example of JSD values around 0.05).

For each REMD run, 100 RMSD density profiles are calculated along the trajectory (0 - t<sub>1</sub>, 0 - t<sub>2</sub>, 0 - t<sub>3</sub>, ..., 0 - t<sub>100</sub>) to capture the simulation time evolution of the convergence, which allowed us to plot the JSD curves shown in Figure 1c.

**2D Projections of the Conformational Landscape.** The conformation landscape is projected here to a 2D map, using two complementary projections: the first one has the radius of gyration of the cyclic peptide backbone on the x-axis
<figure><img src="image_11.png" /><figcaption>**Figure 2.** Convergence of simulations assessed by repeating simulations. The Amber96 force field in an implicit solvent is employed here. Results for two peptides of REMD simulations repeated with five runs. For each run, the density profile of the RMSD to the experimental reference structure is shown on the left. The density profile of the five concatenated trajectories is shown in black. The free energy map is also obtained by the assembly of all five runs. The map is clustered with the *k*-means algorithm, and the centroid structures are superposed to the experimental reference structure (in red). The two first residues are indicated with red and gray balls, respectively. The size of each cluster is indicated below the 3D structures. For the other peptides, the results are shown in Supporting Information Figure S1.</figcaption></figure>

and the RMSD to the experimental reference structure calculated over the backbone atoms on the y-axis. This projection gives a good idea where the highest populated conformations are situated compared with the experimental reference structure. The radius of gyration is a good measure of the overall shape of the cyclic peptide ring. The second projection is obtained through a principal component analysis (PCA) of the cyclic peptide backbone atom positions of the molecular dynamics trajectory (see below for more details). We employed the "MDanalysis" python package⁹¹ for the PCA calculation.

For the identification of low free energy clusters, we employed two techniques: *k*-means clustering⁹² and HDBscan clustering.⁹³ See *Supporting Information* for more details.

**Prediction of the Experimental Reference Structure.**

Although our molecular dynamics simulations yield an ensemble of conformations and not just a single structure, the most populated cluster of conformations is expected to be near the experimental reference structure. This can be evaluated immediately by looking at the 2D maps that project the conformational landscape onto the RMSD to the experimental reference structure, by reporting the RMSD of the centroid structure of the most populated cluster, i.e., the one with the lowest free energy in general (not always, as the HDBscan clustering algorithm can merge together two clusters of low free energy that do not include the lowest free energy peak).

In real cases, the experimental reference structure is in general not available. Therefore, we set up a protocol to predict it without the use of the RMSD/radius of gyration projection but by performing a principal component analysis (PCA) directly from the coordinates of the backbone atoms sampled during the MD simulation (REMD or ST protocol). To perform a PCA, the MD trajectory needs first to be aligned to one reference frame. We simply choose the first frame as the reference frame. From the PCA, we retained as many components as necessary to obtain a cumulated variance of 80% of the total variance. In our case, about four to six principal components (PC) had to be retained. We clustered this multidimensional PCA space with the HDBscan algorithm using the same parameters as described above. The PCA space is projected onto the first two principal components to give the

2D maps shown in Figure S4. The centroid of the most populated and lowest free energy cluster is indicated by a green star. In the rare case that the position of the lowest free energy value in the PCA maps is different from the position of the most populated cluster, we plotted a cyan star for the former and a yellow star for the latter. For each star, its RMSD to the experimental reference structure and radius of gyration are calculated, and this coordinate is indicated again by a star of the same color on the 2D maps shown in Figure 9.

We also tested a variant of this protocol, where the frame number corresponding to the centroid of the most populated cluster from the first PCA run is used as a reference frame to realign the MD trajectory and to perform a second PCA run, which should be less dependent on the arbitrary choice of the first frame as the reference frame from the first run. As the results did not change significantly (results not shown), we retained the first version of our protocol.

**Combining Different Force Fields and Simulation Protocols.** As will be shown in the Results and Discussion section, each force field or simulation protocol yields a slightly different free energy map for the same peptide. Even though different force fields cannot be combined directly, their resulting free energy maps can be combined to obtain consensus free energy maps which should give a more complete view of the conformational space of a given peptide. We combined simulations having either a different force field or a different simulation protocol (REMD vs ST, implicit vs explicit solvent), or both. These combinations were done by concatenating the trajectories of the simulations. As the number of frames of each trajectory is not the same for each simulation, we determined its minimum value and subsampled the other trajectories with the same number of frames, which yielded equal weights of each simulation that contributed to the concatenated trajectory. As before, we constructed the free energy maps from the concatenated trajectories and applied the same protocol as the one described above for the prediction of the experimental reference structure. All possible combinations over the nine simulations were done by joining together two, three, or more simulations.

The performance of each combination in the prediction of the experimental reference structure is quantified by a single
value, the median RMSD of the first clusters among all nine peptides of our data set.

## RESULTS AND DISCUSSION

### Convergence Assessed by Repeating Simulations.

Usually, convergence of molecular dynamics simulations can be assessed by repeating the same simulation with different initial conditions (initial velocities, for example). We assessed the convergence of REMD simulations, as shown for the Amber96 force field for two peptides in Figure 2. For the other peptides, the results are shown in Supporting Information Figure S1. We compared the density profiles of each of the 5 runs of the RMSD to the experimental reference structure to the average density profile. For almost all peptides, the density profiles converged at the end of the simulation (1000 ns, or more exactly 900 ns as we reject the first 10% to avoid bias from initial conditions), as shown in Figure 2 for peptide 8.C (see Figure S1 for the other peptides). In terms of convergence, the worst case was peptide 10.A, but as can be seen from Figure 2, the ranking of the three major clusters is not affected by the slow convergence; only the absolute values of the cluster sizes did not converge at the end of the simulation. The speed of convergence can be assessed quantitatively by tracing the evolution of the profile–profile distance measured by the Jensen–Shannon divergence (JSD) (Figure 3). For the 8C peptide, all five simulations converged

<figure><img src="image_12.png" /><figcaption>**Figure 3.** Speed of convergence is shown with the evolution of the Jensen–Shannon divergence (JSD) of each RMSD density profile to the average RMSD density. Same REMD simulations as in Figure 2. For the other peptides, the time evolution of the JSD values is shown in Figure S2b. Below the threshold of a JSD of 0.05, two RMSD profiles are nearly identical. Convergence is obtained if all curves fall and stay below this threshold.</figcaption></figure>

after only 300 ns, as all Jensen–Shannon divergences are under the threshold that we defined (0.05) and that correspond to very high similarity of the density profiles. Overall, the six smaller peptides 7A to 8C show a much faster convergence than the three larger peptides 9A to 10B, as demonstrated by the JSD curves in Figure S2b and the first line of mean JSD values at the end of the trajectories in Figure 4.

As the comparison of the 8C and 10A peptide results shows (Figure 2), the convergence is not a guarantee of obtaining the

<figure><img src="image_13.png" /><figcaption>**Figure 4.** Individual convergence of REMD and ST simulations, as measured by the mean JSD of all runs at the end of the trajectories (i.e., the mean of the last points of the curves in Figures 3 and S2).</figcaption></figure>

highest population for the cluster with the lowest RMSD to the experimental reference structure. Even though the REMD simulations of the 8C peptide converged very rapidly, all five runs yielded the highest population in a cluster more than 3 Å away from the reference structure. On the contrary, the five nonconverged runs of the 10A peptide all ranked the lowest RMSD cluster at the first position in terms of free energy/ population size.

For the ST simulations, we assessed the convergence by repeating the simulations with the Amber96 force field with two runs, runs 1 and run 2. Figures 4 and S3 show that the convergence is also very good among almost all peptides. Only one case showed a poor convergence (peptide 9A). This case is a bit surprising, as the cluster that accumulated 58% of the trajectory frames in run 1 did not even appear with a small percentage on run 2. It seems that run 1 explored a cluster that has not been explored by run 2, whereas the other two clusters are both explored by each run. Among both runs, we selected run 1 for further analysis which is less favorable in terms of the ranking of the clusters, so the performance of ST with Amber96 is a bit lower rated for 9A. Unlike the REMD results, the convergence speed does not show a clear dependence on peptide size with the ST protocol. As can be seen in Figure S2a, the convergence speed is quite different among the nine peptides and also among the peptides with the same number of amino acids: peptide 7A shows a horizontal line above the 0.05 JSD threshold. 7B converges very rapidly, and 7C needs about 4000 ns simulation time to converge both curves together. For the three larger peptides of eight amino acids, the convergence speed shows a similar heterogeneity. For the largest peptides of ten amino acids (10A and 10B), the convergence is average for 10B (JSD < 0.05 after 4500 ns) and slow for 10A due to a horizontal line like in the 7A case. Nevertheless, a JSD threshold of 0.1 instead of 0.05 would classify both cases, 7A and 10A, as being converged after 6000 ns.

The individual convergence of REMD Amber96 and ST- Amber96 is summarized in Figure 4 with the mean JSD values at the end of the trajectories. The convergence of the REMD Amber14 and RSFF2C simulations is shown in the same figure. The JSD values are higher here than for the Amber96 force field, especially for the larger or more flexible peptides (8C– 10B) but are still relatively contained. We cannot explain why Amber14 and its derived RSFF2C force fields yield a lower convergence among repeated REMD simulations than the Amber96 force field. Maybe the fact that Amber14 and RSFF2C were developed for an explicit solvent plays a role here.

### Agreement across Different Simulation Conditions.

Even though each simulation protocol taken alone shows overall very good convergence properties, the agreement among different simulation protocols and force fields is not as good. A quick comparison of the results shown in Figures 2 and S3 for peptide 8C illustrates this fact: while the most populated cluster of the REMD simulations is at about 3.3 Å from the experimental reference structure, the ST simulations yield a better ranking with the most populated cluster at about 1.7 Å. This is more due to different population sizes of each cluster than to “shape” differences of the free energy maps, i.e., the explored regions. The positions of the clusters of the free energy maps are quite similar among REMD and ST, but each simulation stays a different time in each cluster. To have a more quantitative analysis of the interprotocol agreement, we traced the evolution of RMSD profile–profile JSD distances
<figure><img src="image_14.png" /><figcaption>**Figure 5.** Speed of REMD with ST agreement, as measured by the JSD of the RMSD profile of a simulation protocol (ST in blue and REMD in orange) to the averaged RMSD profile (combining ST and REMD profiles). The results obtained with the Amber96 force field are shown here.</figcaption></figure>

comparing both protocols, as we compared the runs inside a single protocol (intraprotocol convergence). Therefore, we merged all available runs into one large trajectory per protocol and took the same number of frames per protocol by subsampling the trajectory with a higher number of frames. With the data from each protocol, we generated an average RMSD profile which each protocol should reproduce. The JSD curves in Figure 5 show that four peptides (7A, 7C, 8A, and 8B) showed a good agreement between the REMD and ST protocols, whereas the others, mainly larger peptides, do not. For the cases with lower agreement (7B, 8C, 9A, 10A, and 10B), the curves are almost horizontal at the second half of the trajectories, indicating that longer simulation times would not help here to obtain a better agreement among both protocols.

To summarize the number of cases where each protocol converges independently and taken together, we calculated the maximum JSD value of the end points of each JSD curve, as a metric to assess intraprotocol convergence or interprotocol agreement at the end of the simulated trajectories. The bar plot shown in Figure 6 confirms that convergence is obtained more

<figure><img src="image_15.png" /><figcaption>**Figure 6.** Final convergence or agreement, as measured at the end of the trajectories by the maximum JSD of the RMSD profile of a run or simulation protocol to the averaged RMSD profile. The results obtained with the Amber96 force field are shown here. The individual convergence of each simulation protocol (ST in blue and REMD in orange) is compared to the agreement of both taken together ("ST vs REMD" in green). The time evolution of the JSD values is shown in Figure S2.</figcaption></figure>

often by simply repeating a simulation with the same protocol (REMD or ST here) than by doing a simulation with different MD simulation methods, such as REMD and ST here. Nevertheless, both methods agreed quite well, especially for the six smaller peptides (7A–8C). The three larger peptides (9A, 10A, and 10B) challenged intraprotocol convergence and interprotocol agreement.

For the analyses shown up to here, it has to be reminded that a single value, as here the maximum JSD value at the end point, is quite reductive to describe the similarity of two conformational ensembles. As we will see below, the visual comparison of the 2D free energy maps is more informative. As an example, while peptide 7B showed a high maximum JSD of 0.16 in the comparison ST vs REMD in Figure 6, the difference is mainly due to high RMSD clusters present in the REMD simulations but not in the ST ones for the Amber96 force field (see Figure 9). The highest populated clusters were identical for both simulation protocols. On the other hand, peptide 9A had the same maximum JSD value as peptide 7B, but the comparison of the free energy maps shows that the REMD and ST maps of peptide 7B for the Amber96 force field look more similar than the ones of peptide 9A. This visual difference in 2D is not well represented by a single JSD value.

Doing the same analysis as in Figure 6 to quantify the interprotocol and force field agreement for all tested simulation conditions (MD protocol and force field), we obtained nine heatmaps, as shown in Figure 7, one for each peptide. The average and standard deviation over all nine peptides are shown in Figure 8. The individual heatmaps of each peptide should help in the comparison of the free energy maps shown in Figure 9. The order of the columns is the same in Figures 7–9. The first five columns correspond to the REMD protocols, and the last four columns correspond to the ST protocols (as indicated by “ST-” before the force field name, absence of it means the REMD protocol). For each protocol, one simulation is done with an explicit solvent (columns 4 and 8, as indicated by “-exp”), whereas the others are done with an implicit solvent. Figure 7 shows the singular behavior of the Charmm36m force field in the ST protocol with the high JSD values in the corresponding “ST-Charmm36m” row or column.
<figure><img src="image_16.png" /><figcaption>**Figure 7.** Interprotocol and force field agreement, as measured by the mean JSD of the RMSD profiles to the averaged RMSD profile. The mean JSD is calculated here at the end of the trajectories, as a measure of the agreement obtained at the end of the REMD or ST simulations. Agreement between force fields and simulation protocols (ST where indicated, REMD otherwise). “exp” stands for the explicit solvent, and all other simulations were done in an implicit solvent.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 8.** By averaging through the heatmap values of the nine peptides shown in Figure 7, the “average JSD” values are reported here on the left with their respective standard deviations on the right.</figcaption></figure>

For peptide 8A, the two REMD simulations in an implicit solvent with the Amber14 and RSFF2C force fields show quite high JSD values compared to the other simulations, as will be discussed in the next section analyzing the free energy maps of Figure 9. For peptide 7C, the RMSD profile of the explicit solvent REMD simulation with the RSFF2C force field is only similar to the other explicit solvent simulation in our test set: ST-Amber14-exp, even though the force field and the simulation protocol are different.

Overall, Figure 8 in comparison to Figure 4 confirms that agreement across force fields and/or simulation protocols (i.e., interprotocol) is lower than the convergence of repeated simulations under the same conditions (i.e., intraprofile). The lowest JSD values, i.e., best agreement, are observed for

the Amber96 force field, comparing REMD with the ST simulation protocol (JSD = 0.12 ± 0.07 (mean ± SD)). The fact that the RSFF2C force field is derived from the Amber14 force field is confirmed by the second best JSD value of 0.14 ± 0.07, comparing the REMD simulations with Amber14 and RSFF2C in an implicit solvent. The two explicit solvent simulations, REMD RSFF2C-exp and ST-Amber14-exp, hold the third best JSD value of 0.14 ± 0.04, hinting to the visual similarity of the free energy maps obtained in the two corresponding columns in Figure 9.

The switch in the MD protocol (REMD ↔ ST) had overall a lower impact (average JSD values: 0.12 ± 0.07 for Amber96 and 0.18 ± 0.11 for Amber14) than the switch of the force field family from Amber96 to Amber14 (or RSFF2C): average
<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 9.** Free energy maps projected onto the radius of gyration (horizontal axis) and the RMSD to the reference structure. The units of the axes are in Å and the free energy unit is in kJ/mol. The first five columns were obtained with the REMD protocol and the last four columns with the ST protocol. The force field employed is indicated on top of each column. RSFF2C-exp and ST-Amber14-exp are explicit solvent simulations. The other columns are implicit solvent simulations. The clustering results are indicated by the cluster ids and their fraction of the total number of frames, from the most populated cluster (id = 1) to the less populated cluster (clusters with a size below 2% are not shown). The centroid of each cluster is indicated by white dots. The centroid position of the most populated and lowest free energy cluster of the PCA projection (see Figure S4) is indicated by a green star and green dashed lines. In the rare case that the position of the lowest free energy value in the PCA maps is different from the position of the most populated cluster, we plotted a cyan star for the former and a yellow star for the latter. The radius of gyration of the reference experimental structure is indicated by a black cross on the 0.0 RMSD line.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Figure 10.** (a) RMSD to the reference structure for the best ranked cluster for each peptide. Results shown for each peptide. The last bar is the combination of Amber96, Amber14, ST-Amber96, and ST-Amber14. (b) Comparison of the two simulation methods, REMD and ST, for three force fields Amber96 (blue), Amber14 (orange), and Charmm36m (purple) in an implicit solvent. The RMSD values of the lowest free energy cluster obtained in either method are compared by subtracting the RMSD value obtained by REMD from the one obtained by ST. Positive values indicate cases where REMD gives better results, and negative values indicate cases where ST gives better results. The green bars compare the explicit solvent ST-Amber14 simulation to its implicit solvent equivalent. The red bars compare the explicit solvent REMD-RSFF2C simulation to its implicit solvent equivalent. For the green and red bars, positive values indicate cases where the implicit solvent gives better results and negative values indicate cases where the explicit solvent gives better results.</figcaption></figure>

JSD of 0.23 ± 0.09 for REMD and 0.19 ± 0.10 for ST. The switching from an implicit to an explicit solvent had a greater impact on the RSFF2C force field with REMD than on the Amber14 force field with ST, as quantified by the JSD values of 0.24 ± 0.08 (RSFF2C ↔ RSFF2C-exp) and 0.17 ± 0.09 (ST-Amber14 ↔ ST-Amber14-exp).

**Free Energy Maps Obtained with REMD and ST on Several Force Fields.** Figure 9 shows the free energy maps of seven different simulation conditions over all nine peptides of our data set. Overall, the three force fields Amber96, Amber14, and the Amber14 variant, RSFF2C, produce similar free energy maps, whereas the Charmm36m force field gives quite different results.

The observation of multiple free energy minima in almost each map shows that the cyclic peptides of our data set have some conformational flexibility, despite the fact that all of them, except peptide 8C, were designed to have high conformational stability. Only the smallest cyclic peptides of our data set with seven amino acids (7A–C) show in most cases only one wide free energy minimum. It is difficult to prove that the alternative free energy minima with high RMSD distances to the experimental reference structure exist in vitro, as they have not been observed in the NMR studies that produced the reference structures. But some of these alternative free energy minima were observed across several simulation conditions (MD protocols and force fields), which should give strong evidence for their existence. For example, peptide 8B has a low free energy region centered at an RMSD of 1.7 Å and an $R_G$ of 4.0 Å that is observed in all but two (Amber14 and RSFF2C) simulation conditions. The more flexible peptide 8C has a low free energy “band” spanning a large interval of RMSD values from 0.2 Å to 3.5 Å almost without interruption, i.e., white areas in the maps. Several free energy minima are observed in this “band”, across all simulation conditions. The larger peptide 10B shows a similar free energy landscape, with the difference that the high RMSD

regions are shifted to lower $R_G$ values than that of the experimental reference structure. Only the two explicit solvent simulations (RSFF2C-exp and ST-Amber14-exp) do not sample significantly the highest RMSD region above 2.8 Å and the lowest RMSD region below 1 Å; they are both centered in a more restricted RMSD interval of [1, 2.8]Å, compared to the implicit solvent simulations that span a large RMSD interval of [0.5, 3.5]Å. This may be explained by the slowdown of the conformational changes of the peptide due to explicit water molecules, as found already in ref 94. But for the other ten residues peptide 10A, the opposite is observed when comparing ST-Amber14 with ST-Amber14-exp. The latter spans higher RMSD regions.

Looking at the impact of the MD protocol (REMD vs ST), one observes that overall, the maps are quite similar for a given force field, showing that the simulations are well converged and that the sampling method does not greatly influence the obtained free energy maps. Some small differences can be observed nevertheless: For some peptides, regions with high RMSD conformations compared with the experimental reference structure are more populated by the REMD simulations compared to the ST simulations. This can be seen the best on peptides 9A (Amber96 only) and 10A (Amber96 and Amber14). The lowest free energy cluster is very similar for a given force field, except for some cases where this cluster is lower in RMSD for one method compared to the other. REMD yields a lower RMSD for peptides 9A and 10B (Amber14 both) and ST yields a lower RMSD for peptides 8A (Amber14), 8C, and 10B (Amber96), see Figure 10b.

**Prediction of the Experimental Reference Structure.**

Applying the protocol described in the Materials and Methods section to predict the experimental reference structure from our MD simulations, we obtained the best clusters indicated by the green stars and green dashed lines in Figure 9. The green stars correspond to the most populated and lowest free energy cluster in the PCA maps (Figure S4) which are constructed
without knowledge of the experimental reference structure. As can be seen in Figure 9, the positions of the green stars in the RMSD/radius of gyration projection superpose very well with the global free energy minima of these maps. As described in the Materials and Methods section, in the case that the position of the lowest free energy value in the PCA maps is different from the position of the most populated cluster, we plotted a cyan star for the former and a yellow star for the latter. This rare case happened only in one map, namely, for peptide 8B with the RSFF2C-exp force field: the cyan star gave a much lower RMSD compared to the yellow star (0.5 Å vs 1.8 Å). In the further analysis, we retained the position of the yellow star for peptide 8B, as this single case cannot validate one prediction method over the other.

In Figure 10a, we report the RMSD values of the green (or yellow) stars for each MD simulation. To see the full context of these RMSD values, one should also keep a look on the full free energy maps shown in Figure 9. We remind here that the experimental reference structure is the best representative conformer of the NMR ensemble (first model here). In addition, we also calculated the minimum RMSD of the best ranked cluster to the 20 NMR models of each peptide, see Figure S5. As this only reduces globally the RMSD values without changing the profiles of 10a, we kept here the RMSD values to the best representative conformer of the NMR ensemble.

For the smallest cyclic peptides of the set, with seven amino acids (7A–7C), the different simulations converged to similar RMSD values (Figure 10a). Comparing the RMSD values among these three peptides 7A to 7C, one can observe that peptide 7A has a 3-fold higher RMSD than peptide 7C (1.5 Å vs 0.5 Å). This shows that the notion of “low” RMSD is specific to each peptide.

For the other peptides, only peptide 10A shows very good agreement among the force fields and simulation protocols: All simulations, except the two using Charmm36m, yield the same global free energy minimum at only 0.4–0.5 Å (Figure 10a), despite the relatively large size of peptide 10A of 10 amino acids, which allows a large range of conformations with up to 4 Å of RMSD to the reference structure (Figure 9).

The performance on peptide 10A is even more impressive if it is compared to the second peptide with 10 amino acids in the set: peptide 10B. Here, the lowest RMSD is only 1.0 Å obtained with two simulations: ST-Amber96 and REMD Amber14 (Figure 10a). Nevertheless, all simulations sampled with a significant frequency conformations with a low RMSD near 1.0 Å (see Figure 9). The differences among the simulations in the population sizes of the clusters may come from insufficient convergence, as REMD and ST show quite different results, even if the same force field is used (Figure 10b). It may come also simply from a relatively flat free energy landscape for peptide 10B, which does not give a clear preference for the experimental reference conformation.

The results obtained on peptide 8A are quite similar to those for peptide 7C, except for the REMD Amber14 and RSFF2C simulations in implicit solvent. The other seven simulations of peptide 8A give a very low RMSD of 0.4–0.6 Å (Figure 10a) and do not sample much the higher RMSD regions (Figure 9). The REMD Amber14 and RSFF2C simulations sample the same low RMSD regions as the other simulations, but at a lower frequency than the unique cluster at 2.2 Å (Figure 9), which is not sampled at all by the other simulations. We investigated this high RMSD cluster of peptide 8A in depth,

see Supporting Information Section S8. In summary, the high RMSD cluster seems to be due to the ψ backbone dihedral force field parameters of Amber14 and by extension also RSFF2C which is derived from Amber14. A joint swap of the ψ angles of the first two residues (1ASP and 2ASP) has to happen to bring the high RMSD structures nearer to the reference structure (see Figure S10). We replaced in Amber14 the ψ angle force field parameters with the ones from Amber96, allowing the “Amber14 mod” REMD simulation to escape the high RMSD basin with the same speed as with Amber96 (see Figure S10a,e). Alternatively, higher maximum temperatures in REMD allow the simulations using the original Amber14 force field to escape also the high RMSD basin (see Figure S15).

For peptide 8B, the Amber14 force field and its derivative RSFF2C give the best results in terms of low RMSD with a good value of 0.9 Å. The use of an explicit solvent degraded this good result here for both force fields (Figure 10b). Here, the two simulation protocols, REMD and ST, give exactly the same results (Figure 10b), and also the free energy maps look very similar (Figure 9).

Peptide 8C is particular as it contains only L-amino acids and is therefore probably more flexible than the other peptides that contain a mixture of L- and D-amino acids. In fact, the free energy maps of 8C show a large landscape from 0.5 Å to 3.5 Å of RMSD (Figure 9). For such a flexible peptide, it is difficult to rank the lowest RMSD cluster in the first position. The Amber14 force field picks up the right cluster, with a very low RMSD of 0.4 Å. The use of an explicit solvent degrades this performance, as the ST-Amber14 explicit solvent simulations give a high RMSD cluster of 2.5 Å as the most populated cluster. Interestingly, for the RSFF2C force field, it is the opposite case (Figure 10b): the explicit solvent improved the ranking with the lowest observed RMSD for 8C of 0.3 Å, instead of a high RMSD of 2.7 Å for the implicit solvent simulation. The use of the Amber96 force field gives here even worse results, with 3.4 Å using the REMD protocol. The ST protocol also yields a quite low free energy basin at 3.4 Å, but the basin at 1.5 Å has a slightly lower free energy, explaining the great difference of the two protocols using the Amber96 force field for peptide 8C (Figure 10a). Both protocols result in a quite similar free energy map (Figure 9, peptide 8C), showing that the differences observed in Figure 10a should not be overinterpreted; a look on the full free energy maps is still necessary.

The last peptide to be analyzed is peptide 9A: Here, the REMD Amber96 simulation gives the best results of 0.6 Å. The Amber14 force field give a doubled RMSD value (1.2 Å) for the REMD protocol and for the ST protocol, but only in an explicit solvent. The ST-Amber14 implicit solvent simulation ranks the higher RMSD cluster of 2 Å at the first position. The RSFF2C force field gives almost identical free energy maps as with the Amber14 force field using the REMD protocol. The explicit solvent RSFF2C simulation populates the 0.6 Å cluster a bit more than its implicit solvent counterpart (Figure 9).

**Explicit vs Implicit Solvent.** The Amber14 force field has been optimized for the use of an explicit solvent; therefore, we expected better results using this force field with an explicit solvent than with an implicit solvent. Especially the comparison of the ST-Amber14 implicit solvent simulations with its explicit solvent counterpart (ST-Amber14-exp column in Figure 9) should show the improvements by the use of an explicit solvent vs an implicit solvent. But surprisingly, this is
<figure><img src="image_21.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(b)</figcaption></figure>

**Figure 11.** (a) For each simulation, the RMSD values of the nine peptides shown in Figure 10a are assembled into a box plot. The simulations are sorted here with increasing median RMSD values, shown with green lines. (b) All 15 combinations of four implicit solvent simulations without ST-Charmm36m.

observed only for peptides 7A (lowest free energy cluster does not improve, but lower RMSD areas are explored; see Figure 9) and 9A (negative green bar in Figure 10b). If the comparison for peptide 9A is done against the other implicit solvent simulations, the lowest RMSD cluster is as well or better ranked with the three REMD simulations (Amber96, Amber14, and RSFF2C), showing that an explicit solvent is not a requirement to obtain a good ranking. On the contrary, we observed more cases where the use of an explicit solvent degrades the ranking for ST-Amber14, namely, 8B, 8C, and 10B (positive green bars in Figure 10b).

**RSFF2C Force Field.** Finally, we turn to looking at the free energy maps obtained with RSFF2C which is a derivative of the Amber14 force field. The MD simulations are done in implicit and explicit solvents using the REMD protocol (RSFF2C and RSFF2C-exp, respectively). For most peptides, the results obtained with RSFF2C are quite similar to the results obtained with Amber14 in an implicit solvent for both sets, REMD and ST (see Figure 9). Peptide 8A is a special case, as only the two REMD simulations with Amber14 and RSFF2C found a large cluster at a high RMSD of 2.2 Å. This special case underlines the similarity of the two force fields

Amber14 and RSFF2C, at least in implicit solvent conditions. The RSFF2C simulations in an explicit solvent (RSFF2C-exp) gave better results for three peptides (8A, 8C, and 10B; see negative red bars in Figure 10b) and gave only for two peptides, 7C and 8B (positive red bars in Figure 10b), worse results compared to the implicit solvent simulation (RSFF2C). This is the opposite result for what we found when comparing implicit vs explicit solvent for the Amber14 force field (see above), despite the fact that the RSFF2C force field is derived from the Amber14 force field. It seems therefore more important to have an explicit solvent with the RSFF2C force field than with the Amber14 force field.

**Combining Different Force Fields and Simulation Protocols.** As described in the Materials and Methods section, we combined different simulations to achieve more robust results, especially in the prediction of the experimental reference structure. The rationale behind this approach is the observation that while the low RMSD clusters are often the same among the different simulations, the high RMSD clusters are more spread. Combining the simulations should increase the populations of the low RMSD clusters while decreasing those of high RMSD clusters and therefore increasing the
<figure><img src="image_23.png" /><figcaption>**Figure 12.** Combining simulations to improve the predictions of the reference structure. See Supporting Information figures for more details. (a) Predictions for all possible combinations of simulations. One point corresponds to one simulation (first column), two combined simulations (second column) or any higher number of combined simulations (columns three to nine). The y-value of each point is the median RMSD to the reference structure calculated over all nine peptides. The violin plots show the distribution of the points along the y-axis. (b) Same plot as in (a), but without the ST-Charmm36m simulation. (c) Same plot as in (a), but without the ST-Charmm36m simulation and with implicit solvent simulations only.</figcaption></figure>

chance to rank the low RMSD clusters at the first position. As we have already seen in Figure 10a, none of the simulations gives a good prediction for all nine peptides of our data set. This can also be seen in Figure 11a, where the RMSD values of the nine peptides shown in Figure 10a are assembled into a box plot per simulation. The upper limit extends often above 2.0 Å. The single simulations are sorted here by the median RMSD of the nine RMSD values of each peptide.

The median RMSD values of Figure 11a are reported in the first column in Figure 12a. If only a second simulation is added to one of the nine simulations, the median RMSD values are significantly improved, as can be seen by comparing the violin plots of the first two columns in Figure 12a. Figure S18 gives the details of the 36 possible combinations of twins among the nine simulations (binomial coefficient: $\binom{n}{k} = \frac{n!}{k!(n-k)!} = 36$ with $n = 9$ and $k = 2$).

Looking at the median RMSD values in Figure S18, represented by green bars, the top 3 twin combinations are Amber96 + RSFF2C-exp, RSFF2C-exp + ST-Amber14, and RSFF2C + ST-Amber14. As many twin combinations are clustered around 0.9 Å for the median RMSD (see also Figure 12a), other criteria might be used in addition to rank or select twin combinations, for example, the upper and lower limits of the box plots in Figure S18. Here, the Amber14 + ST-Amber14 twin shows a good compromise with a low upper limit and a not too high lower limit in the RMSD scale. But due to the high number of 36 twin combinations and the similarity of the values, this selection or the one of the top 3 should not be used as a general rule. What is more important is the result that the combination of two simulations is for the majority better than a single simulation taken alone. Combining more than two simulations continues to improve the RMSD distributions in Figure 12a, especially by reducing the upper limit of these distributions. Combining all nine simulations yielded a median RMSD as low as the best single simulation (here ST-Amber14). As in another set of peptides, the ST-Amber14 might not be the best simulation; the combination of diverse simulations yields robustness in the choice of the right simulations.

Instead of trying to identify the single best simulation among our set of simulations, we can exclude the worst simulation, namely, ST-Charmm36m which has an important gap in median RMSD with respect to the other eight simulations (see

Figure 11a). Unsurprisingly this improves the RMSD distributions for all simulation combinations (see Figure 12b).

Producing eight REMD/ST simulation trajectories per peptide is quite costly, especially for the two explicit solvent simulations. Fortunately, discarding both in addition to ST-Charmm36m does not significantly degrade the RMSD distributions (see Figure 12c). There are essentially less lowest and highest RMSD values, probably due to the reduced number of combinations. For example, for four combined simulations, there are 70 combinations among eight possible simulations (Figure 12b) but only 15 combinations among six possible simulations (12c). Producing four implicit solvent REMD and ST simulations may be a good compromise in terms of robustness vs calculation cost, as the median RMSD values are all below 1.0 Å.

All 15 combinations of implicit solvent simulations are listed in Figure 11b in the order of increasing median RMSD values. The lower and upper bounds of the second best combination are lower than the combination ranked in the first position. Despite its slightly higher median RMSD value, it is one of the best choices among the 15 combinations. The second best combination is Amber96, Amber14, ST-Amber96, and ST-Amber14. The combination at the seventh position is composed of only REMD simulations: Amber96, Amber14, RSFF2C, and Charmm36m. It yields an even lower upper bound and may be a viable choice to reduce the time to set up two different simulation protocols, namely, REMD and ST. Even if RSFF2C and Charmm36m are ranked lower as single simulations (see Figure 11a) than ST-Amber96 and ST-Amber14, they nevertheless enrich the diversity in free energy maps, necessary to obtain robust predictions of the reference structure.

Overall, this shows that a limited number of implicit solvent simulations with the two Amber force fields, Amber96 and Amber14, and with the two simulation methods REMD and ST yield already robust results, predicting the experimental reference structure accurately for 8 out of 9 cyclic peptides. Only for peptide 9A, the best cluster of the Amber96 + Amber14 + ST-Amber96 + ST-Amber14 combination is at a higher RMSD of 2.0 Å, while the other eight peptides are below 1.0 Å or at 1.5 Å for peptide 7A, see the last cyan bars of "10: combination" in Figure 10a. Just looking at the RMSD values of the selected best clusters in Figure 10a does not explain the high RMSD value for peptide 9A. It has to be
reminded here that the RMSD values of the four combined simulations are not obtained by simple averaging over the RMSD values of the four individual simulations but by combining the MD trajectories and from the PCA analysis of the combined trajectory. Looking at the free energy maps of the individual and combined simulations gives a clearer picture, see Figure S6. The two Amber14 simulations (REMD + ST) have an important population on the cluster of 2.0 Å to the reference structure, especially for ST with 91%. The two Amber96 simulations (REMD + ST) show both a cluster at low RMSD, especially for REMD with 70%, but this is still lower than that for the high RMSD cluster. At the end, the high RMSD cluster has 31% and the low RMSD cluster 21% in the combination of the four simulations. But at least the difference in population between both clusters is much lower than that for both Amber14 simulations, where the same low RMSD cluster is almost missing completely.

Finally, we mention that as the peptides are not simulated in interaction with a binding partner, there is no need here to use explicit solvent simulations to generate accurate free energy landscapes. Three or four implicit solvent simulations will take less computational time than a single explicit solvent simulation and will more crucially give robust and accurate results, which are not achieved with a single explicit solvent simulation.

**Validation on Four Additional Peptides.** We selected four additional peptides as a test set with the following PDB codes: 6UD9, 6UCX, 6UDZ, and 6UDW. The peptides listed in Table S1 are of the same nature as the nine initial peptides of the study, as they have the same number of amino acids, two of 8 A.A. and two of 10 A.A., and they have mixed chirality. All of them were designed peptides from the David Baker group⁹⁵ but solved by X-ray instead of NMR.

We applied our recommendation to use the four simulations, Amber96 + Amber14, ST-Amber96, and ST-Amber14, combined. We produced 1000 ns for each of the eight replicas of the REMD simulations and 10 µs for each ST simulation. Figure S7 shows the 16 free energy maps of each simulation at 300 K, as well as the combined free energy maps for each of the four peptides. One can observe that the two REMD simulations had difficulties in exploring the low RMSD region for 6UD9 and 6UDW, but not for 6UDZ nor for the Amber14 simulation of 6UCX. ST performed better here, except for 6UCX with the Amber96 force field. The ST-Amber14 simulations performed well, but populated an higher RMSD cluster at the first position for 6UDW. The RMSD values of the best ranked clusters are reported in Figure S8. The combination of all four simulations is still beneficial, as the low RMSD region is populated for every peptide, see Figure S7. The low RMSD population is highest for three of the four peptides; only for 6UDW, a higher RMSD cluster is more populated (42%) than the lowest RMSD cluster (17%).

This validation of the combined use of four simulations shows that even if only half of the simulations have a cluster with a low RMSD, the combined simulation can still recover it, as was seen here on the two peptides 6UD9 and 6UCX. This stems from the fact that the high RMSD clusters are often different in their 3D structures and are therefore not clustered together in the PCA maps (Figure S9).

## CONCLUSION

Our study explored various molecular dynamics “ingredients” for the cartography of the conformational landscape of cyclic

peptides and the prediction of their reference experimental structures. Three types of ingredients were investigated: the force field, the solvent model, and the simulation method. For the force fields, no single best force field could be observed; older force fields like Amber96 and Amber14 performed as well as the recent RSFF2C force field. The use of an explicit solvent instead of an implicit solvent showed a mixed picture: for RSFF2C, the prediction of the native structure worked better in an explicit solvent, while for Amber14, we observed the opposite. Under the same force field and solvent model, both simulation methods, REMD and ST, give very similar free energy maps, especially for the lowest free energy clusters. Some differences could be observed in higher free energy clusters, and REMD tended to produce more high free energy clusters than ST. No clear advantage for one method vs the other in the prediction of the native structure could be observed. The intraprotocol convergence of repeated simulations under the same conditions has been observed to be higher than the interprotocol agreement of simulations using different conditions, i.e., force field, solvent model, or simulation method. This finding is important to underline, as convergence of simulations is often only evaluated by repeating simulations under the same conditions, which seems to be an insufficient convergence criterion.

Among the nine peptides, the native structure is predicted with a median backbone RMSD ranging from 0.8 to 1.8 Å with a single MD simulation. The ST-Amber14 (implicit solvent) simulation yielded the best result here. To improve these values, we combined systematically all nine simulations, merging two, three, or more simulations together. The average median RMSD is gradually reduced from 1.15 to 0.85 Å with the number of merged simulations. More importantly, the maximum median RMSD is also greatly reduced from 1.8 to 0.85 Å. This can be explained by the observation that high RMSD clusters are less conserved among the seven simulations than the lowest RMSD cluster. The merging of the simulations enforces therefore the populations of the lowest RMSD cluster, which becomes the most populated cluster more often. The two explicit solvent simulations performed here are not necessary to obtain good predictions, the six implicit solvent simulations are sufficient for this. As the ST-Charmm36m simulation in an implicit solvent had the highest median RMSD of 1.8 Å with an important gap of 0.6 Å to the next best simulation, we excluded it. Among the 15 combinations of four simulations in an implicit solvent, one combination gave a good median RMSD of 0.85 Å with a maximum RMSD below 2.0 Å: REMD Amber96 + Amber14 and ST Amber96 + Amber14. With a slightly higher median RMSD of 0.9 Å, but with a maximum RMSD of only 1.5 Å, the combination of four REMD simulations is also an interesting alternative: Amber96 + Amber14 + RSFF2C + Charmm36m.

The combination of different “ingredients” of MD simulations, namely, the force field and the simulation method, seems to yield more robust results with a better native structure prediction performance. Three or four implicit solvent simulations are already sufficient and will take less computational time than a single explicit solvent simulation. In addition to the combination of complementary MD simulations, our native structure prediction method is based on a PCA projection of the free energy landscape and robustly yields good predictions for the nine peptides. It would be interesting to apply our method also to the structure prediction
of proteins and protein-peptide complexes, as well as cyclic
peptides with chemical modifications, like N-methylation.³⁴

■ ASSOCIATED CONTENT

**Data Availability Statement**

ST simulation protocol: https://github.com/samuelmurail/SST2/tree/main/bin. Trajectory data and python scripts for analysis: 10.5281/zenodo.15112950.

Supporting Information

The Supporting Information is available free of charge at
https://pubs.acs.org/doi/10.1021/acs.jctc.5c01123.

Information about the clustering algorithms; Figures of
the convergence of REMD and ST simulations; PCA
projections of the free energy maps; Results on four
additional peptides; Section S8 with the detailed results
on the peptide 8A; Boxplots of the MD combinations
tested (PDF)

■ AUTHOR INFORMATION

**Corresponding Author**

Dirk Stratmann – Université Paris Cité, CNRS UMR 8251, INSERM ERL U1133, Unité de Biologie Fonctionnelle et Adaptative, BFA, F-75013 Paris, France; Sorbonne Université, Faculté des Sciences et Ingénierie, UFR 925, MNHN, UMR CNRS 7590, Institut de Minéralogie de Physique des Matériaux et de Cosmochimie, IMPMC, F-75005 Paris, France; orcid.org/0000-0002-7373-0181; Email: dirk.stratmann@u-paris.fr

Authors

Samuel Murail – Université Paris Cité, CNRS UMR 8251, INSERM ERL U1133, Unité de Biologie Fonctionnelle et Adaptative, BFA, F-75013 Paris, France

Jaysen Sawmynaden – Sorbonne Université, Faculté des Sciences et Ingénierie, UFR 925, MNHN, UMR CNRS 7590, Institut de Minéralogie de Physique des Matériaux et de Cosmochimie, IMPMC, F-75005 Paris, France

Akli Zemirli – Université Paris Cité, CNRS UMR 8251, INSERM ERL U1133, Unité de Biologie Fonctionnelle et Adaptative, BFA, F-75013 Paris, France

Maud Jusot – Sorbonne Université, Faculté des Sciences et Ingénierie, UFR 925, MNHN, UMR CNRS 7590, Institut de Minéralogie de Physique des Matériaux et de Cosmochimie, IMPMC, F-75005 Paris, France

Fabio Pietrucci – Sorbonne Université, Faculté des Sciences et Ingénierie, UFR 925, MNHN, UMR CNRS 7590, Institut de Minéralogie de Physique des Matériaux et de Cosmochimie, IMPMC, F-75005 Paris, France; orcid.org/0000-0002-4892-2667

Jacques Chomilier – Sorbonne Université, Faculté des Sciences et Ingénierie, UFR 925, MNHN, UMR CNRS 7590, Institut de Minéralogie de Physique des Matériaux et de Cosmochimie, IMPMC, F-75005 Paris, France

Pierre Tufféry – Université Paris Cité, CNRS UMR 8251, INSERM ERL U1133, Unité de Biologie Fonctionnelle et Adaptative, BFA, F-75013 Paris, France

Complete contact information is available at:
https://pubs.acs.org/10.1021/acs.jctc.5c01123

Notes

The authors declare no competing financial interest.

■ ACKNOWLEDGMENTS

We would like to thank Fan Jiang and Wei Kang from Peking
University, Beijing, for sending us their RSFF2C force field
parameters files. This project was provided with computer and
storage resources by GENCI at CINES thanks to the grant
2021-A0100707641 on the supercomputer Occigen. The
RPBS and iPOP-UP platforms also contributed to calculation
time.

■ REFERENCES

(1) Robinson, J. A.; DeMarco, S.; Gombert, F.; Moehle, K.; Obrecht, D. The design, structures and therapeutic potential of protein epitope mimetics. *Drug discovery today* **2008**, *13*, 944–951.

(2) Olmez, E. O.; Akbulut, B. S. Protein-peptide interactions revolutionize drug development. In *Binding Protein*, InTech, 2012.

(3) Gao, M.; Cheng, K.; Yin, H. Targeting protein- protein interfaces using macrocyclic peptides. *Pept. Sci.* **2015**, *104*, 310–316.

(4) Wang, L.; Wang, N.; Zhang, W.; Cheng, X.; Yan, Z.; Shao, G.; Wang, X.; Wang, R.; Fu, C. Therapeutic peptides: current applications and future directions. *Signal Transduction Targeted Ther.* **2022**, *7*, 48.

(5) Ji, X.; Nielsen, A. L.; Heinis, C. Cyclic Peptides for Drug Development. *Angew. Chem., Int. Ed.* **2024**, *63*, No. e202308251.

(6) Lamers, C. Overcoming the shortcomings of peptide-based therapeutics. *Future Drug Discovery* **2022**, *4*, FDD75.

(7) Cowell, S.; Lee, Y.; Cain, J.; Hruby, V. Exploring Ramachandran and chi space: conformationally constrained amino acids and peptides in the design of bioactive polypeptide ligands. *Curr. Med. Chem.* **2004**, *11*, 2785–2798.

(8) Workman, R. J.; Gorle, S.; Pettitt, B. M. Effects of Conformational Constraint on Peptide Solubility Limits. *J. Phys. Chem. B* **2022**, *126*, 10510–10518.

(9) Tran, H. L.; Lexa, K. W.; Julien, O.; Young, T. S.; Walsh, C. T.; Jacobson, M. P.; Wells, J. A. Structure–Activity Relationship and Molecular Mechanics Reveal the Importance of Ring Entropy in the Biosynthesis and Activity of a Natural Product. *J. Am. Chem. Soc.* **2017**, *139*, 2541–2544.

(10) Bhardwaj, G.; Mulligan, V. K.; Bahl, C. D.; Gilmore, J. M.; Harvey, P. J.; Cheneval, O.; Buchko, G. W.; Pulavarti, S. V.; Kaas, Q.; Eletsky, A.; et al. Accurate de novo design of hyperstable constrained peptides. *Nature* **2016**, *538*, 329–335.

(11) Yongye, A. B.; Li, Y.; Giulianotti, M. A.; Yu, Y.; Houghten, R. A.; Martínez-Mayorga, K. Modeling of peptides containing D-amino acids: implications on cyclization. *J. Comput.-Aided Mol. Des.* **2009**, *23*, 677–689.

(12) Hosseinzadeh, P.; Bhardwaj, G.; Mulligan, V. K.; Shortridge, M. D.; Craven, T. W.; Pardo-Avila, F.; Rettie, S. A.; Kim, D. E.; Silva, D.-A.; Ibrahim, Y. M.; Webb, I. K.; Cort, J. R.; Adkins, J. N.; Varani, G.; Baker, D. Comprehensive computational design of ordered peptide macrocycles. *Science* **2017**, *358*, 1461–1466.

(13) Hickey, J. L.; Sindhikara, D.; Zultanski, S. L.; Schultz, D. M. Beyond 20 in the 21st Century: Prospects and Challenges of Non-canonical Amino Acids in Peptide Drug Discovery. *ACS Med. Chem. Lett.* **2023**, *14*, 557–565.

(14) Gentilucci, L.; De Marco, R.; Cerisoli, L. Chemical modifications designed to improve peptide stability: incorporation of non-natural amino acids, pseudo-peptide bonds, and cyclization. *Curr. Pharm. Des.* **2010**, *16*, 3185–3203.

(15) Mas-Moruno, C.; Beck, J. G.; Doedens, L.; Frank, A. O.; Marinelli, L.; Cosconati, S.; Novellino, E.; Kessler, H. Increasing αvβ3 Selectivity of the Anti-Angiogenic Drug Cilengitide by N-Methylation. *Angew. Chem., Int. Ed.* **2011**, *50*, 9496–9500.

(16) Jwad, R.; Weissberger, D.; Hunter, L. Strategies for Fine-
Tuning the Conformations of Cyclic Peptides. *Chem. Rev.* **2020**, *120*,
9743–9789.

(17) Driggers, E. M.; Hale, S. P.; Lee, J.; Terrett, N. K. The exploration of macrocycles for drug discovery—an underexploited structural class. *Nat. Rev. Drug Discovery* **2008**, *7*, 608.
(18) Saha, I.; Dang, E. K.; Svatunek, D.; Houk, K. N.; Harran, P. G. Computational generation of an annotated gigalibrary of synthesizable, composite peptidic macrocycles. *Proc. Natl. Acad. Sci. U.S.A.* **2020**, *117*, 24679–24690.

(19) Li, J.; Yanagisawa, K.; Sugita, M.; Fujie, T.; Ohue, M.; Akiyama, Y. *CycPeptMPDB*: A Comprehensive Database of Membrane Permeability of Cyclic Peptides. *J. Chem. Inf. Model.* **2023**, *63*, 2240–2250.

(20) Merz, M. L.; Habeshian, S.; Li, B.; David, J.-A. G. L.; Nielsen, A. L.; Ji, X.; Il Khwildy, K.; Duany Benitez, M. M.; Phothirath, P.; Heinis, C. De novo development of small cyclic peptides that are orally bioavailable. *Nat. Chem. Biol.* **2023**, *20*, 624–633.

(21) Costa, L.; Sousa, E.; Fernandes, C. Cyclic Peptides in Pipeline: What Future for These Great Molecules? *Pharmaceuticals* **2023**, *16*, 996.

(22) Barman, P.; Joshi, S.; Sharma, S.; Preet, S.; Sharma, S.; Saini, A. Strategic Approaches to Improvise Peptide Drugs as Next Generation Therapeutics. *Int. J. Pept. Res. Ther.* **2023**, *29*, 61.

(23) Hellinger, R.; Sigurdsson, A.; Wu, W.; Romanova, E. V.; Li, L.; Sweedler, J. V.; Süssmuth, R. D.; Gruber, C. W. Peptidomics. *Nat. Rev. Methods Primers* **2023**, *3*, 1–21.

(24) Rossino, G.; Marchese, E.; Galli, G.; Verde, F.; Finizio, M.; Serra, M.; Linciano, P.; Collina, S. Peptides as Therapeutic Agents: Challenges and Opportunities in the Green Transition Era. *Molecules* **2023**, *28*, 7165.

(25) Moulahoum, H.; Ghorbanizamani, F.; Tok, K.; Zihnioglu, F. On the Cyclization of Non-cyclic Peptides for Biological Applications: Inspiration from Naturally Cyclic Peptides. *ChemistrySelect* **2023**, *8*, No. e202301335.

(26) Adiga, U. Trends in Peptide Drug Discovery: Unleashing the Therapeutic Potential of Peptides: A Review. *Int. J. Chem. Biol. Sci.* **2023**, *24*, 171.

(27) Ramadhani, D.; Maharani, R.; Gazzali, A. M.; Muchtaridi, M. Cyclic Peptides for the Treatment of Cancers: A Review. *Molecules* **2022**, *27*, 4428.

(28) Muttenthaler, M.; King, G. F.; Adams, D. J.; Alewood, P. F. Trends in peptide drug discovery. *Nat. Rev. Drug Discovery* **2021**, *20*, 309–325.

(29) Hill, T. A.; Shepherd, N. E.; Diness, F.; Fairlie, D. P. Constraining Cyclic Peptides To Mimic Protein Structure Motifs. *Angew. Chem., Int. Ed.* **2014**, *53*, 13020–13041.

(30) Malde, A. K.; Hill, T. A.; Iyer, A.; Fairlie, D. P. Crystal Structures of Protein-Bound Cyclic Peptides. *Chem. Rev.* **2019**, *119*, 9861–9914.

(31) Kessler, H. Conformation and Biological Activity of Cyclic Peptides. *Angew. Chem., Int. Ed. Engl.* **1982**, *21*, 512–523.

(32) Damjanovic, J.; Miao, J.; Huang, H.; Lin, Y.-S. Elucidating Solution Structures of Cyclic Peptides Using Molecular Dynamics Simulations. *Chem. Rev.* **2021**, *121*, 2292–2324.

(33) Leman, J. K.; et al. Macromolecular modeling and design in Rosetta: recent methods and frameworks. *Nat. Methods* **2020**, *17*, 665–680.

(34) Jusot, M.; Stratmann, D.; Vaisset, M.; Chomilier, J.; Cortés, J. Exhaustive exploration of the conformational landscape of small cyclic peptides using a robotics approach. *J. Chem. Inf. Model.* **2018**, *58*, 2355–2368.

(35) Sugita, Y.; Okamoto, Y. Replica-exchange molecular dynamics method for protein folding. *Chem. Phys. Lett.* **1999**, *314*, 141–151.

(36) Earl, D. J.; Deem, M. W. Parallel tempering: Theory, applications, and new perspectives. *Phys. Chem. Chem. Phys.* **2005**, *7*, 3910–3916.

(37) Voelz, V. A.; Dill, K. A.; Chorny, I. Peptoid conformational free energy landscapes from implicit-solvent molecular simulations in AMBER. *Pept. Sci.* **2011**, *96*, 639–650.

(38) Butterfoss, G. L.; Yoo, B.; Jaworski, J. N.; Chorny, I.; Dill, K. A.; Zuckermann, R. N.; Bonneau, R.; Kirshenbaum, K.; Voelz, V. A. De novo structure prediction and experimental characterization of folded peptoid oligomers. *Proc. Natl. Acad. Sci. U.S.A.* **2012**, *109*, 14320–14325.

(39) Wakefield, A. E.; Wuest, W. M.; Voelz, V. A. Molecular Simulation of Conformational Pre-Organization in Cyclic RGD Peptides. *J. Chem. Inf. Model.* **2015**, *55*, 806–813.

(40) Damas, J. M.; Filipe, L. C.; Campos, S. R.; Lousa, D.; Victor, B. L.; Baptista, A. M.; Soares, C. M. Predicting the Thermodynamics and Kinetics of Helix Formation in a Cyclic Peptide Model. *J. Chem. Theory Comput.* **2013**, *9*, 5148–5157.

(41) Yu, H.; Lin, Y.-S. Toward structure prediction of cyclic peptides. *Phys. Chem. Chem. Phys.* **2015**, *17*, 4210–4219.

(42) Geng, H.; Jiang, F.; Wu, Y.-D. Accurate structure prediction and conformational analysis of cyclic peptides with residue-specific force fields. *J. Phys. Chem. Lett.* **2016**, *7*, 1805–1810.

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

(44) Pietrucci, F.; Marinelli, F.; Carloni, P.; Laio, A. Substrate Binding Mechanism of HIV-1 Protease from Explicit-Solvent Atomistic Simulations. *J. Am. Chem. Soc.* **2009**, *131*, 11811–11818.

(45) Cossio, P.; Marinelli, F.; Laio, A.; Pietrucci, F. Optimizing the Performance of Bias-Exchange Metadynamics: Folding a 48-Residue LysM Domain Using a Coarse-Grained Model. *J. Phys. Chem. B* **2010**, *114*, 3259–3265.

(46) McHugh, S. M.; Rogers, J. R.; Yu, H.; Lin, Y.-S. Insights into How Cyclic Peptides Switch Conformations. *J. Chem. Theory Comput.* **2016**, *12*, 2480–2488.

(47) McHugh, S. M.; Yu, H.; Slough, D. P.; Lin, Y.-S. Mapping the sequence–structure relationships of simple cyclic hexapeptides. *Phys. Chem. Chem. Phys.* **2017**, *19*, 3315–3324.

(48) Slough, D. P.; McHugh, S. M.; Lin, Y.-S. Understanding and designing head-to-tail cyclic peptides. *Biopolymers* **2018**, *109*, No. e23113.

(49) Slough, D. P.; McHugh, S. M.; Cummings, A. E.; Dai, P.; Pentelute, B. L.; Kritzer, J. A.; Lin, Y.-S. Designing Well-Structured Cyclic Pentapeptides Based on Sequence–Structure Relationships. *J. Phys. Chem. B* **2018**, *122*, 3908–3919.

(50) Slough, D. P.; Yu, H.; McHugh, S. M.; Lin, Y.-S. Toward accurately modeling N-methylated cyclic peptides. *Phys. Chem. Chem. Phys.* **2017**, *19*, 5377–5388.

(51) Cummings, A. E.; Miao, J.; Slough, D. P.; McHugh, S. M.; Kritzer, J. A.; Lin, Y.-S. β-Branched Amino Acids Stabilize Specific Conformations of Cyclic Hexapeptides. *Biophys. J.* **2019**, *116*, 433–444.

(52) Fonseca Lopez, F.; Miao, J.; Damjanovic, J.; Bischof, L.; Braun, M. B.; Ling, Y.; Hartmann, M. D.; Lin, Y.-S.; Kritzer, J. A. Computational Prediction of Cyclic Peptide Structural Ensembles and Application to the Design of Keap1 Binders. *J. Chem. Inf. Model.* **2023**, *63*, 6925–6937.

(53) Miao, J.; Descoteaux, M. L.; Lin, Y.-S. Structure prediction of cyclic peptides by molecular dynamics + machine learning. *Chem. Sci.* **2021**, *12*, 14927–14936.

(54) Hui, T.; Descoteaux, M. L.; Miao, J.; Lin, Y.-S. Training Neural Network Models Using Molecular Dynamics Simulation Results to Efficiently Predict Cyclic Hexapeptide Structural Ensembles. *J. Chem. Theory Comput.* **2023**, *19*, 4757–4769.

(55) Jumper, J.; et al. Highly accurate protein structure prediction with AlphaFold. *Nature* **2021**, *596*, 583–589.

(56) Rettie, S. A.; Campbell, K. V.; Bera, A. K.; Kang, A.; Kozlov, S.; De La Cruz, J.; Adebomi, V.; Zhou, G.; DiMaio, F.; Ovchinnikov, S.; Bhardwaj, G. Cyclic peptide structure prediction and design using AlphaFold. *bioRxiv* **2023**, *529956*.

(57) Zhang, C.; Zhang, C.; Shang, T.; Wu, X.; Duan, H. Highfold: accurately predicting cyclic peptide monomers and complexes with AlphaFold. *bioRxiv* **2023**, *554979*.

(58) Kosugi, T.; Ohue, M. Design of Cyclic Peptides Targeting Protein–Protein Interactions Using AlphaFold. *Int. J. Mol. Sci.* **2023**, *24*, 13257.

(59) Abramson, J.; Adler, J.; Dunger, J.; Evans, R.; Green, T.; Pritzel, A.; Ronneberger, O.; Willmore, L.; Ballard, A. J.; Bambrick, J.; et al. Accurate structure prediction of biomolecular interactions with AlphaFold 3. *Nature* **2024**, *630*, 493–500.
(60) Miao, J.; Ghosh, A. P.; Ho, M. N.; Li, C.; Huang, X.; Pentelute, B. L.; Baleja, J. D.; Lin, Y.-S. Assessing the Performance of Peptide Force Fields for Modeling the Solution Structural Ensembles of Cyclic Peptides. *J. Phys. Chem. B* **2024**, *128*, 5281–5292.

(61) Fisher, M. F.; Payne, C. D.; Chetty, T.; Crayn, D.; Berkowitz, O.; Whelan, J.; Rosengren, K. J.; Mylne, J. S. The genetic origin of evolidine, the first cyclopeptide discovered in plants, and related orbitides. *J. Biol. Chem.* **2020**, *295*, 14510–14521.

(62) Meng, E. C.; Goddard, T. D.; Pettersen, E. F.; Couch, G. S.; Pearson, Z. J.; Morris, J. H.; Ferrin, T. E. UCSF ChimeraX: Tools for structure building and analysis. *Protein Sci.* **2023**, *32*, No. e4792.

(63) Case, D. A.; et al. AmberTools. *J. Chem. Inf. Model.* **2023**, *63*, 6183–6191.

(64) Sousa da Silva, A. W.; Vranken, W. F. ACPYPE-Antechamber python parser interface. *BMC Res. Notes* **2012**, *5*, 367.

(65) Onufriev, A.; Case, D. A.; Bashford, D. Effective Born radii in the generalized Born approximation: the importance of being perfect. *J. Comput. Chem.* **2002**, *23*, 1297–1304.

(66) Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. *SoftwareX* **2015**, *1*, 19–25.

(67) Bussi, G.; Donadio, D.; Parrinello, M. Canonical sampling through velocity rescaling. *J. Chem. Phys.* **2007**, *126*, 014101.

(68) Sindhikara, D.; Meng, Y.; Roitberg, A. E. Exchange frequency in replica exchange molecular dynamics. *J. Chem. Phys.* **2008**, *128*, 01B609.

(69) Cecchini, M.; Rao, F.; Seeber, M.; Caflisch, A. Replica exchange molecular dynamics simulations of amyloid peptide aggregation. *J. Chem. Phys.* **2004**, *121*, 10748–10756.

(70) Huang, J.; Rauscher, S.; Nawrocki, G.; Ran, T.; Feig, M.; de Groot, B. L.; Grubmüller, H.; MacKerell, A. D. CHARMM36m: an improved force field for folded and intrinsically disordered proteins. *Nat. Methods* **2017**, *14*, 71–73.

(71) Kollman, P. A. Advances and continuing challenges in achieving realistic and predictive simulations of the properties of organic and biological molecules. *Acc. Chem. Res.* **1996**, *29*, 461–469.

(72) Maier, J. A.; Martinez, C.; Kasavajhala, K.; Wickstrom, L.; Hauser, K. E.; Simmerling, C. ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99SB. *J. Chem. Theory Comput.* **2015**, *11*, 3696–3713.

(73) Jiang, F.; Wu, H.-N.; Kang, W.; Wu, Y.-D. Developments and Applications of Coil-Library-Based Residue-Specific Force Fields for Molecular Dynamics Simulations of Peptides and Proteins. *J. Chem. Theory Comput.* **2019**, *15*, 2761–2773.

(74) Feng, J.-J.; Chen, J.-N.; Kang, W.; Wu, Y.-D. Accurate Structure Prediction for Protein Loops Based on Molecular Dynamics Simulations with RSFF2C. *J. Chem. Theory Comput.* **2021**, *17*, 4614–4628.

(75) Parrinello, M.; Rahman, A. Polymorphic transitions in single crystals: A new molecular dynamics method. *J. Appl. Phys.* **1981**, *52*, 7182–7190.

(76) Nosé, S.; Klein, M. Constant pressure molecular dynamics for molecular systems. *Mol. Phys.* **1983**, *50*, 1055–1076.

(77) *Temperature Generator for REMD-Simulations*. https://virtualchemistry.org/remd-temperature-generator/.

(78) Patriksson, A.; Spoel, D. v. d. A temperature predictor for parallel tempering simulations. *Phys. Chem. Chem. Phys.* **2008**, *10*, 2073–2077.

(79) Eastman, P.; Swails, J.; Chodera, J. D.; McGibbon, R. T.; Zhao, Y.; Beauchamp, K. A.; Wang, L.-P.; Simmonett, A. C.; Harrigan, M. P.; Stern, C. D.; Wiewiora, R. P.; Brooks, B. R.; Pande, V. S. OpenMM 7: Rapid development of high performance algorithms for molecular dynamics. *PLoS Comput. Biol.* **2017**, *13*, No. e1005659.

(80) Chow, K.-H.; Ferguson, D. M. Isothermal-isobaric molecular dynamics simulations with Monte Carlo volume sampling. *Comput. Phys. Commun.* **1995**, *91*, 283–289.

(81) Åqvist, J.; Wennerström, P.; Nervall, M.; Bjelic, S.; Brandsdal, B. O. Molecular dynamics simulations of water and biomolecules with

a Monte Carlo constant pressure algorithm. *Chem. Phys. Lett.* **2004**, *384*, 288–294.

(82) Zhang, Z.; Liu, X.; Yan, K.; Tuckerman, M. E.; Liu, J. Unified Efficient Thermostat Scheme for the Canonical Ensemble with Holonomic or Isokinetic Constraints via Molecular Dynamics. *J. Phys. Chem. A* **2019**, *123*, 6056–6079.

(83) Dorgham, K.; Murail, S.; Tuffery, P.; Savier, E.; Bravo, J.; Rebollo, A. Binding and Kinetic Analysis of Human Protein Phosphatase PP2A Interactions with Caspase 9 Protein and the Interfering Peptide C9h. *Pharmaceutics* **2022**, *14*, 2055.

(84) Park, S.; Pande, V. S. Choosing weights for simulated tempering. *Phys. Rev. E* **2007**, *76*, 016703.

(85) Nguyen, P. H.; Okamoto, Y.; Derreumeaux, P. Communication: Simulated tempering with fast on-the-fly weight determination. *J. Chem. Phys.* **2013**, *138*, 061102.

(86) Nguyen, H.; Roe, D. R.; Simmerling, C. Improved Generalized Born Solvent Model Parameters for Protein Simulations. *J. Chem. Theory Comput.* **2013**, *9*, 2020–2034.

(87) 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., Bioinf.* **2004**, *55*, 383–394.

(88) McGibbon, R. T.; Beauchamp, K. A.; Harrigan, M. P.; Klein, C.; Swails, J. M.; Hernández, C. X.; Schwantes, C. R.; Wang, L.-P.; Lane, T. J.; Pande, V. S. MDTraj: A Modern Open Library for the Analysis of Molecular Dynamics Trajectories. *Biophys. J.* **2015**, *109*, 1528–1532.

(89) Scherer, M. K.; Trendelkamp-Schroer, B.; Paul, F.; Pérez-Hernández, G.; Hoffmann, M.; Plattner, N.; Wehmeyer, C.; Prinz, J.-H.; Noé, F. PyEMMA 2: A Software Package for Estimation, Validation, and Analysis of Markov Models. *J. Chem. Theory Comput.* **2015**, *11*, 5525–5542.

(90) Lin, J. Divergence measures based on the Shannon entropy. *IEEE Trans. Inf. Theory* **1991**, *37*, 145–151.

(91) Naughton, F. B.; Alibay, I.; Barnoud, J.; Barreto-Ojeda, E.; Beckstein, O.; Bouysset, C.; Cohen, O.; Gowers, R. J.; MacDermott-Opeskin, H.; Matta, M.; Melo, M. N.; Reddy, T.; Wang, L.; Zhuang, Y. MDAnalysis 2.0 and beyond: fast and interoperable, community driven simulation analysis. *Biophys. J.* **2022**, *121*, 272a–273a.

(92) MacQueen, J. others Some methods for classification and analysis of multivariate observations. In *Proceedings of the fifth Berkeley symposium on mathematical statistics and probability*, 1967, pp 281–297.

(93) McInnes, L.; Healy, J.; Astels, S. hdbscan: Hierarchical density based clustering. *J. Open Source Softw.* **2017**, *2*, 205.

(94) Anandakrishnan, R.; Drozdetski, A.; Walker, R. C.; Onufriev, A. V. Speed of Conformational Change: Comparing Explicit and Implicit Solvent Molecular Dynamics Simulations. *Biophys. J.* **2015**, *108*, 1153–1164.

(95) Mulligan, V. K.; Kang, C. S.; Sawaya, M. R.; Rettie, S.; Li, X.; Antselovich, I.; Craven, T. W.; Watkins, A. M.; Labonte, J. W.; DiMaio, F.; Yeates, T. O.; Baker, D. Computational design of mixed chirality peptide macrocycles with internal symmetry. *Protein Sci.* **2020**, *29*, 2433–2445.