Best Practices in Constant pH MD Simulations : Accuracy and Sampling

This is a sel -a chi ed e sion o an o iginal a icle. This e sion may di e om he o iginal in pagina ion and ypog aphic de ails. Au ho (s): Ti le: Yea : Ve sion: Copy igh : Righ s: Righ s u l: Please ci e he o iginal e sion: CC BY 4.0 h ps://c ea i ecommons.o g/licenses/by/4.0/ Bes P ac ices in Cons an pH MD Simula ions : Accu acy and Sampling © 2022 The Au ho s. Published by Ame ican Chemical Socie y Published e sion Buslae , Pa el; Aho, Noo a; Jansen, An on; Baue , Paul; Hess, Be k; G oenho , Ge i Buslae , P., Aho, N., Jansen, A., Baue , P., Hess, B., & G oenho , G. (2022). Bes P ac ices in Cons an pH MD Simula ions : Accu acy and Sampling. Jou nal o Chemical Theo y and Compu a ion, 18(10), 6134-6147. h ps://doi.o g/10.1021/acs.jc c.2c00517 2022 Bes P ac ices in Cons an pH MD Simula ions: Accu acy and Sampling Pa el Buslae ,* ,# Noo a Aho, # An on Jansen, Paul Baue , Be k Hess,*and Ge i G oenho * J. Chem. Theo y Compu . 2022,10.1021/acs.jc c.2c00516 Ci e This: h ps://doi.o g/10.1021/acs.jc c.2c00517 Read Online ACCESS Me ics & Mo e A icle Recommenda ions * sı Suppo ing In o ma ion ABSTRACT: Va ious app oaches ha e been p oposed o include he e ec o pH in molecula dynamics (MD) simula ions. Among hese, he λ-dynamics app oach p oposed by B ooks and co-wo ke s [Kong, X.; B ooks III, C. L. J. Chem. Phys. 1996,105, 2414−2423] can be pe o med wi h li le compu a ional o e head and h o each ypeence be used o ou inely pe o m MD simula ions a mic osecond ime scales, as shown in he accompanying pape [Aho, N. e al. J. Chem. Theo y Compu . 2022, DOI: 10.1021/acs.jc c.2c00516]. A such ime scales, howe e , he accu acy o he molecula mechanics o ce ield and he pa ame iza ion becomes c i ical. He e, we add ess hese issues and p o ide he communi y wi h guidelines on how o se up and pe o m long ime scale cons an pH MD simula ions. We ound ha ba ie s associa ed wi h he o sions o side chains in he CHARMM36m o ce ield a e oo high o eaching con e gence in cons an pH MD simula ions on mic osecond ime scales. To a oid he high compu a ional cos o ex ending he sampling, we p opose small modi ica ions o he o ce ield o selec i ely educe he o sional ba ie s. We demons a e ha wi h such modi ica ions we ob ain con e ged dis ibu ions o bo h p o ona ion and o sional deg ees o eedom and hence consis en pKaes ima es, while he sampling o he o e all con igu a ional space accessible o p o eins is una ec ed as compa ed o no mal MD simula ions. We also show ha he esul s o cons an pH MD depend on he accu acy o he co ec ion po en ials. While hese po en ials a e ypically ob ained by i ing a low-o de polynomial o calcula ed ee ene gy p o iles, we ind ha highe o de i s a e essen ial o p o ide accu a e and consis en esul s. By esol ing p oblems in accu acy and sampling, he wo k desc ibed in his and he accompanying pape pa es he way o he widesp ead applica ion o cons an pH MD beyond pKap edic ion. ■INTRODUCTION Thanks o imp o emen s in algo i hms, o ce ields, and compu e ha dwa e, molecula dynamics (MD) simula ions ha e become a e sa ile ool o in es iga ing he con o ma- ional landscape o complex biomolecula sys ems a he a omic le el. 1−5 An impo an algo i hmic imp o emen has been he explici inclusion o pH in MD simula ions, 6−12 as pH is an impo an expe imen al pa ame e ha a ec s he s uc u e and dynamics o biomolecules. To p o ide he use s o he GROMACS MD package 13 wi h access o simula ions a cons an pH, we ha e implemen ed he λ-dynamics-based cons an pH app oach by B ooks and co- wo ke s. 9,14 In con as o he p e ious implemen a ion in a o k o GROMACS 3.3, 10 he new implemen a ion, which is desc ibed in an accompanying pape , 15 is e icien and cons an pH MD simula ions can be pe o med wi h abou 25% o CPU, 30−40% o CPU + GPU compu a ional o e head compa ed o no mal MD simula ions i espec i e o he numbe o i a able si es in he sys em. In he accompanying me hodological pape , 15 we also demons a e ha he me hod can be success ully applied o calcula e pKa alues o i a able si es in a p o ein. The pu pose o his pape is o p o ide use s wi h guidelines and ecommenda ions on how o se up and pe o m cons an pH MD simula ions, including he necessa y pa ame iza ion s eps. In cons an pH MD simula ions, i a able g oups can dynamically change hei p o ona ion s a e. These changes a e d i en by in e ac ions be ween he g oup and he chemical en i onmen (modeled wi h a o ce ield) and he aqueous p o on concen a ion (modeled wi h a pH po en ial). Because a he o ce ield le el a numbe o con ibu ions o he ee ene gy o (de)p o ona ion a e no included explici ly, i.e., quan um mechanical in e ac ions associa ed wi h bond b eak- age and o ma ion as well as he ac ual p o on pa icle, co ec ions o he o ce ield a e needed in λ-dynamics-based Recei ed: May 18, 2022 A iclepubs.acs.o g/JCTC © XXXX The Au ho s. Published by Ame ican Chemical Socie y A h ps://doi.o g/10.1021/acs.jc c.2c00517 J. Chem. Theo y Compu . XXXX, XXX, XXX−XXX Downloaded ia UNIV OF JYVASKYLA on Sep embe 19, 2022 a 12:16:48 (UTC). See h ps://pubs.acs.o g/sha ingguidelines o op ions on how o legi ima ely sha e published a icles. cons an pH MD. In GROMACS, hese co ec ions a e implemen ed as analy ical unc ions, VMM(λj), i ed o he ee ene gy p o ile associa ed wi h he dep o ona ion o a i a able esidue ja he o ce ield le el. 15 The accu acy o such ee ene gy p o iles depends no only on how closely he o ce ield model ep esen s he ue po en ial ene gy su ace bu also on he con e gence o sampling o all o he deg ees o eedom in he sys em. The e o e, whe eas in no mal MD he accu acy o he dynamics depends solely on he quali y o he o ce ield, he accu acy o λ-dynamics-based cons an pH MD depends addi ionally on whe he all ele an deg ees o eedom a e sampled su icien ly in he simula ions equi ed o pa a- me izing he co ec ion po en ials. We ound ha insu icien sampling o he dihed al deg ees o eedom in he amino acid side chains can lead o poo con e gence in he dep o ona ion ee ene gy p o iles, as also obse ed by Klimo ich and Mobley in simula ions wi hou cons an pH. 16 We aced he lack o he dihed al sampling o he ba ie s ha sepa a e he minima in he o sion po en ials. These ba ie s a e oo high o each a con e ged sampling o he dihed al ee ene gy landscape on he ime scales o ypical cons an pH MD simula ions. Because he in e ac ion be ween he i a able g oup and he en i onmen depends c i ically on he dihed al angles o he side chain, a lack o con e gence in hese dihed al angles also a ec s he sampling o he p o ona ion s a es. Ra he han inc easing he ime scale o he MD simula ions o ob ain con e ged dihed al and p o ona ion s a e dis ibu- ions o in oducing enhanced sampling echniques, 12,17−21 we p opose educing he ba ie s o dihed al o a ions in a sys ema ic way. We will demons a e ha such op imized dihed al o ce ield pa ame e s imp o e pKaes ima es o amino acids wi hou comp omising he o e all con o ma ional sampling o he p o ein. Wi h a highe accu acy o he unde lying dep o ona ion ee ene gy p o iles, we ound ha he co ec sampling o p o ona ion s a es also depends c i ically on he o de o he polynomial i used o ob ain an analy ical o m o he co ec ion po en ial. We show ha he commonly accep ed i s -o de i , 12 al hough i mly based on linea esponse heo y, 22 is no su icien ly accu a e and can lead o e oneous p o ona ion dynamics in cons an pH MD simula ions. Because he dominan ene ge ic con ibu ion o p o on a ini y comes om elec os a ic in e ac ions, 22 i is o u mos impo ance o use an accu a e desc ip ion o such in e ac ions. Cons an pH MD simula ions ha e been pe o med wi h a ious elec os a ic models, including gene alized Bo n, 9 shi ed cu o , 23 and pa icle mesh Ewald (PME). 10,12 O hese me hods, he Ewald summa ion-based PME me hod 24,25 is gene ally conside ed o p o ide he mos accu a e desc ip ion o he elec os a ic in e ac ions in pe iodic biomolecula sys ems. 26 Because Ewald summa ion can only p o ide accu a e esul s i he simula ion box emains neu al, 27 he cha ge luc ua ions associa ed wi h he dynamic p o o- na ion and dep o ona ion in cons an pH MD simula ions need o be compensa ed o p e en a i ac s. Ti a able si es can be di ec ly coupled o special pa icles, modeled as ions o wa e molecules, 21,28 such ha cha ge is ans e ed di ec ly be ween he i a able si e and ha pa icle. Al e na i ely, all si es can be coupled collec i ely o a su icien ly la ge numbe o bu e pa icles. 29 The la e app oach has he ad an age ha spon aneous luc ua ions in he in e ac ion o he bu e pa icles wi h hei en i onmen a ec all i a able si es o he same ex en . The disad an age is ha he se up and pa ame iza ion o he bu e app oach a e mo e in ol ed, as hese equi e selec ing he numbe o bu e s and pa ame izing hei in e ac ion wi h he es o he sys em. To acili a e he use o bu e s in cons an pH MD, we p o ide a pa a- me iza ion s a egy aimed a p e en ing bu e clus e ing, bu e binding o i a able si es, and bu e pe mea ion in o hyd ophobic egions. We demons a e ha bu e s pa a- me ized wi h his s a egy also a oid ini e-size e ec s associa ed wi h he pe iodici y o small simula ion boxes. 30,31 ■METHODS He e, we go h ough all o he impo an aspec s o he simula ion se ups. Simula ed Sys ems. We pe o med s anda d and cons an pH MD simula ions o he sys ems lis ed in Table 1. The o iginal and modi ied (desc ibed in de ail below) CHARMM36m 32,33 o ce ields we e used in all simula ions. Table 1. Table o Simula ed Sys ems a sys em box size (nm3) no. o wa e s no. o ions o ce ield BUF15×5×5∼4000 11 Na, 11 Cl, 2 Bu CHARMM36m BUF23×3×3∼3000 1 o 2 Bu CHARMM36m ADA 5 ×5×5∼4000 11 Na, 11 Cl, 1 o 10 Bu CHARMM36m, CHARMM36m-cph ADA33×3×3∼850 2 Na, 2 Cl, 10 Bu CHARMM36m-cph ADA77×7×7∼11100 31 Na, 31 Cl, 10 Bu CHARMM36m-cph ADAlow sal 5×5×5∼4000 4 Na, 4 Cl, 10 Bu CHARMM36m-cph ADAhigh sal 5×5×5∼4000 38 Na, 38 Cl, 10 Bu CHARMM36m-cph AEA 5 ×5×5∼4000 11 Na, 11 Cl, 1 o 10 Bu CHARMM36m, CHARMM36m-cph AKA 5 ×5×5∼4000 11 Na, 12 Cl, 1 o 10 Bu CHARMM36m, CHARMM36m-cph AHA 5 ×5×5∼4000 11 Na, 12 Cl, 1 o 10 Bu CHARMM36m, CHARMM36m-cph AAA15×5×5∼4000 11 Na, 12 Cl, 1 o 10 Bu CHARMM36m, CHARMM36m-cph AAA25×5×5∼4000 11 Na, 11 Cl, 1 o 10 Bu CHARMM36m, CHARMM36m-cph 1CVO16.8 ×7.8 ×7.1 ∼12400 35 Na, 48 Cl, 20 Bu CHARMM36m, CHARMM36m-cph 1CVO27.1 ×7.1 ×7.1 ∼11400 13 Cl, 150 Bu CHARMM36m MEMB15.9 ×5.9 ×8.1 ∼4200 8 Na, 8 Cl, 50 Bu CHARMM36m MEMB25.9 ×5.9 ×8.1 ∼4200 8 Na, 8 Cl, 1 Bu CHARMM36m SOL 4.3 ×4.3 ×4.3 ∼2600 50 Bu CHARMM36m a We deno e he modi ied CHARMM36m o ce ield as CHARMM36m-cph. Jou nal o Chemical Theo y and Compu a ion pubs.acs.o g/JCTC A icle h ps://doi.o g/10.1021/acs.jc c.2c00517 J. Chem. Theo y Compu . XXXX, XXX, XXX−XXX B The able also p esen s he box size, he numbe o CHARMM36 TIP3P wa e molecules, 34−36 he numbe o ions, and bu e pa icles included in each sys em. The i ing coe icien s o he VMM co ec ion po en ial o he bu e pa icles we e ob ained wi h sys em BUF1. To ind he op imal cha ge ange and Lenna d−Jones pa ame e s o he bu e pa icles, enhanced sampling simula ions wi h he accele a ed weigh his og am (AWH) me hod we e pe o med on sys em BUF2. Sys ems ADA, AEA, AKA, and AHA a e alanine ipep ides wi h capped e mini and as he cen al esidue aspa ic, glu amic, lysine, and his idine amino acids, espec i ely. AAA1and AAA2sys ems a e alanine ipep ides wi h p o ona ed e mini. C- and N- e mini we e made i a able in AAA1and AAA2sys ems, espec i ely. Two se s o simula ions o he ca dio oxin V p o ein we e pe o med. Sys em 1CVO1was used o calcula e he pKa alues o i a able esidues, while he la ge sys em 1CVO2was used o compu e he adial dis ibu ion unc ion o he bu e pa icles a ound he p o ein. The memb ane sys ems MEMB1and MEMB2con ained 106 1-palmi oyl-2-oleoyl-glyce o-3-phos- phocholine (POPC) lipids. The s a ing coo dina es and opologies o all sys ems a e p o ided as Suppo ing In o ma ion. Simula ion Se up. Pe iodic bounda y condi ions we e applied in all sys ems. Elec os a ic in e ac ions we e modeled wi h he pa icle mesh Ewald me hod, 24,25 while an de Waals in e ac ions we e modeled wi h Lenna d−Jones po en ials which we e smoo hly swi ched o ze o in he ange om 1.0 o 1.2 nm. Simula ions we e pe o med a a cons an empe a u e o 300 K, main ained by he - escale he mos a , 37 wi h a ime cons an o 0.5 ps−1and a a cons an p essu e o 1 ba , main ained by he Pa inello−Rahman ba os a , 38 wi h a pe iod o 2.0 ps. The leap og in eg a o wi h an in eg a ion s ep o 2 s was used. Bond leng hs o hyd ogens in he solu e we e cons ained wi h he LINCS algo i hm, 39 while he in e nal deg ees o he CHARMM TIP3P wa e molecules 35,36 we e cons ained wi h he SETTLE algo i hm. 40 P io o he cons an pH MD simula ions, he ene gy o all sys ems was minimized using he s eepes descen me hod, ollowed by a 1 ns equilib a ion. Cons an pH MD Simula ion Se up. In he cons an pH MD simula ions, he mass o λ-pa icles was se o 5 a omic uni s and he empe a u e was kep cons an a 300 K wi h a sepa a e - escale he mos a o he λdeg ees o eedom 37 wi h a ime cons an o 2.0 ps−1. The single-si e ep esen a ion, de ined and desc ibed in he accompanying pape , 15 was used o Asp, Glu, Lys, C- e , and N- e , whe eas he mul isi e ep esen a ion, also desc ibed in ha pape , was used o His. 15,41,42 The mul isi e ep esen a ion models each physical s a e o g oups wi h chemically coupled i a able si es wi h an independen λ-coo dina e and akes chemical coupling in o accoun by applying he linea cons ain on hese λ-coo dina es equi ing =1 ii g oup g oup . 15 The same pH and biasing po en ials we e used as in Aho e al. 15 In he sampling simula ions o single i a able esidues, he pH was se equal o he pKaand he ba ie heigh o he biasing po en ial was se o ze o. The i a ion o he ca dio oxin V (PDB ID 1CVO 43 ) p o ein was pe o med by unning 10 independen eplicas o 100 ns each o 15 equidis an ly spaced pH alues in he ange om 1.0 o 8.0 using bo h he o iginal and a modi ied CHARMM36m o ce ield. In he i a ion simula ions, he ba ie heigh o he biasing po en ial was se o 7.5 kJ mol−1 o g oups modeled wi h a single-si e ep esen a ion and o 5 kJ mol−1 o g oups modeled wi h a mul isi e ep esen a ion. Re e ence Simula ions. The cons an pH simula ions equi e a co ec ion po en ial VMM(λj) o each i a able esidue j. These co ec ion po en ials a e he in eg als o polynomial i s o he expec a ion alue o ⟨∂V/∂λ⟩λin e e ence s a e simula ions a ixed λ- alues. 15 Thus, a e in eg a ion, an n h-o de polynomial i o ⟨∂V/∂λ⟩λyields an (n+ 1) h-o de polynomial unc ion ha ep esen s VMM(λ). Howe e , in ou implemen a ion, he i o ⟨∂V/∂λ⟩λwas used, a he han VMM(λ). We hus e e o he i ing o de as he o de o he polynomial i o ⟨∂V/∂λ⟩λ. We pe o med he e e ence simula ions as ollows: The pa ial cha ges in he ipep ide sys ems we e linea ly in e pola ed be ween λ=−0.1 and λ= 1.1 wi h a s ep o 0.05. Fo His, all h ee λcoo dina es we e changed unde he cons ain λ1+λ2+λ3= 1. Fo each se o λ- alues, called a g id poin , we pe o med an 11 ns MD simula ion in which he ∂V/∂λjwe e sa ed e e y picosecond and accumula ed. The o al cha ge o he sys em was kep neu al by simul aneously changing he cha ge o a single bu e pa icle. The i ing p ocedu e is desc ibed in ull de ail in he accompanying pape . 15 Dihed al F ee Ene gy P o iles. Because he sampling o p o ona ion s a es is igh ly coupled o he sampling o side chain dihed al deg ees o eedom, we compu ed he ee ene gy p o iles associa ed wi h he o a ion o he dihed als in he side chain o he cen al amino acid in he capped ipep ide sys ems (Table 1) by means o umb ella sampling. 44,45 As he i s s ep, we pe o med 20 ns MD simula ions wi h a ime-dependen po en ial on he dihed al angle wi h a o ce cons an o 418.4 kJ mol−1 ad−2. The cen e o his po en ial was mo ed om 0° o 360°wi h a a e o 18° ns−1. F om hese simula ions, ames wi h dihed al angles closes o 0°, 10°, 20°, e c., we e selec ed as e e ences o he umb ella eplicas. The di e ence be ween he dihed al angle in he selec ed ames and he a ge angle was always below 0.1°. Then, we pe o med 36 umb ella sampling simula ions o 11 ns wi h a ha monic es aining po en ial cen e ed a he e e ence dihed al angle and a o ce cons an o 418.4 kJ mol−1 ad−2. We used he WHAM p ocedu e, 46 implemen ed in GROMACS, 47 o unbias hese umb ellas and ob ain ee ene gy p o iles associa ed wi h he ull o a ion o he dihed al angle. Dihed al Po en ial Ene gy P o iles a he QM and MM Le els. To check he alidi y o he p oposed o ce ield modi ica ions, we compu ed he po en ial ene gy p o iles o he N−Cα−Cβ−Cγdihed al o aspa ic acid wi h capped esidues. The p o iles we e compu ed a bo h quan um mechanical (QM) and molecula mechanical (MM) le els. The QM p o iles we e compu ed a he MP2/6-31+G*le el o heo y using he Fi e ly QC package, 48 which is pa ially based on he GAMESS (US) 49 sou ce code. The MM p o iles we e compu ed o bo h he o iginal and he modi ied CHARMM36m o ce ields. The po en ial ene gy was compu ed o he N−Cα−Cβ−Cγdihed al angle wi h 10° inc emen s. Fo each dihed al alue, he s uc u es we e ene gy minimized p io o po en ial ene gy calcula ion. Accele a ed Weigh His og am Alchemical Simula- ions. Bu e pa icles a e used in cons an pH MD o main ain he neu ali y o he simula ed sys em. Ideally, bu e s should no in oduce any a i ac s due o binding o i a able g oups, binding o each o he , o pene a ing in o hyd ophobic egions. Jou nal o Chemical Theo y and Compu a ion pubs.acs.o g/JCTC A icle h ps://doi.o g/10.1021/acs.jc c.2c00517 J. Chem. Theo y Compu . XXXX, XXX, XXX−XXX C To p e en such beha io , we op imized he cha ge ange and Lenna d−Jones pa ame e s o he bu e s. To his end, we pe o med a se ies o enhanced sampling simula ions wi h he accele a ed weigh his og am me hod (AWH). 50,51 In one se o simula ions wi h wo bu e s in he simula ion box (BUF2, Table 1), we quan i ied he sampling e iciency om he ic ion me ic 50,51 as a unc ion o he absolu e cha ge pe bu e pa icle. In hese simula ions, he cha ge o one bu e was changed om 0 o +0.8 while simul aneously he cha ge o he o he bu e was changed om 0 o −0.8 in o de o main ain neu ali y. The CHARMM36m Lenna d−Jones pa ame e s o sodium we e used o he bu e s in hese simula ions. In he o he se o simula ions, we compu ed he ee ene gy di e ence be ween in oducing a neu al bu e in wa e (BUF2) and inside he hyd ophobic egion o a POPC bilaye sys em (MEMB2,Table 1) o a ious alues o he Lenna d−Jones pa ame e s o he bu e . In hese simula ions, he Lenna d−Jones in e ac ions be ween he bu e and he es o he sys em we e inc eased om nonin e ac ing a λ= 0 o ully in e ac ing (λ= 1) in 10 disc e e s eps. Fo all sys ems, we simula ed 10 eplicas o 10 ns, om which he ee ene gy di e ences and ic ion coe icien s we e ob ained as a e ages o e he eplicas. Dihed al Analysis. Dis ibu ions o side chain dihed al angles in p o eins we e de i ed om publicly sha ed MD ajec o ies o p o eins wi h PDB IDs 1U19, 52 2RH1, 53 2Y02, 54 and 5UEN 55 ob ained om he GPCRMD 56 and SARS-CoV-2 da abases (h ps://co id.molssi.o g/). 57 Fo each ajec o y, he dis ibu ions o he ollowing dihed al angles we e calcula ed: (1) N−Cα−Cβ−Cγin aspa ic acid (2) N−Cα−Cβ−Cγin glu amic acid (3) Cα−Cβ−Cγ−Cδin glu amic acid (4) N−Cα−Cβ−Cγin his idine (5) H−Oϵ2−C−Oϵ1in aspa ic and glu amic acids. In his wo k, we also compu ed he dis ibu ions o hese dihed als om s anda d MD ajec o ies o ca dio oxin V (PDB ID 1CVO 43 ). Compa isons o Dis ibu ions. To compa e wo dis ibu ions Pi(x) and Pj(x), wi h he co esponding cumula i e dis ibu ion unc ions Fi(x) and Fj(x), we used Kolmogo o −Smi no s a is ics (KSS) 58 = [ ]F F F x F xKSS( , ) sup ( ) ( ) i j x i j (1) The la ge he KSS, he less simila he dis ibu ions a e. The dis ibu ions we e conside ed consis en when KSS(Fi,Fj) < 0.03. The KSS was compu ed using he sc ip om he PCAlipids package. 59,60 Ti a ion. To es ima e he pKa alues o i a able g oups om mul iple simula ions a a ious pH alues, we compu ed he a e age ac ion o dep o ona ed ames (Sdep o ) o e all eplicas a each pH alue. Fo a g oup wi h a single i a able si e, his a e age was ob ained as = + SN N N (pH) dep o dep o p o dep o (2) whe e Np o and Ndep o a e he o al numbe o ames in which he si e is p o ona ed and dep o ona ed, espec i ely. Fo i a able si es modeled wi h he single-si e ep esen a ion, we Figu e 1. Dis ibu ions o he λ-coo dina e (A and E) and dihed als (B−D and F−H) in cons an pH MD simula ion o Asp wi h he o iginal (A−D) and modi ied CHARMM36m (E−H) o ce ield. While he simula ions we e pe o med o an ADA ipep ide, only he cen al aspa ic acid is shown o cla i y in he inse o A. (A and D) Dis ibu ions o hi d-o de i s o ⟨∂V/∂λ⟩λob ained wi h he o iginal o ce ield. Di e en colo s co espond o independen eplicas. Dis ibu ions o he λ-coo dina e (A) as well as dis ibu ions o N−Cα−Cβ−Cγ(B) and O 1 -Cγ- O 2 - H 2 (D) dihed als a e no iden ical. Righ column (E−H) shows he dis ibu ions o he modi ied CHARMM36m o ce ield wi h hi d-o de i o ⟨∂V/∂λ⟩λ. Dis ibu ions a e iden ical. Jou nal o Chemical Theo y and Compu a ion pubs.acs.o g/JCTC A icle h ps://doi.o g/10.1021/acs.jc c.2c00517 J. Chem. Theo y Compu . XXXX, XXX, XXX−XXX D conside ed he si e p o ona ed i λis below 0.2 and dep o ona ed i λis abo e 0.8. Fo si es ha a e desc ibed wi h he mul isi e desc ip ion, we conside ed a s a e p o ona ed i he λassocia ed wi h he p o ona ed o m o he esidue is abo e 0.8 and dep o ona ed i he λassocia ed wi h he dep o ona ed o m o he esidue is abo e 0.8. To es ima e he mac oscopic pKa alue o his idine, which con ains wo i a able si es Nϵand Nδ, we calcula ed o each pH alue he a e age ac ion o ames in which he esidue is dep o ona ed a ei he o he wo si es = + + SN NNN (pH) 1 mac o dep o p p (3) whe e Np , N , and N a e he numbe o ames in which λp> 0.8, λϵ> 0.8, and λδ> 0.8. The a e aged ac ions a each pH alue we e i ed o he Hende son−Hasselbalch equa ion = + S1 10 1 K dep o (p pH) a (4) which yielded he pKa alues as i ing pa ame e s. ■RESULTS AND DISCUSSION Fi s , we demons a e ha a lack o sampling o he ele an dihed al deg ees o eedom in amino acid side chains wi h he CHARMM36m o ce ield educes he accu acy o he co ec ion po en ials o λ-dynamics. To o e come hese con e gence p oblems, we modi y he o ce ield by educing he ba ie s in he o sion po en ial and show ha his signi ican ly imp o es he accu acy o he co ec ion po en ials and hence he esul s o cons an pH simula ions, including pKaes ima es, wi hou a ec ing he p o ein con o ma ional dynamics. A e he alida ion o he modi ied o ce ield pa ame e s, we show how he bu e pa icles ha main ain he neu ali y o he simula ion box ha e o be pa ame ized o p e en ini e-size e ec s on p o on a ini ies due o pe iodici y. 30,31 Sampling. Klimo ic and Mobley ha e shown ha calcula ed hyd a ion ee ene gies o single amino acids depend on he s a ing con o ma ion. 16 Because a ew picoseconds ypically su ice o sample bond and angle deg ees o eedom in he amino acid as well as he o a ional deg ees o eedom o he wa e molecules, we specula e ha hei obse a ion implies a lack o sampling in he dihed al deg ees o eedom in he amino acid side chain. The e o e, we sys ema ically analyzed he con e gence o bo h λ-coo dina es and dihed al angles du ing cons an pH simula ions o single amino acids in wa e . We pe o med 100 ns cons an pH MD simula ions a pH = pKao sys ems ADA, AEA, AKA, AHA, AAA1, and AAA2 (Table 1). To enhance he sampling o he λ-coo dina e in hese sys ems, we an he simula ions wi hou a ba ie in he biasing po en ial (Vbias(λ), eq 5 in Aho e al. 15 ). The co ec ion po en ial (VMM(λ), eq 5 in Aho e al. 15 ) was ob ained by i ing a hi d-o de polynomial unc ion o he ⟨∂V/∂λ⟩λ alues o he e e ence ajec o ies. We will show la e ha o accu a e and ep oducible cons an pH MD esul s, a highe o de i is equi ed. Ne e heless, in spi e o i s limi ed accu acy, using he same hi d-o de i o all sys em su ices o sys ema ically compa e he dis ibu ions o he ele an deg ees o eedom and assess hei con e gence. Figu e 2. Dis ibu ions o dihed al angles o which he o sion po en ials we e co ec ed om s anda d MD simula ions. (A) Dihed al dis ibu ions om publicly a ailable ajec o ies 56,57 ( o al simula ion ime ≈10 μs). P obabili y densi y is signi ican only a ound local minima. (B) Dihed al dis ibu ions ob ained om simula ions o ca dio oxin V wi h bo h he o iginal and he modi ied ba ie s o di e en p o ona ion s a es o he i a able esidues. Local minima a e p ese ed, and no addi ional con igu a ions a e obse ed. Jou nal o Chemical Theo y and Compu a ion pubs.acs.o g/JCTC A icle h ps://doi.o g/10.1021/acs.jc c.2c00517 J. Chem. Theo y Compu . XXXX, XXX, XXX−XXX E In Figu e 1A, we show he dis ibu ions o he λ-coo dina e in i e cons an pH MD eplicas o he ADA sys em wi h he o iginal CHARMM36m o ce ield pa ame e s. Dis ibu ions o he λ-coo dina es in he o he sys ems (AEA, AKA AHA, and AAA1and AAA2) a e shown in he Suppo ing In o ma ion (SI, Figu es S1−S5). The dissimila i y be ween he λ-dis ibu ions in he eplicas (maximum KSS be ween eplicas o 0.29, 0.11, 0.04, and 0.095 o ADA, AEA, AHA, and AAA1, espec i ely) indica es a lack o con e gence. In addi ion, he dis ibu ions o he dihed al angles, shown in Figu e 1B−D, a e also no iden ical o all eplicas. Because he e is no ba ie om he biasing po en ial o he λ-coo dina e, we conclude ha he lack o con e gence in λis due o insu icien sampling o he dihed al deg ees o eedom. Fo ce Field Modi ica ions. To o e come he lack o sampling, one can inc ease he simula ion ime o enhance sampling by means o special algo i hms such as eplica exchange MD. 61,62 Replica exchange has been applied equen ly in he con ex o cons an pH MD wi h exchanges be ween eplicas a di e en empe a u e (T-REMD), pH (pH-REMD), o bo h. 20,63−66 Howe e , because he p o o- na ion s a e has li le e ec on he o sion ba ie heigh (Figu e 3), changing he pH ac oss he pH ladde would do li le o enhance he sampling o he dihed al deg ees o eedom. Fu he mo e, REMD me hods a e compu a ionally mo e demanding han pe o ming a single MD simula ion and also p e en access o he dynamical p ope ies o he sys em due o he jumps be ween eplicas. As o some applica ions he dynamical p ope ies may be highly needed, we conside i impo an ha all ele an p o ona ion s a es can be sampled co ec ly in a single cons an pH MD ajec o y. Because he sampling o λ-coo dina es is igh ly coupled o he sampling o he dihed al angles, con e gence wi hin a single ajec o y can in p inciple also be eached by lowe ing he ba ie s o he o sion po en ials. P e iously, modi ica ions o he o ce ield ha e been in oduced o ca boxyl g oups. To imp o e he sampling o he syn and an i con o ma ions o he ca boxyl p o on in Glu, Asp, and he C- e minus, B ooks and co-wo ke s educed he ba ie o his o a ion by a ac o o 8 and also scaled he ca boxyl oxygen adii by 0.95. 9 In con as , G ubmulle and co- wo ke s modi ied his o sion po en ial o p e en sampling he an i-con o ma ion al oge he . 67 Howe e , acco ding o ou analysis (Figu e 1), he e is a lack o con e gence no only in he ca boxyl dihed al angle bu also in he o he side chain dihed al angles. Reducing he o sional ba ie s wi hou a ec ing he o e all sampling o he con o ma ional space is possible only i he egions nea such ba ie s a e spa sely sampled. We, he e o e, analyzed he dis ibu ions o he dihed al angles in he side chains o i a able amino acids in he publicly a ailable ajec o ies o G-p o ein-coupled ecep o s and o SARS-CoV-2 p o eins. 56,57 The dis ibu ions o hese dihed al angles, plo ed in Figu e 2A, e lec he shape o he unde lying o sion po en ials wi h maxima coinciding wi h local minima o he po en ial p o iles. The low densi y nea ba ie s sugges s ha hese ba ie s a e a he high and migh be educed wi hou a ec ing he dihed al dis ibu ions. To achie e con e gence o bo h dihed al angles and λ-coo dina es on a 100 ns ime scale, we hus al e ed he maxima o he o sion po en ials by adding a small dihed al Figu e 3. Modi ica ion o he o sional ba ie s in he Asp side chain. (A) Aspa ic acid and i s a omic nomencla u e. (B and C) Co ec ions added o dihed al o sions o he o iginal o ce ield. Co ec ions wi h wo (B) and h ee (C) local minima we e used o o sions wi h wo and h ee local minima, espec i ely. Heigh s o he co ec ions we e selec ed h ough he i e a i e p ocess, aimed a achie ing consis en λ-dis ibu ions wi hou in oducing addi ional local minima in he ee ene gy p o iles. (D−F) O iginal and modi ied o sional ba ie s o Asp o bo h p o ona ed (H+) and dep o ona ed (H−) s a es. Jou nal o Chemical Theo y and Compu a ion pubs.acs.o g/JCTC A icle h ps://doi.o g/10.1021/acs.jc c.2c00517 J. Chem. Theo y Compu . XXXX, XXX, XXX−XXX F angle (ϕ)-dependen co ec ion o he CHARMM36m o ce ield =V n n( ) cos( ) ( 1) 1 4cos(2 ) i i i n i i Ä Ç Å Å Å Å Å Å Å Å É Ö Ñ Ñ Ñ Ñ Ñ Ñ Ñ Ñ (5) whe e niis he mul iplici y o he o sion angle i(i.e., he numbe o minima) wi h ni= 2 o conjuga ed bonds and ni= 3 o alipha ic bonds. The pa ame e ϵiis an empi ical coe icien ha is op imized such ha he ba ie s a e low enough o con e ge he dis ibu ion o ϕiwi hou in oducing addi ional minima on he po en ial ene gy su ace. Fo each side chain dihed al angle o he i a able amino acids, he coe icien ϵiwas op imized in an i e a i e ashion: A e an ini ial guess, we compu ed he ee ene gy p o iles associa ed wi h o a ion o he dihed al as well as i e unbiased 100 ns ajec o ies a pH = pKawi h di e en s a ing condi ions and a biasing po en ial (Vbias(λ), eq 5 in Aho e al. 15 ) wi hou ba ie . P io o hese cons an pH MD simula ions, we ecompu ed he co ec ion po en ial, VMM(λ), by i ing a hi d-o de polynomial o he ⟨∂V/∂λ⟩λ alues ob ained om he modynamic in eg a ion simula ions pe o med wi h he cu en alue o ϵi. F ee ene gy p o iles we e inspec ed isually o a i icial minima, while dis ibu ions o bo h dihed al angles and λ-coo dina es we e compa ed be ween he i e unbiased eplica uns based on hei simila i y. The coe icien ϵiwas g adually inc eased un il he dis ibu ions in he di e en eplicas we e su icien ly simila (KSS < 0.03), while a he same ime no addi ional minima appea ed in he ee ene gy p o iles. In Figu e 3 we show he op imized o sion po en ials as well as hei e ec on he ee ene gy p o iles o he dihed al angles in Asp. The co ec ions and ee ene gy p o iles o Glu, His, and he C- e minus a e shown in Figu es S6−S8 o he SI. Wi h he excep ion o he Cβ−Cγ−Oδ−H dihed al, he co ec ions in oduce no addi ional minima on he ee ene gy p o ile o hese o sions. Fu he mo e, as shown on he igh panels o Figu e 1, he dis ibu ions o he dihed als and λ-coo dina es a e nea ly indis inguishable o all eplicas a e co ec ion. Because wi h he co ec ed po en ials he Kolmogo o −Smi no s a is ics o Asp, Glu, His, and he C- e minus a e 0.028, 0.015, 0.027, and 0.022, we conclude ha he co ec ions imp o e he con e gence o bo h he λ-coo dina es and he dihed al deg ees o eedom in cons an pH simula ions. No e ha al hough he dis ibu ions o he λ-coo dina e a e su icien ly simila , he sampling o hese coo dina es is no ye uni o m (Figu e 1E). We will show below ha his disc epancy is due o he low o de o he polynomial i used o ob aining he co ec ion po en ial VMM(λ). Wi h he excep ion o he O 1 −Cγ− O 2 −H dihed al, he co ec ions we p opose he e lead o changes in he o sion ba ie o a mos 16 kJ mol−1(i.e., o he Glu N−Cα−Cβ−Cγ o sion). Fo many biomolecula o ce ields, he pa ame e s o he o sion po en ials a e ob ained by i ing sui able pe iodic unc ions o ene gies e alua ed a he MP2 le el o heo y. 68−71 The pa ame e s o each ype o o sional po en ial a e simul aneously i ed o mul iple amino acids. The e o e, he a e age oo mean squa ed (RMS) di e ence be ween he o sional ene gy a he CHARMM36m le el and ha a he MP2 le el o heo y is on he o de o 10 kJ mol−1. The RMS de ia ion be ween he modi ied and he o iginal o sion po en ials is a mos 8 kJ mol−1, and he RMS de ia ion be ween he ab ini io po en ial a he MP2/6-31+G*le el and he N−Cα−Cβ−Cγ o sion po en ial in ASP is educed om 4 kJ mol−1 o he o iginal CHARMM36m o ce ield o 3.5 kJ mol−1 o he modi ied CHARMM36m o ce ield (Figu e S9). The e o e, we conclude ha wi h he co ec ions o he o sion po en ials, he modi ied o ce ield p o ides an equally good i o QM po en ial p o iles as he o iginal CHARMM36m o ce ield. 33,71 We also pe o med s anda d MD simula ions and simula ed i e eplicas o 100 ns o he wo p o ona ion s a es o he Asp ipep ide in wa e using bo h he o iginal and he modi ied CHARMM36m o ce ield pa ame e s. Wi hou he modi ica ions, he local minima a e no consis en ly sampled in all eplicas (Figu e S10). In con as , wi h he co ec ions, iden ical dis ibu ions o he dihed al angles a e ob ained also in s anda d MD simula ions. In addi ion, he modi ica ions a e essen ial o sample bo h syn and an i con o ma ions o he ca boxyl p o on, in line expe imen . 72 We no e, howe e , ha he co ec ion equi ed o sampling bo h o hese con- o ma ions signi ican ly al e s he shape o he ba ie (Figu e 3F). Ne e heless, because o hei low mass, p o on can unnel h ough such ba ie s, and he e o e, we conside he shape and heigh o he o sional ba ie less ele an o his speci ic dihed al han o he o he dihed als. Finally, we demons a e ha he modi ica ions do no al e he dis ibu ions o he dihed al angles in p o ein simula ions. We pe o med MD simula ions o he 1CVO1sys em bo h wi h and wi hou he modi ica ions o he o sion po en ials o i a able amino acids wi h ei he (i) all o hese esidues p o ona ed, (ii) all dep o ona ed, o (iii) all Asp esidues dep o ona ed and all o he esidues p o ona ed. In Figu e 2B, we plo he dis ibu ions o he dihed al angles o which co ec ions we e in oduced. The high simila i y be ween he dis ibu ions sugges s ha he co ec ions do no lead o he sampling o di e en dihed al dis ibu ions, e en i he ela i e weigh s o he minima a e sligh ly al e ed, in pa icula , o he H−O−C−O dihed al. We conclude, he e o e, ha he co ec ions in oduced o acili a e sampling o he dihed al and λ-coo dina es do no signi ican ly al e he p o ein con o ma ional landscape and can hence be used o pe o m bo h no mal and cons an pH MD simula ions. Quali y o he Co ec ion Po en ial VMM.Re e ence Po en ial. To e i y ha he modi ied o sion po en ials o e come he con e gence p oblems, we pe o med cons an pH simula ions a pH = pKaand wi hou a ba ie in he biasing po en ial. Because wi h such a se up he po en ial ene gy p o ile o he λ-pa icle should be la , we expec ed a uni o m λ-dis ibu ion, p o ided ha he dihed al deg ees o eedom a e su icien ly sampled. Howe e , as shown in Figu e 1E, he dis ibu ions a e iden ical be ween eplicas bu no uni o m despi e he co ec ions o he o sion po en ials. Because bo h he pH-dependen po en ial VpH(λ) and he biasing po en ial a e la by cons uc ion a pH = pKa, he de ia ions mus o igina e om disc epancies be ween he co ec ion po en ial VMM(λ) and he unde lying ee ene gy p o ile associa ed wi h dep o ona ion. The co ec ion po en ial is ob ained as a polynomial i o he ⟨∂V/∂λ⟩λ alues om he modynamic in eg a ion simula ions. Because linea e- sponse (LR) heo y p edic s a linea dependence be ween he hyd a ion ee ene gy and he magni ude o a (poin ) cha ge, a i s -o de i has o en been used o ob ain he co ec ion po en ial o cons an pH MD. 9,12 Howe e , e en i he change in he cha ge domina es he ee ene gy o changing he Jou nal o Chemical Theo y and Compu a ion pubs.acs.o g/JCTC A icle h ps://doi.o g/10.1021/acs.jc c.2c00517 J. Chem. Theo y Compu . XXXX, XXX, XXX−XXX G p o ona ion s a e, hyd ogen-bond ea angemen s can con ib- u e as well. Because he e ec s due o such s uc u al ea angemen s a e neglec ed in LR heo ies, we hypo hesized ha highe o de i s may be necessa y o ob aining su icien ly accu a e co ec ion po en ials. To es ou hypo hesis, we in es iga ed he accu acy o he polynomial i o he co ec ion po en ial. In Figu e 4A, we show he mean e o o he co ec ion po en ial wi h espec o he compu ed ee ene gy di e ence associa ed wi h dep o o- na ion as a unc ion o he i ing o de . Fo he LR app oxima ion he i ing e o s a e highe han 10 kJ/mol. In he wo s -case scena io, such e o s could lead o de ia ions in p edic ed pKa alues o mo e han one pKauni . Wi h an e o o 4 kJ mol−1, he hi d-o de i , used abo e o add ess he con e gence issues, does no yield a su icien ly accu a e ep esen a ion o he unde lying ee ene gy p o ile. Inc easing he o de o he polynomial i educes his e o , and as shown in Figu e 4B, a leas a se en h-o de i is equi ed o p o ide a uni o m dis ibu ion o he λcoo dina e o he Asp ipep ide in cons an pH MD simula ions a pH = pKa. Also, o ca boxyl g oups in he side chains o Glu and in he C- e minus, a polynomial i o ⟨∂V/∂λ⟩λo a leas se en h- o de is needed o p o ide a su icien ly accu a e co ec ion po en ial (Figu e 4 and Figu es S11 and S14 in SI). Fo he imidazole ing o His wi h h ee coupled i a able si es, a se en h-o de i su ices as well (Figu e S13), while o he amino bases in he side chain o Lys and he N- e minus, a leas an eigh h-o de i is equi ed (Figu es S12 and S15 in SI). We specula e ha he highe o de i is needed o he la e si es due o he la ge change in he cha ge on he cen al ni ogen a om om −0.3 e o −0.96 eupon dep o ona ion. The change in he cha ge o he ca boxylic oxygen om 0.55 e o −0.76 eis smalle as a e he changes on he ni ogen a oms o he imidazole ing o His ( om −0.36 e o −0.7 e). Pa ame e iza ion o Bu e Pa icles. A change in he p o ona ion s a e a ec s he o al cha ge o he simula ed sys em, which can lead o a i ac s when Ewald summa ion is used o ea he elec os a ic in e ac ions. 27,31 In ou implemen a ion o cons an pH MD, we a oid his p oblem by in oducing i a able bu e s in o he simula ion box ha compensa e o he cha ge luc ua ions o he i a able esidues. 29 In he o iginal implemen a ion o cons an pH MD in GROMACS, 10 he bu e s we e hyd onium molecules ha compensa ed o he o e all cha ge luc ua ions by changing hei cha ge be ween 0 and +1 e. To p e en sampling cha ges beyond his in e al, a biasing po en ial wi h s eep edges a λ= 0 and 1 was in oduced o es ic he ange o λ- alues. Howe e , his po en ial s ill in oduces addi ional o ces a he edges o he λ-in e al. To a oid he e ec s o such o ces, we use a comple ely la biasing po en ial o he bu e s, also ou side o he cha ge in e al. Because changing he cha ge o a bu e pa icle in solu ion induces local ea angemen s o he hyd ogen-bonding ne - wo k ha in u n could a ec he p o on a ini y o a nea by i a able g oup, we wan o minimize he impac o cha ging he bu e pa icles. To de e mine he cha ge ange in which he bu e s do no cause signi ican hyd ogen-bond ne wo k ea angemen , we an AWH simula ions wi h wo ions, he cha ges o which a e changed simul aneously in opposi e di ec ions (BUF2sys em). F om he ic ion me ic a ailable in he AWH me hod, 50 we es ima ed he local di usion coe icien , which is ela ed o he e iciency o sampling: The highe he ic ion, he slowe he dynamics and he mo e sampling is equi ed o each con e gence. We calcula ed he ic ion coe icien (Figu e 5A) o he coo dina e associa ed wi h changing he cha ge on he bu e . Fo cha ges highe han 0.5 e, he ic ion was mo e han 50% highe han ha o ze o cha ges, e lec ing longe co ela ion imes and hence slowe dynamics. We, he e o e, conclude ha he op imal ange o he bu e cha ge is be ween −0.5 eand 0.5 e. The collec i e λ-coo dina e o he bu e pa icles is no es ic ed o a ixed in e al by a wall-like po en ial. To a oid he bu e cha ge exceeding he op imal ange, mul iple bu e pa icles a e needed in he simula ion box. The op imal numbe o bu e s can be calcula ed based on he analysis o cha ge luc ua ions pe o med by Donnini e al. 29 Wi h a small cha ge, a bu e pa icle is apola . To p e en clus e ing o such apola pa icles in wa e , pe mea ion in o hyd ophobic a eas, such as memb ane in e io s, o in e ac ions wi h he p o ein, he Lenna d−Jones pa ame e s (σand ϵ) o he bu e s we e chosen such ha he bu e s ha e only epulsi e in e ac ions wi h all o he a oms, excep wa e . A e expe imen ing wi h he pa ame e s o he bu e pa icles, we se led on a σo 0.25 nm and an ϵo 4 kJ mol−1. This choice leads o dec eased clus e ing o bu e s, low bu e concen- a ions in he p oximi y o i a able si es, and educed pene a ion in o hyd ophobic egions (Figu e 5). The esul ing ee ene gies o neu al bu e inse ion in o wa e and he hyd ophobic egion o he memb ane a e −2.09 ±0.07 and 1.2 ±0.6 kJ mol−1compa ed o 8.45 ±0.05 and 7.8 ±0.3 Figu e 4. Quali y o he VMM(λ) co ec ion po en ial as a unc ion o he o de o he polynomial i o ⟨∂V/∂λ⟩λ o Asp. (A) Fi ing e o as a unc ion o i ing o de (black line). G ay dashed line shows he a e age e o in he calcula ed ⟨∂VMM/∂λ⟩. (B) Dis ibu ions o λ-coo dina es o he hi d- and se en h-o de polynomial i s o ⟨∂V/∂λ⟩λ. Whe eas wi h he lowe o de i he dis ibu ion is signi ican ly ugged, he dis ibu ion becomes nea ly la and uni o m on he [0, 1] in e al o he λ-coo dina e i a se en h-o de i is used. Jou nal o Chemical Theo y and Compu a ion pubs.acs.o g/JCTC A icle h ps://doi.o g/10.1021/acs.jc c.2c00517 J. Chem. Theo y Compu . XXXX, XXX, XXX−XXX H