An Eﬃcient GaMD Multi-Level Enhanced Sampling Strategy for Polarizable Force Fields Simulations of Large Biological Systems

We introduce a novel multi-level enhanced sampling strategy grounded on Gaussian accelerated Molecular Dynamics (GaMD). First, we propose a GaMD multi-GPUs-accelerated implementation within Tinker-HP. For the speciﬁc use with the ﬂexible AMOEBA polarizable force ﬁeld (PFF), we introduce the new "dual–water" GaMD mode. By adding harmonic boosts to the water stretching and bonding terms, it accelerates the solvent-solute interactions while enabling speedups with fast multiple–timestep integrators. To further reduce time-to-solution, we couple GaMD to Umbrella Sampling (US). The GaMD—US/dual–water approach is tested on the 1D Potential of Mean Force (PMF) of the CD2–CD58 system (168000 atoms) allowing the AMOEBA PMF to converge within 1 kcal/mol of the experimental value. Finally, Adaptive Sampling (AS) is added enabling AS–GaMD capabilities but also the introduction of the new Adaptive Sampling–US–GaMD (ASUS–GaMD) scheme. The highly parallel ASUS– GaMD setup decreases time to convergence by respectively 10 and 20 compared to GaMD–US and US.

Understanding interactions within biomolecules is crucial for many topics such as drug discovery.Some structural modifications, sometime undetected by experiment, can drastically change the nature of the physics ruling interacting complex systems.2][3][4][5][6] It requires accurate models able to capture the true potential energy hyper-surface and long simulations to both access the large biological processes time-scale and satisfy the ergodicity principle. 7][10][11] Beside these developments, several additional strategies have been pursued overs the years to further speed up the simulations.[15] Alternatively, an intensive algorithmic work has been undertaken, introducing techniques such as multiple-time-step integrator schemes 16,17 or collective variables-driven molecular dynamics methods. 18,191][22][23][24][25] Although such methods are powerful as they can estimate free energies of binding or the stability of secondary and quaternary structures of proteins, 26,27 the free energy estimations can suffer from biases either generated by the initial choice of the collective variable (CV) or by the existence of multiple CV within the mechanism process (e.g dual mechanisms). 28For these reasons, collective variable-free methods have become increasingly popular. 29Among them, the recent Gaussian accelerated Molecular Dynamics (GaMD) has shown great promises due to its high sampling acceleration, its user-friendly tunable parameters and its minor additional computational cost. 30GaMD accelerates conformational sampling by adding a harmonic boost to the potential energy.Coupled with the second order cumulant expansion, GaMD allows us to compute unbiased properties by using an accurate reweighting procedure through cumulant expansion to the second order.
Although new generation many-body polarizable force fields (PFFs) are more accurate in describing biomolecular interactions, [31][32][33] they are computationally more challenging than traditional approaches.Therefore, to overcome these limitations, we provide here a novel multi-level enhanced sampling strategy to accelerate PFFs simulations with the AMOEBA force field.To do so, we combine together the Tinker-HP massively parallel multi-GPUs platform 15 to a highly scalable GaMD implementation (level 0 ) and to additional enhanced sampling techniques based on recent developments of the field.As a first speedup, we propose an extension of the GaMD formalism with a new GaMD mode enabling the use of the AMOEBA PFF 34,35 and fast multiple-time-step integrators 17 (level 1 ).We then discuss the explicit coupling of such PFF-oriented GaMD approach to Umbrella Sampling (US) 36 and Adaptive Sampling (AS) 6 techniques (level 2 ).To demonstrate their applicability to PFF, these physics-based hybrid enhanced sampling strategies are then applied to the Potential of Mean Force (PMF) study of a large biological complex CD2-CD58 interacting via salt bridges.Finally, we combine all together within the Adaptive Sampling-US-GaMD method (ASUS-GaMD) scheme (level 3 ).
Introducing the GaMD "dual water" mode.GaMD is a potential-biasing method for unconstrained enhanced sampling without the need to set predefined CV.It smooths the potential energy surface by adding a harmonic boost potential as described in the seminal paper 11 .Its general framework makes it suitable for the development of hybrid schemes and variants, such as replica-exchange umbrella sampling GaMD (GaREUS), 37 Ligand GaMD (LiGaMD) 38 and Peptide GaMD (Pep-GaMD). 39 the system potential energy is lower than a threshold energy E, a harmonic potential energy boost is applied to smooth the potential energy surface.By denoting q ∈ R 3N the configurations, when the system potential energy U (q) is lower than a threshold energy E, a boost, which depends on U (q) is added: with ∆U GaM D (q) the external harmonic potential boost and k the harmonic force constant.The two adjustable GaMD parameters k and E are automatically determined following the original procedure described in ref 30 .The boost intensity can be managed thorough a user-specified uper limit labeled as σ 0 (e.g 10k B T ) predefined before the simulation.To ensure accurate reweighting with the cumulant expansion the ∆U GaM D standard deviation, σ ∆V , should satisfy σ ∆V < σ 0 30,40,41 .GaMD provides different modes: the boost is either applied on the total potential GaMD-pot, on the dihedral potential GaMD-dih, or on both at the same time GaMD-dual 42,43 .Recently, another mode was introduced:LiGaMD which adds the boost to a ligand non bonded interactions, 38 accelerating the sampling of ligand-protein interactions.It is known that interactions involving water are essential for such systems and that protein stability processes are controlled by water-protein interactions. 6,44,45To accelerate these interactions, one would like to use the GaMD-dual mode on the non bonded interactions of water molecules.However, such boost requires the evaluation of the complete non bonded energies and, in the context of multi-timestep integrators such as BAOAB-RESPA1 17 where they are split between short and long range, these are only available at outer (large) timestep.But the fluctuations of the associated bias are such that it has to be evaluated at shorter timesteps, so that the whole procedure is not compatible with multi-timestep integrators such as BAOAB-RESPA1.For similar reasons, the GaMD-dual mode with a bias applied to the complete potential energy is not compatible with even simple RESPA integrators in which the potential energy is split between bonded and non bonded terms.Therefore, GaMD-dual mode becomes rapidly limited by the simulation time.To overcome this issue, we developed a new mode, GaMDdualwater (denoted GaMD-dualwat), which adds a boost to the protein dihedral potential energy term and the water stretching and bending terms, this time fully compatible with RESPA and RESPA1 like integrators, allowing water molecule to be more flexible and thus favoring their conformational changes.
∆U GaM D−dw (U (q)) = ∆U dihedral protein (U (q)) + ∆U stretch water (U (q)) + ∆U bend water (U (q)) This mode is enabled by the flexibility of the AMOEBA 03 water model 34 but is not compatible with the TIP3P rigid water model used in most of non-polarizable force fields such as CHARMM and AMBER. 46This framework allows to further reduce the computational cost gap between PFFs and nPFFs.This new mode, in addition to the other GaMD-dih and GaMD-dual modes, is now available within the Tinker-HP software. 12,15In the following, we first tested its GPU scalability and performance on the STMV system ( 1 066 624 atoms), and its sampling efficiency is demonstrated on simulations of the alanine dipeptide and the CD2-CD58 complex.
Level 0: Efficiency and GPU scalability.The GaMD implementation is such that only a small computational and communication (in parallel) overhead is added compared to cMD.
The GaMD-dih and GaMD-dualwat have been considered on the STMV system (1 066 624 atoms) with the AMOEBA PFF and the 10 fs outer time-step HMR BAOAB-RESPA1 multiple-time-step integrator. 17V100 GPUs from the national Jean Zay supercalculator have been used for all the benchmark computations.Similar scalability studies have been performed on the Jean Zay CPUs (Figure S1).On 1 and 2 GPUs (Figure 1), the GaMD data communications are negligible, 1%.On 4 GPUs, the communications are increasing and the performance decreases by 7%.Overall, the use of GaMD only slightly alters the performance.This high scalability opens the door to simulate at a high-accuracy large complex biomolecular systems with PFFs.We used the many-body AMOEBABIO18 PFF. 47,48The system was minimized with a RMS of 1 kcal/mol and sampled within the NPT thermodynamic ensemble with the Bussi thermostat 49 and a MonteCarlo barostat 50 at 300 K and 1 atmosphere.We used the Velocity Verlet integrator and a 1 fs time-step. 51Smooth Particle Mesh Ewald (SPME) algorithm was employed to compute non-covalent interactions 52 with a real space cutoff equal to 7 Å and a Van der Waals cutoff set to 9 Å.For AMOEBA the convergence criteria for multipoles was set to 10 −5 .After short testing simulations, we found an optimal value of 3 kcal/mol for GaMD-dih and GaMD-dual σ 0 , in accordance with ref 11 , and 4 kcal/mol for GaMD-dualwat (see Figure S2, Table S1 and S2).We ran 3 independent simulations of 60 ns for each mode.The different sampled basins are also compared to a 1 µs cMD AMOEBA reference.
Reweighted, see Technical Appendix, free energy surfaces obtained from these simulations are depicted on Figure 2 and show that GaMD-dual captures well the α r (50°,25°), α L (-75°,-25°) and P II (-75°,150°) basins.These results are consistent with the 1 µs cMD trajectory (see Figure S3 in SI) depicting these three basins.While the GaMD-dih mode capture the α r basin after 150 ns, the GaMD-dualwat capture it in 100 ns (SI Figure S4).We also observe a sampling acceleration between GaMD-dual and GaMD-dualwat compared to the reference 1 microsecond cMD.To characterize the GaMD boost harmonicity, its distribution anharmonicity γ is calculated as in. 30γ serves as an indicator of the sampling convergence and reweighting procedure accuracy.Depicted on Figure S5 in SI GaMD-dih as well as GaMD-dual depicts high anharmonicity with respectively 0.252 and 0.016 compared to GaMD-dualwat with 0.005.Additionally, we see a steep anharmonicity convergence to less than 10 −3 for GaMD-dualwat while being relatively stable at 2 × 10 −1 for GaMD-dih (SI Figure S4).In comparison the anharmonicity is about 0.001 with GaMD-dih and AM-BER99SB.PFFs thus increase the statistical noise and stress the importance of using low anharmonicity GaMD modes.In that sense, GaMD-dualwat appears more suitable than GaMD-dual for PFFs simulations with an anharmonicity equal to 0.0005.As stated before, another advantage of GaMD-dualwat is that it can be coupled to multiple-time-step such as BAOAB-RESPA1 17 in contrast to the GaMD-dual mode that remains limited to single timestep integrators.Comparative results of GaMD-dualwat with both integrators can be found in Figure S6 of the SI.Its coupling with multiple-time-step clearly compensates the slightly lower sampling performance compared with GaMD-dual.The sampling enhance-ment brought by the GaMD-dualwat can be partly related to how it affects the diffusion of water: we report in Table S3 of the SI the self diffusion coefficients of bulk water computed within a same setup (same size of box and same integrator) and observe that it is increased with the GaMD-dualwat mode compared to the simple GaMD-dih one, favoring global conformationnal changes due to water reorganization.While the added sampling efficiency is already significant for the alanine dipeptide, we expect it to be larger on more complex and larger biological systems such as CD2CD58 where water reorganization plays a bigger role.
Combined with a highly-parallel GPUs infrastructures and multiple-time-step integrators, the GaMD-dualwat should allow to help reaching very high-resolution conformational space of large molecular systems.In addition to the sampling acceleration it provides, the low associated anharmonicity drastically reduces the statistical noise associated with reweighting.
Level 2: Speeding-up simulations with the parallel AS-GaMD scheme.We further coupled our newly introduced GaMD mode to additional enhanced sampling strategies.Recently, we developed a new adaptive sampling technique (AS) which was shown to allow massive sampling of the SARS-CoV-2 Main Protease conformational space. 53We coupled these two methodologies together, yielding the AS-GaMD method.The principle is similar to the AS, the only modification being that each cMD at each iteration is now a GaMD simulation.The double bias coming from both AS and GaMD implies that a suitable and careful reweigthing scheme has to be introduced to reconstruct unbiased free energy surface.All mathematical tools for the reweigthing scheme are provided in the Technical Appendix.We applied this methodology to the same system, the alanine dipeptide using the same simulation protocol.
At each iteration, we projected the structures on the two main dihedral angles space.The probability law for selection of new structures was taken as the inverse of the square of the probability density on this reduce space.This choice aims at improving the sampling of undiscovered region of the given space.
Phase space sampling obtained from AS-GaMD simulation after both GaMD and AS reweigthing are depicted on Figure 3. 5 iterations of AS-GaMD with 5 GaMD simulations of 5 ns have been performed for a total simulation time of 110 ns.All the known region of space captured by the GaMD simulations are captured by AS-GaMD but with more spread basins, meaning that the GaMD-dualwat sampling has been enhanced by its coupling with AS.Furthermore, we also observed that the α r regions are rapidly captured during the first iteration, i.e with only 25 ns (SI Figure S7).It thus represents an important gain in simulation time compared to the GaMD-dualwat alone, and is thus promising to sample conformational space of larger and more complex biomolecular systems.
Level 3: Pushing the limit of PMF convergence with GaMD-US and ASUS-GaMD.
5][56] In addition to the choice of the CVs, it is also difficult to estimate the PMF convergence as it is system dependent.Good indicators to check if convergence is reached are: the overlap between neighboring windows and the evolution of the PMF curve as a function of the simulation time per window.To accelerate the sampling within each window, Oshima et al. recently combined GaMD with replica-exchange and US. 37Here, we first only applied a GaMD boost in each US window in order to enhance the sampling in the orthogonal space.
To demonstrate the PMF convergence acceleration, we studied the dissociation of the salt bridges interface within the CD2CD58 complex.This system, made of several salt bridges and hydrogen bondings interactions, was already studied by some of us. 28Although it has been shown that PFFs allow a better description of the salt bridges interactions, their computational cost has long hindered the study of such large system.Since the portability of Tinker-HP on multi-GPU and the global acceleration of the PFFs, reaching such system is now easily achievable.To start this study we took the same CD2CD58 complex as in our previous study (e.g 1226 atoms) but we solvated it in a waterbox of 100 × 100 × 100 Å.
Counterion were added to neutralize the system.We used the AMOEBABIO18 PFF. 47,48e system was minimized with a RMS of 1 kcal/mol in the NVT thermodynamic ensemble with the Bussi thermostat. 49Temperature was set to 300 K while pressure was set to 1 atmosphere.We used the multiple-time-step BAOAB-RESPA1 with a 10 fs timestep with the Hydrogen Mass Repartionning scheme (HMR) 17 and Smooth Particle Mesh Ewald (SPME) algorithm to compute electrostatic and polarization interactions 52 with a real space cutoff of 7 Å and a Van der Waals cutoff of 9 Å.The convergence criteria for polarization was set to 10 −5 .39 US windows were generated, ranging from 1 to 20 Å with a width of 0.5 Å between them.CV was chosen as the distance between the center of mass formed by the interfacial residues isolated by Bayas et al. on CD2 and CD58 (Table 1 in ref 57 ).A spring constant of 10 kcal/mol.Å 2 was employed to restrain the system along the chosen CV.Each window was run for 5 ns for equilibration and then for 50 ns.Histogram overlap as well as the PMF curve as a function of the simulation time allocated per window were employed to check the convergence of the simulations (SI Figure S9).The final US PMF show a slow decrease of the free energy barrier with the simulation time, suggesting a slow convergence to about 12.5 kcal/mol.Binding affinity was found to be experimentally around 7.1 ± 0.03 kcal/mol, suggesting that our simulations are not converged. 57In order to improve sampling within each window a new US was performed similar to the previous US protocol but now with an additional GaMD-dualwat potential applied in each window.The GaMD parametrization protocol and reweighting procedure are described in the Technical Appendix and in SI (Figure S8 and Table S4).The optimized GaMD-dualwat parameters σ 0 are equal to 1 and 3 kcal/mol for respectively dihedral and dual water modes.Figure 4 shows the difference between standard US and GaMD-US.The GaMD-US PMF and boost harmonicity converge at 40 ns per window (SI Figure S10 A/, B/ and C/).The predicted free energy barrier is now within the 1 kcal/mol of the experiment.It shows that GaMD-dualwat, even without presence of Replica Exchange, could considerably improve PMF convergence of large systems.
It also demonstrates that salt bridges and, more generally, protein-protein interactions are well described with PFFs while it was demonstrated that non-PFFs fail to describe these interactions. 58 further push the sampling, we coupled together GaMD, AS and US (ASUS-GaMD).We provide two reweighting schemes that either use modified Multistate Bennett Acceptance Ratio (MBAR) equations or the Rao-Blackwell estimator.The mathematical expressions are general and can be used with any weighted dynamics.Starting from an initial US simulation ( few ns), each window is decomposed in several AS independent trajectories with an additional GaMD-dualwat potential boost (GaMD).Here, we ran 2 iterations of 5×5 ns GaMD-US per window.The PMF evolution can be found in SI (Figure S11) while the resulting PMF is depicted on Figure 4. We observe that ASUS-GaMD reach GaMD-US in one iteration, showing the sampling acceleration impact provided by the AS part within ASUS-GaMD.Although a careful reweighting is needed for the different AS, GaMD and US layers, the overall ASUS-GaMD approach inherits the strong adaptive sampling advantages of being pleasantly parallelizable and considerably accelerates the PMF convergence.
Combined with the use of modern GPUs, these sampling techniques allow to crush time to solution in PFFs evalution of PMFs.Although it is difficult to truly quantify the final speedup (i.e. a PMF convergence remains partially system-dependant), one can see in Figure S12(SI) that if we extrapolate the US convergence, ASUS-GaMD converges 1.4 times faster.
Thanks to the native parallelism inherited from AS, the PMF evaluation can be done in one fifth of the simulation time yielding a speedup of 7. If we consider that convergence was already reached with a 25ns per window setup, this factor grows to 14. ASUS-GaMD can thus reduce to days computation that would have taken months.This work also allows to invoke any variant of the combined approaches, offering therefore access to GPU-accelerated GaMD-adaptive sampling (AS-GaMD) simulations that will be helpful to further extend conformational space studies of proteins 6 To conclude, these methodologies will contribute further to allow high-resolution sampling of large biological systems up to millions of atoms using polarizable force fields.

Technical Appendix
We denote by ξ(q) the reaction coordinate along which we performed the US simulation and q the configuration.Here, a configuration means the positions q ∈ R 3N of all the atoms of the system.The imposed US bias potential is with K the force constant.
We combined the AS, US and GaMD such that each US window j ∈ [[1, ..., M ]], ξ 1 , ..., ξ M , is parallelized and accelerated by adaptive sampling replicas and GaMD boost potential: We denote by (q j,n ) n∈1,N the N configurations generated by the AS replicas of US window j and (ω j,n ) n∈1,N their respective AS weights.These weights are normalized such that with v j,n the unnormalized AS weights.The canonical average of an observable ϕ is estimated by ϕ j = ϕ(q)e −βU j (q) dq e −βU j (q) dq In practice, to get a smooth reweighted PMF, the reaction coordinate ξ is discretized in K bins around values x 1 , ..., x K .We want to estimate for each k ∈ [[1, ..., K]] its free energy, up to an additive constant, As q is distributed according to the density probability law e −βU e −βU , with ϕ k = 1 ξ(q)∈Bin(x k ) .
1st step: GaMD with cumulant expansion We, first, remove the GaMD bias.Here, we want to find a relation between ϕ and ϕ where the prime average represents the canonical average over the potential U = U +U GaM D .
Starting from the canonical average, we notice : ϕ = ϕ(q)e −βU (q) dq e −βU (q) dq = ϕ(q)e βU GaM D (q) e −βU (q) dq e βU GaM D (q) e −βU (q) dq = ϕe βU GaM D e βU GaM D By applying this with ϕ = ϕ k , where C is a constant and F (x k ) is the free energy F (x k ) = − 1 β ln ϕ k .To reduce the estimator variance, we used the cumulant expansion to the second order, By combining with equation 10, the free energy is rewritten as 2nd step: AS modified MBAR

Modified MBAR
Let's define c j = e −βU j (q) dq, The prime comes from the use of the MBAR on the reference energy U of the previous section.The starting point is to use the MBAR identity (ref 59 equation 5) and notice c i e βU j α i,j i = e −βU j (q) e −βU i (q) α i,j (q) dq = c j e βU i α i,j j which holds for arbitrary functions q −→ α ij (q) with i, j ∈ [[1, ..., M ]].Notice that each window generated the same number of configurations N .The MBAR estimator has been proven to be optimal by using and by summing over j M k=1 e −βU k (q) /c k j (17)   Using equation 6 we obtain the estimators ω j,n e −βU i (q j,n ) M k=1 e −βU k (q j,n ) /ĉ k (18)  and finally with eq 13 ω j,n e −βU i (q j,n ) M k=1 e β(F k −U k (q j,n )) (19)   which must be solve self-consistently.
The RB theorem characterizes the transformation of a crude estimator into a better estimator that has smaller mean-squared-error w.r.t to the dataset.Assuming M equilibrium states sampled independently, with potential U i .The biased free energy F i is equal to Where N = M i=1 N i and b i are unknown biased energies added to state i to adjust the relative weight to match the free energy.b i were introduced to make equation 20 valid which is the requirement for using the RB estimator.The RB estimator for this ensemble is ω j,n e −β(U i (q j,n )+b i ) M k=1 e −β(U k (q j,n )+b i ) Combining with equation 20: ω j,n e −β(U i (q j,n )+b i ) M k=1 e −β(U k (q j,n )+b i ) Thus, the unbiased free energy F * i can be calculated using equation 20 after solving 22 for b i .Equation 22 has major interests: (1)it is more stable, (2)reduces the number of floating point operations and (3) the problem is reduced to minimizing a convex function.Indeed, if with r i,n the weight of configuration q(i, n) r i,n = M ω i,n M k=1 e −βU k (q i,n ) /ĉ k (27)   c 0 is unknown but does not depends of k ∈ [[1, ..., K]] so equation 12 can be rewritten as with

Figure 3 :
Figure 3: AS-GaMD sampling obtained on the alanine dipeptide after reweithing the 5 ( of AS iterations) × 5 ( of independant GaMD simulations per iteration) × 5 (simulation time of one GaMD simulation) ns of simulation.

Figure 4 :
Figure 4: PDB 1QA9 CD2CD58 representation with CD2 and CD58 subcomplexes represented respectively in blue and red using the newribbons representation.Residues at the interface considered in the COM distance between the two subcomplexes are represented in blue and red for respective basic and acid residues using the CPK representation.VMD software was employed to generate the structure.

N
n=1 ω j,n = N and ω j,n = v j,n N m=1 v j,m

Finally, we want
to express ϕ w.r.t the AS weights in each US window j ∈ [[1, ..., M ]].This can be done in two ways either using the MBAR or the Rao-Blackwell estimator.

e
−βU j /c j M k=1 e −βU k /c k i = M j=1 c j e −βU i /c j M k=1 e −βU k /c k j(16)We obtain a set of M equations for all i ∈ [[1, ..., M ]]c i = M j=1e −βU i (q)