{'fraction_non_alphanumeric': 0.026535131870344272, 'fraction_numerical': 0.019367827662572983, 'mean_word_length': 4.763750290090508, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Combining Enhanced Sampling with Experiment Directed Simulation of the GYG peptide

Experiment directed simulation is a technique to minimally bias molecular dynamics simulations to match experimentally observed results. The method improves accuracy but does not address the sampling problem of molecular dynamics simulations of large systems. This work combines experiment directed simulation with both the parallel-tempering and parallel-tempering well-tempered ensemble replica-exchange methods to enhance sampling of experiment directed simulations. These methods are demonstrated on the GYG tripeptide in explicit water. The collective variables biased by experiment directed simulation are chemical shifts, where the set-points are determined by NMR experiments. The results show that it is possible to enhance sampling with either parallel-tempering and parallel-tempering well-tempered ensemble in the experiment directed simulation method. This combination of methods provides a novel approach for both accurately and exhaustively simulating biological systems.


Introduction
Molecular simulations are connected to experiments through ensemble properties, called observables. Often experimental and simulated values disagree, due to challenges such as limited accuracy of force-fields and discrepancies between how experiments can be represented in a simulation. This disagreement can be improved by adding an auxiliary biasing potential energy function. 1,2 For example, a simulation can be improved by incorporating experimental results as structural restraints in simulations. 3 Often these restraints are a type of harmonic energy penalty that brings conformational ensembles closer to the experimentally observed values. 4,5 Visiting states that are outside structural restraints is unfavored with harmonic biases. This introduces some ambiguity in how strong these energetic restraints should be.
To address the ambiguity in using energetic restraints to match experimental data, a maximum entropy principle 6 can be applied to find an energy bias that is minimal. 7,8 This bias acts only in the dimension of the collective variable which corresponds to the experimental observable. 9 Experiment directed simulation (EDS) is an implementation of maximum entropy biasing. EDS uses an adaptive biasing method to converge to the maximum entropy bias that achieves agreement between molecular simulation collective variables and experimental observables. 10 This bias is also provably the smallest change that can be made to match experimental data. 8,10,11 EDS utilizes a single replica, whereas most other maximum entropy methods utilize multiple replicas. 1 Other biasing methods, such as replica-exchange chemical shift restrained molecular dynamics simulations 4 and restrained-ensemble molecular dynamics simulation method 12 that utilize multiple replicas require multiple replicas in order to match simulated results to some reference values. This increases the required computational effort, although there is an improvement in sampling due to the inherent increase of samples with the replicas. EDS, being a single-replica technique, provides no improvement to sampling. EDS has been applied to classical molecular dynamics, 10 ab-initio molecular dynamics, 13 and coarse-graining. 14 EDS was able to give more accurate results on dynamic properties such as the self-diffusion coefficient of lithium in electrolyte solutions, 10 and the proton-diffusivity in ab initio water 13 as well.
EDS is a biasing method and this is designed to match the ensemble average of some collective variable to a certain reference value. EDS lets the underlying simulation determine the fluctuations and underlying distribution of the collective variable. Only the expected value is shifted to match the experimental value. Because EDS is only a method of simulation, it cannot correct or mitigate any error in the reference value, such as experimental error, human error, or instrumentation error. If the experimental value is instead a distribution or becomes a distribution because a Bayesian inference approach with prior belief is desired, other techniques allow matching to that desired distribution like experiment directed metadynamics 11 or metainference. 15 Aside from accuracy, another challenge in molecular dynamics simulations is sufficient sampling of phase space. The so-called "sampling problem" 16 is especially true in solvated macromolecules, where systems can become trapped in energy minima. EDS introduces a bias in the potential to match experimental averages; it is not an enhanced sampling technique. Thus, there is a need to use an enhanced sampling techniques with EDS. Multiple enhanced sampling methods have been developed to alleviate this problem, 17 often by biasing sampling to occur along a certain collective variable. 18 This presents a complication with EDS if both an enhanced sampling technique and EDS are applied on the same set of collective variables. For example, metadynamics 19 will add Gaussian bias to the potential energy with respect to a collective variable, thereby modifying the potential of mean force (PMF) of that CV. EDS cannot be applied to a collective variable with a biased PMF, because EDS requires canonical samples. For example, if metadynamics is applied to Φ dihedral angle and EDS is applied Ψ dihedral angle, the Ψ angle needs to be re-weighted to achieve canonical sampling due to bias on Φ and the inherent correlation between the Φ and Ψ. EDS needs the Ψ variable to be canonically sampled. Hence, EDS and Metadynamics cannot be applied simultaneously in a simulation. The well-tempered ensemble method (WTE) provides an alternative approach because it biases the fluctuations in energy, and with an assumption of the underlying potential energy distribution, does not affect the canonical sampling of other CVs. WTE can be further combined with parallel-tempering replica-exchange to improve sampling further. [20][21][22] In this work we consider the parallel-tempering well-tempered ensemble (PT-WTE) combined with EDS 22 to improve both sampling and accuracy on an example system. Like in parallel-tempering replica-exchange, N replicas of the simulation system are run concurrently. The coldest replica is the system of interest, yet it can escape free energy minima by swapping configurations with higher temperatures. The goal of PT-WTE is to reduce the number of replicas while maintaining similar replica efficiency. 22 It has been shown that PT-WTE metadynamics decreases the number of replicas from 100 to 10 while maintaining the similar energy landscape of solvated tryptophan-cage protein. 23 It has also been shown to improve speed of convergence to an accurate FES in model and real systems. 22,23 The technique can be applied to larger systems such as proteins to efficiently sample their large conformational space. 24 The tri-peptide glycine-tyrosine-glycine (GYG) was used as a model system to demonstrate an implementation. GYG is found in structural and transport proteins such as collagen 25 and ion channels. 26 Either mutations or complete deletions of this sequence from proteins showed either complete or partial loss of a function in those proteins, 27-29 so it has specific functional relevance. The experimental data utilized here is from backbone chemical shifts 30 of GYG peptide. These are used to assess accuracy of the simulation and the free energy surface along the Φ-Ψ dihedral angles is used to assess sampling.
Metadynamics metainference 31 and replica-averaged metadynamics 32 are the nearest examples of the method proposed here in literature. Metadynamics metainferences is a combination of the enhanced sampling technique of metadynamics 19 with the ability to match experimental data of metainference. 31 Replica-averaged metadynamics is the combination of the replica-averaging approach to matching experimental data 32 with the enhanced sampling of metadynamics. A distinguishing difference between our approach and these is that the biasing force from EDS can be calculated with an enhanced-sampling method first and then be applied to a second simulation without replicas. This capability allows dynamic properties like autocorrelation functions or diffusivities to be computed.

Theory
In the WTE method, 22 Newton's equation of motion becomes: where m is the mass vector, R is the position vector, U (R) is the system potential energy and V (·) is a timedependent function of the instantaneous potential energy which biases the system. V (·) is defined with the following time-derivative:V where ω and ∆T are tuned constants. The intuition of this equation is to add a Guassian at each time t centered at U (t), just like in the metadynamics method. 19,33 The key parameter in WTE is the bias factor, like in well-tempered metadynamics, 33 which is defined relative to system temperature and ∆T as γ = (T + ∆T )/T Simulating according to this bias will asymptotically lead to the following distribution of potential energy, assuming normality of the potential energy probability distribution: 22 where U and ∆U 2 are the expected value and variance of the potential energy of the unbiased canonical sampling. This expression shows that under the WTE bias, the average energies are the same but the variance of distribution is increased by √ γ. A bias factor of 1 thus leads to canonical sampling. Increasing the bias factor increase the fluctuations in potential energy. EDS modifies the potential energy of a system with a maximum entropy linear bias in a collective variable. 9 After equilibration, EDS will have the following biased potential energy: where α is a Lagrange multiplier determined from an equilibration process described in White et al. 10 and f (R) is a differentiable collective variable. The α is selected such that f =f , wheref is a pre-determined value for what the expected value of the collective variable should be. In this work, we have experimentally determined NMR chemical shifts forf . The limitation of EDS is that finding α requires calculating f with good sampling of the underlying canonical distribution. The PT-WTE method can ensure we have sufficient sampling of the distribution of f values. The WTE process only needs to be modified to replace U (R) with U ′ (R) and the previous analysis holds.
The benefit of WTE is that the increased energy fluctuations arising from the choice of bias factor can improve replica-exchange efficiency in a parallel-tempering setting. The probability of exchange for replicas i and j becomes: 34 where β i = 1/k B T i and U ′ i is the potential energy of the ith replica including the EDS bias. Note that EDS will have different α for each replica because f is a function of temperature. This exchange rule ensures that the resulting simulation will arrive at a canonical distribution without re-weighting. The result of this method is both an enhanced-sampled ensemble and a value for α.

Methods
A short peptide GYG was simulated with EDS 10 and PT-WTE 22 to enhance sampling and improve fit to experimental NMR chemical shifts from Platzer et al. 30 To compare the performance of EDS with enhanced sampling to other simulations, EDS with PT-WTE was compared to a single replica unbiased simulation, a single replica biased simulation, a parallel-tempering replica-exchange molecular dynamics (PT) simulation, and replica unbiased PT-WTE simulation.
0.008 mg/cm 3 GYG peptide in 10 mM NaCl counter-ions was first energy minimized then annealed and equilibrated in the NVT ensemble using the Canonical Sampling through Velocity Rescaling (CSVR) thermostat. 35 The CHARM27 36 force field and TIP3P 37 water model were used throughout. Electrostatic forces were calculated with the particle mesh Ewald method, 38 and dispersion forces were calculated with shifted Van der Waals potentials with a cutoff distance of 10Å. The covalent hydrogen bonds were constrained using the LINCS algorithm 39 to enable a 2 fs time step. All simulations were performed in the Gromacs-5.1.4 simulation engine. 40 To automatically run Gromacs commands in python, GromacsWrapper was used. 41 To generate an initial extended structure of GYG, the python package PeptideBuilder was used. 42 An equilibration is required for the PT-WTE method. 16 replicas of PT-WTE were tuned for 400 ps to find parameters that gave increased replica efficiency. The replica temperatures used for the 16 replica systems were 293, 299, 304, 310, 316, 323, 329, 336, 342, 349, 355, 363, 370, 377, 385, 392 K, as generated by the method of Nemoto et al. 43 The exchanges were attempted every 250 time steps. The WTE parameters were chosen following Barducci et al. 33 and are shown in Table 1.
After the aforementioned short PT-WTE equilibration step, EDS was begun with the static bias from PT-WTE equilibration step. This combined PT-WTE with EDS had 16 replicas with the same set of temperatures as mentioned above and was 40 ns long. To control against the presence of either EDS bias, or PT-WTE bias, or both biases, 16 replica 40 ns PT simulations with and without EDS bias in the absence of WTE bias and 16 replica 40 ns PT-WTE simulations with and without EDS bias were run and as shown in Table 1. The large amount of time is not necessary for convergence of EDS, but was chosen to get good sampling statistics. The EDS collective variables biased were chemical shifts 30 calculated via the Plumed2 molecular simulation plugin. 44 The EDS parameters chosen are shown in Table 1.
An 8 replica PT-WTE simulation was also done with and without EDS. These simulations demonstrate how PT-WTE can use fewer replicas than PT and still achieve good efficiency. The parameters are in Table 1. The replicas ran at the following temperatures: 293, 306, 320, 335, 350, 366, 383, and 400 K. Similar to 16 replica PT-WTE, 8 replica PT-WTE initially consisted of 400 ps equilibration step, followed by a static production step of 40 ns. There were 3620 atoms in the simulation with 1197 tip3p water molecules and 38 peptide atoms. We used bias factor of 10 for 16 replica PT-WTE and 20 for 8 replica PT-WTE. Deighen et al. studied the effect of bias factor in 10 replica simulation that studied 39 amino acids long Trp-cage protein. 23 They found that increasing the bias factor from 10 to 24 decreased the the time it takes to converge to the reference FES. 23 Hence, we used larger bias factor for 16 replica PT-WTE compared to 8 replica PT-WTE.
Overview of simulation systems and their parameters.  Table 1 provides an overview of the simulations run for this work. Multiple combinations of enhanced sampling methods with and without EDS were simulated. Each simulation ran for 40 ns. The 16 replica PT-WTE simulations were equilibrated for 400 ps to find the biasing potentials that increased the replica-exchange efficiency, prior to enabling EDS. This PT-WTE tune step was run long enough to ensure the efficiency exceeded 30%. Parallel-tempering replica-exchange (PT) with 16 replicas has 23% replica efficiency, as shown in Table 1.
The replica efficiency is defined here as the minimum acceptance rate of configuration swapping between replicas. In other words, the lowest among the 16 replicas is taken to be the efficiency. This ensures that all replicas exchange at certain threshold rate prohibiting a different independent system of simulations that are exchanging only with replicas that have closer temperatures. The averages are also shown for the 40ns simulations. Replicaexchange of 20% is the minimum for a well-sampled system, 40 and thus 16 was the minimum number of replicas for PT. The choice of PT-WTE simulation parameters achieved good exchange with the smaller, 8 replica system, despite spanning the same temperature range as the 16 replica systems (293K to 400K). Since 8 replica system had a different set of WTE parameters, EDS parameters were modified as well, which included the reduction of the speed of convergence of EDS (lower range parameter) and the increase the speed of convergence of the PT-WTE method (larger hill height). 22 The bias factor was correspondingly decreased to reduce the strength of bias more quickly. PT-WTE with 16 replicas had smaller replica efficiency than 8 replica PT-WTE in Table 1. This may be the result of a different set WTE and EDS parameters used in 8 replica and 16 replica simulations.
The EDS method moved the calculated backbone chemical shifts to the experimental values when applied as seen in Table 1. There is good agreement between the average chemical shifts in the EDS simulations and experiments with any of the enhanced sampling methods applied. Even with the slower parameters used in the 8 replica system, the EDS method still achieved the correct values. Figure 1 shows the cumulative calculated mean backbone chemical shifts over time during 16 replica PT-WTE with EDS bias (a), 16 replica PT-WTE without EDS (b), 16 replica PT with EDS (c), 16 replica PT without EDS, 8 replica PT-WTE with EDS (e), 8 replica PT-WTE without EDS (f), single replica EDS (g), and a standard MD (h). Chemical shifts were modified to show the absolute difference between cumulative averages and the corresponding experimental chemical shifts. Therefore, exact match between simulation and experiment is indicated by the convergence of cumulative calculated mean chemical shifts to 0. As expected, EDS biased average chemical shifts converge to the reference value with 16 replica PT-WTE (Fig. 1a), with 16 replica PT (Fig. 1c), and with 8 replica PT-WTE (Fig. 1e). EDS does not make instantaneous values of biased chemical shifts converge with the reference value, but makes the ensemble average to approach the reference value as shown in fig. S1 of the Supporting Information. In all of the EDS cases ( fig. 1a, 1c, 1e, and 1g), EDS converges within 1 ns. In all EDS bias simulations ( fig. 1a, 1c, 1e, and 1g), average chemical shifts of the hydrogen connected to backbone nitrogen atom had the highest deviation from predicted value compared to other atoms. The Camshift root mean square deviation for H connected to N is 0.56 ppm, whereas Camshift predicts deviations to be highest for N with RMSD of 3.01 ppm in a 28 protein test set. 45 During 40 ns simulation, the cumulative average backbone HN chemical shift was within 1 % of the reference chemical shift. Since PT-WTE and PT swaps configurations between different temperature replicas, chemical shifts also change whenever the conformations change as indicated by the greater fluctuations between mean chemical shifts and reference values in PT-WTE and PT with EDS (Fig. 1a, 1c, and 1d ) simulation compared to single replica EDS (Fig. 1g). In the absence of EDS bias, average chemical shifts in various simulations, shown in figures 1b, 1d, 1e, and 1f, do not converge to 1 because the simulation is less accurate. When EDS was absent in both 16 replica PT-WTE (Fig. 1b) and 16 replica PT (Fig. 1d), average chemical shifts fluctuated more drastically in Fig. 1b compared to simulation without PT-WTE as shown in Fig. 1d. This is due to the better exploration of phase-space with the PT-WTE method. The convergence under a variety of enhanced sampling methods shows that EDS still works well when combined with enhanced sampling.  To see the effect of enhanced sampling and the effect of biased chemical shifts to other collective variables in a simulation, Φ and Ψ dihedral angles were compared to regular unbiased MD and to Ramachandran plots from X-ray crystallography. 46 The simulated free energy landscapes of Tyrosine in GYG peptide along Φ and Ψ dihedral angles are displayed in Fig. 2. Dihedral angles were computed during the simulations and their histograms were used to compute the FES. EDS predicted the same global minimum regardless of presence ( Fig. 2a) or absence of PT-WTE (Fig. 2c) at Φ=-1.2 and Ψ=-0.9. The free energy surface shows different global minimum with EDS because the free energy distribution is biased to satisfy experimental conditions. As shown in figures 2a and 2c, EDS simulations show the global minimum to be at Φ=-1.2 and Ψ=-0.9 radians, different from the one observed by enhanced sampling alone (Fig. 2b). Enhanced sampling technique found a different global minimum at Φ=1 and Ψ=1, in the absence of EDS than in the presence of EDS. Also, PT-WTE without EDS enabled the exploration of a region at Φ=1 and Ψ in interval [-3.14,3.14] radians in figure 2b which is similar to Tyrosine Ramachandran map from Vitalini et al; 47 neither of EDS alone (Fig. 2c), EDS with PT-WTE (Fig. 2a) and regular MD (Fig. 2d) could sample around that region continuously. Although EDS in the presence of PT-WTE (Fig. 2a), did not show aberrant behavior from absence of PT-WTE in figure 2c, it showed improved sampling of global minimum observed in PT-WTE alone at Φ around 1 radians and at Ψ=1 radians compared to the absence of PT-WTE (Fig. 2c) and regular MD (Fig. 2d). Overall, PT-WTE improved sampling and EDS modified the FES to achieve agreement with the NMR chemical shifts.
Ramachandran plots from Ting et al. 46 of Tyrosine in the sequence GYG were compared to our results. The Ting et al. 46 results are derived from a bioinformatics analysis of X-ray crystallographic data of whole proteins. Tyrosine dihedral angle plot of Tyrosine in YG peptide from Ting et al., 46 was compared to our results, since YG was closest GYG. Ting 46 also reported the global minima to be on the negative Φ angle axis. Without EDS, enhanced sampling alone could not label the experimentally observed global minimum as global minimum at Φ=-1 and Ψ=2.8 in Fig. 2b. Ψ global minimum was around zero for both Ting et al. 46 and EDS results ( Fig. 2a and 2c). A local minimum at around Φ=1 and Ψ=1 radians reported by Ting et al., 46 was also observed as a minimum Fig. 2a, 2b, and 2c. This local minimum was not predicted by standard MD at all within 20 ns, as seen in Fig. 2d. EDS did not show drastic deviation in unbiased structural properties, such as a dihedral angle of Tyrosine. Indeed dihedral angles in EDS simulations with or without enhanced sampling showed similar results as X-ray crystallography did.
Another major question about disagreement between simulations and experiments is whether the underlying cause is lack of sampling or inaccuracy of the force-field/system. The 8 replica PT-WTE without EDS shows that even with good sampling, the simulation matches neither the chemical shifts nor the FES from the EDS simulation, as shown in Fig. 2b. This is not necessarily due to a deficiency in the CHARMM force field, but instead could be due to a difference in concentration, ions, and termini between the simulations and the work of Platzer et al., 30 from which the chemical shifts were measured.

Conclusion
Simultaneous enhanced sampling and experiment directed simulation (EDS) was demonstrated on the GYG peptide in explicit solvent. EDS improved the accuracy of the simulation by minimally biasing the backbone chemical shifts of the peptide to match experimental data from Platzer et al. 30 Single replica simulations with and without EDS bias ran for comparable time. Both parallel-tempering replica-exchange (PT) and paralleltempering well-tempered ensemble (PT-WTE) were used with EDS to improve sampling. Compared to absence of enhanced sampling in a control MD, system has significantly benefited from enhanced sampling. PT and PT-WTE provided improved sampling and did not interfere with EDS, with EDS converging within 1 ns. On the same computer architecture, 16 replica PT-WTE simulations and 16 replica PT with and without EDS completed after approximately the same period of wall-clock time. PT-WTE with 8 replicas spanning 293K to 400K was also demonstrated with EDS, showing that the PT-WTE method is able to reduce the number of replicas required relative to PT while not jeopardizing enhanced sampling. We have demonstrated that PT-WTE and PT can be applied with EDS simultaneously without interfering with each other in this system. We will exploit this new method on larger system of proteins and peptides, to take advantage of both enhanced sampling as well as biasing with experiments.
Supporting Information for Combining Enhanced Sampling with Experiment Directed Simulation of the GYG peptide Average of Chemical shifts is converging to 0, whereas moving window average of last 0.4 ns window is shown in less opaque shade of the same color. EDS does not try to match instantaneous values of chemical shifts to reference, but rather makes the cumulative average match the reference. Both moving averages and cumulative averages were subtracted by the corresponding backbone chemical shifts from experiments, which EDS is biasing to match. Perfect match is when cumulative average approaches 0.