Thermodynamic equilibrium of biological macromolecules under mechanical constraints
aa r X i v : . [ c ond - m a t . s o f t ] F e b T HER MODYNAMIC EQUILIBRIUM OF B IOLOGICALMAC ROMOLECULES UNDER MEC HANICAL C ONSTRAINTS
A P
REPRINT
Ashkan Shekaari
Department of PhysicsK. N. Toosi University of TechnologyTehran, Iran [email protected]
Mahmoud Jafari ∗ Department of PhysicsK. N. Toosi University of TechnologyTehran, Iran [email protected]
February 11, 2021 A BSTRACT
Equilibrating proteins and other biomacromolecules is cardinal for molecular dynamics simulationof such biological systems in which they perform free dynamics without any externally-applied me-chanical constraint, until thermodynamic equilibrium with the surrounding is attained. However, insome important cases, we have to equilibrate the system of interest in the constant presence of certainconstraints, being referred to as constrained equilibration in the present work. A clear illustration ofthis type is a single amyloid β -strand or RNA, when the reaction coordinate is defined as the distancebetween the two ends of the strand and we are interested in carrying out replica-exchange umbrellasampling to map the associated free energy profile as the dependent quantity of interest. In suchcases, each sample has to be equilibrated with the two ends fixed. Here, we introduced a simulationtrick to perform this so-called constrained equilibration using steered molecular dynamics. We thenapplied this method to equilibrate a single, stretched β -strand of an amyloid beta dodecamer fibrilwith fixed ends. Examining the associated curves of the total energy and the force exerted on thepractically-fixed SMD atom over the total timespan broadly supported the validity of this kind ofequilibration. K eywords Constrained equilibration · Biomacromolecules · Replica-exchange umbrella sampling
Examining the underlying free energy profile (FEP) is still of foremost importance in fully understanding a largenumber of chemical processes such as protein-ligand binding [1], or drug transportation across the cell membrane [2],which cannot be accurately predicted, for instance in the field of computer-aided rational drug design [3], withouthaving the associated free energy behaviors [4]. A rich set of numerical tools for reliable determination of freeenergy changes based on the fundamental principles of statistical mechanics [5] has so far been developed and is nowwithin reach. This so-called methodological set, in conjunction with recent advancements in computational power, hasunprecedentedly enhanced the fields of free energy calculations, and therefore computational structural biology [6],widening their areas of application as well.The theoretical framework for free energy calculations was established a long time ago and several approximationshave so far been accordingly developed [7, 8, 9, 10, 11, 12]. However, central to accurate determination of freeenergy difference between initial and final states of a system is to sufficiently explore the configurational space of theinitial state so that relevant, low-energy configurations of the final state could adequately be sampled. Conventionalmethods, particularly molecular dynamics (MD) [13] and Monte Carlo (MC) [14], are not indeed successful in thisrespect [15] based on the fact that the system of interest could be trapped in a few conformational states (i.e., local ∗ Corresponding author hermodynamic equilibrium of biological macromolecules under mechanical constraints
A P
REPRINT minima) during the simulation and the resulting potential of mean force (PMF) would accordingly be dependent onthe starting conformation of the system due to the so-called unrepresentative sampling.To circumvent the problem, since the late 1960s, a number of sophisticated techniques and strategies have so far beenproposed and advanced based on performing random walks using non-Boltzmann probability weight factors [16].Replica-exchange umbrella sampling (REUS) [17] is a nearly recent of them, which combines replica-exchange(RE) [18] and umbrella sampling (US) [19] techniques, and therefore, has both merits of RE and US. The formerbelongs to the so-called generalized-ensemble algorithms [20], while the latter is indeed a particular application ofthe more general importance sampling [21] used in statistics. In contrast to RE, all replicas in REUS have the sametemperature value and this is the potential energy that is exchanged.Based on the fact that US is a quasi-equilibrium method—in contrast to non-equilibrium approaches such as steeredmolecular dynamics (SMD) [22], or targeted molecular dynamics (TMD) [23, 24]—each replica in REUS has to besimulated and therefore equilibrated in the canonical ensemble independently. The problem arises when the reactioncoordinate is defined as the distance between the two ends of the system, such as a single RNA or amyloid β -strand,and the associated PMF is to be calculated along that collective variable. In such a case, we may have several replicas,each of which has to be equilibrated independently in a way that the related end-to-end distance (or equivalently thetwo ends of the strand) must be kept fixed during the equilibration simulation. To overcome this problem, we carriedout such a kind of constrained equilibration on a single, stretched β -strand of an amyloid beta dodecamer fibril usingthe constant-velocity protocol of SMD in that a differential value of the order of − Å (cid:14) fs was set to the pullingvelocity of the SMD atom, which, in turn, kept this atom practically fixed over the total timespan. Examining thetime-dependent total energies and forces acting on the SMD atom for all the replicas eventually verified the accuracyof the constrained equilibration. The computational setup was also described in Sec. 2 in detail. The initial atomic positions of a single amyloid β -strand was taken from the RCSB [25] PDB file of amyloid betaA β − dodecamer fibril with the entry code 2MXU [26]. The simulations were carried out on Debian style [27]Linux [28] systems supported by MPICH (version 2-1.4) [29, 30] using NAMD (version 2.14b2) [31] computer soft-ware with the July 2018 update of CHARMM36 force fields [32]. The VMD program (version 1.9.4a9) [33] was alsoused for post-processing. The fibril was solvated in a box of water under periodic boundary conditions with a unit-cellpadding of about 2 nm to decouple the periodic interactions. Minimization was then carried out for 50000 conjugategradient steps (100 ps) followed by a 2-ns free-dynamics equilibration at 310 K and 1.01325 bar in order to reach theequilibrium conformation of the fibril. Langevin forces with a Langevin damping constant of 2.5 (cid:14) ps along with theLangevin piston pressure control were applied for the NPT equilibration to keep fixed the temperature and pressure ofthe system. The integration timestep was also 1 fs.After equilibrating the whole fibril in water, a single β -strand was taken and extended by SMD simulations in vacuumin a way that the distance between the first atom (Glu :C α ; where C α is carbon, 11 is the residue number, and Gluis the associated amino acid) and the last one (Ala :C α ), defined as the reaction coordinate (Fig. 1), increased from31.4 Å (the initial value) to about 42.4 Å according to the function . i ( i = 0 , ..., . Therefore, 12 independentreplicas numbered consecutively from 0 to 11 were generated. Each replica equilibrated for 1 ns at 310 K in vacuumduring an SMD simulation in that Glu :C α and Ala :C α were respectively defined as the fixed and SMD atoms,and the latter was pulled along the reaction coordinate ( + z ; one-dimensional pulling) with an infinitesimal, constantvelocity of v p = 10 − Å (cid:14) fs, which is practically zero over the total timespan. The force acting on the SMD atomwas described as f ( r , t ) = −∇ U ( r , t ) = −∇ (cid:18) k (cid:2) v p t − ( r − r ) . ˆ n (cid:3) (cid:19) , (1)where U is the harmonic potential (spring), k is the spring constant, t is time, r and r are respectively the instantaneousand initial position vectors of the SMD atom, and ˆ n = (0 , , is the normalized pulling vector (along + z ). Theharmonic biasing potential of the form U ( ξ, i ) = k ( ξ − ξ i ) (cid:14) was added to the Hamiltonian of the system, where k = 2 kcal (cid:14) (mol.Å )—corresponding to a thermal fluctuation of the SMD atom of p k B T /k = 0 . Å at 310 K,with k B the Boltzmann constant— ξ is the collective variable, and ξ i = (31 . i ) Å defines the center of the i th biaspotential. 2hermodynamic equilibrium of biological macromolecules under mechanical constraints A P
REPRINT (a) (b)
Figure 1: (a) The tertiary structure of a single β -strand taken from an S-shaped amyloid beta A β − dodecamer fibriland the associated reaction coordinate (collective variable) defined as the distance between Glu :C α and Ala :C α —rendered in VMD using Tachyon parallel (cid:14) multi-processor ray tracing system [34]. (b) Variations in the initial confor-mation of the β -strand during the 11-Å traction in a way that each strand corresponds with the initial conformation ofone replica. The fixed and SMD atoms at the two ends are shown by the blue and red balls, respectively. Fig. 2 illustrates the initial and equilibrated conformations of each replica shown in cyan (at t = 0 ) and magenta (at t = 1 ns), respectively. As is seen, from replica (cid:14) mol with long tails, demonstrating that the corresponding equilibrium states have been attained. Itcould also be observed that the larger the end-to-end distance is stretched, the higher the value to which the totalenergy converges. Further examining the energy–time diagrams revealed that the specific pattern of the total energywas uniquely determined by that of the electrostatic energy based on the observation that the other energy terms in theCHARMM potential energy function—including bond, angle, dihedral, improper dihedral, and van der Waals (vdW)—fluctuated about constant values over the total timespan. This could be inferred from Fig. 4(a), which depicts the timedependence of these various energy contributions associated to the replica v p = 10 − Å (cid:14) fs. The average value of the force over the second half of the trajectory is also clearly zero due toequilibration as well as the infinitesimal pulling velocity. We proposed a simulation trick to carry out equilibration simulations on biological macromolecules, which are con-stantly under the action of externally-applied, mechanical constraints such as stretching or contraction. Such casescould frequently be encountered in replica-exchange umbrella sampling simulations particularly where the reactioncoordinate is explicitly related to the mechanical deformation of the system and the potential of mean force (PMF) isthe dependent quantity of interest. We accordingly applied our approach to equilibrate a single, stretched β -strand,taken from an amyloid beta dodecamer fibril, where the reaction coordinate was defined as the distance between thetwo ends of the strand, and was accordingly extended by 11 Å in a way that every 1 Å of extension generated theinitial conformation of one replica, leaving 12 replicas as well. We then simulated each replica independently in thecanonical ensemble at 310 K using constant-velocity SMD in a way that the SMD atom was being pulled with a differ-ential velocity of magnitude − Å (cid:14) fs along the reaction coordinate, making this atom to be practically fixed during3hermodynamic equilibrium of biological macromolecules under mechanical constraints A P
REPRINT (a) replica (b) replica (c) replica (d) replica (e) replica (f) replica (g) replica (h) replica (i) replica (j) replica (k) replica (l) replica
Figure 2: Conformational transition of each replica from the associated initial state at t = 0 (cyan) to the equilibratedconformation at t = 1 ns (magenta). The right ends of the two strands for each replica exactly coincide, showing thatthe SMD atom is practically fixed during the simulations.the simulations. As a result, the SMD simulation turned into an equilibrating process and each replica consequentlyequilibrated using free dynamics plus two mechanical constraints corresponding with the two fixed ends. The validityof this kind of constrained equilibration was strongly approved by examining the time-dependent variations of the totalenergy and the force exerted on the SMD atom associated to each replica. References [1] Woo H -J, Roux B 2005
PNAS
Methods Mol. Biol.
J. Med. Chem. Free Energy Calculations: Theory and Applications in Chemistry and Biology
Springer-Verlag: Berlin, Heidelberg.[5] Pathria R K, Beale P D 2011
Statistical Mechanics
Third ed. Butterworth-Heinemann: Oxford.[6] Nussinov R, Tsai C J, Shehu A, Jang H 2019
Molecules J. Chem. Phys. J. Chem. Phys. A P
REPRINT time (ns) t o t a l e n e r g y ( kca l/ m o l ) replica Figure 3: The time-dependent total energy curves of the replicas during constrained equilibration—rendered inGnuplot (version 5.2) [35]. The long-tailed converging patterns are evident, demonstrating that all the replicas arefinally in equilibrium.[9] Widom B 1963
J. Chem. Phys. J. Chem. Phys. J. Comput. Phys. J. Chem. Phys. Molecular Dynamics Simulation: Elementary Methods
First ed. Wiley-Interscience.[14] Landau D P, Binder K 2000
A Guide to Monte Carlo Methods in Statistical Physics
Cambridge University Press:Cambridge.[15] Owicki J C, Scheraga H A 1977
J. Am. Chem. Soc. J. Chem. Phys.
J. Chem. Phys.
J. Phys. Soc. Japan J. Comput. Phys. J. Mol. Graph. Model. Monte Carlo Strategies in Scientific Computing
Springer Science & Business Media.[22] Gao M, Wilmanns M, Schulten K 2002
Biophys. J. J. Mol. Graphics Mol. Simul. A P
REPRINT −4−2 0 2 4 6 8 10 0 0.5 1 e n e r g y ( kca l/ m o l ) time (ns) bondangledihedralimproperelectrostaticvdWtotal (a) −4−2 0 2 4 0 0.2 0.4 0.6 0.8 1 f o r ce ( p N ) time (ns) (b) Figure 4: (a) The various energy contributions of the CHARMM potential energy function including bond, angle,dihedral, improper dihedral, electrostatic, van der Waals (vdW), and total energy, and (b) the associated instantaneousforce exerted on the SMD atom, as functions of time over the total timespan of simulating the replica
Nat. Struct. Mol. Biol. Commun. ACM Using MPI: Portable Parallel Programming with the Message-PassingInterface
MIT Press.[31] Phillips J C, Braun R, Wang W, Gumbart J, Tajkhorshid E, Villa E, Chipot C, Skeel R D, Kale L, Schulten K2005
J. Comput. Chem. J. Phys. Chem. B
J. Mol. Graph. ∼∼