A generic implementation of replica exchange with solute tempering (REST2) algorithm in NAMD Sunhwan Jo 5/26/2016 NAMD Developer Workshop U of Chicago 

Background: Replica Exchange Algorithm Simulation Time Replicas T1 T2 T3 T4 Simulation X = ⇣ x [i(1)] m(1), x [i(2)] m(2), . . . , x [i(M)] m(M) ⌘ 

Replicas T1 T2 T3 T4 Exchange Attempt Simulation Time Background: Replica Exchange Algorithm X = ⇣ x [i(1)] m(1), x [i(2)] m(2), . . . , x [i(M)] m(M) ⌘ X0 = ⇣ x [i(1)] m(2), x [i(2)] m(1), . . . , x [i(M)] m(M) ⌘ 

Replicas T1 T2 T3 T4 Exchange Attempt Simulation Time Background: Replica Exchange Algorithm X = ⇣ x [i(1)] m(1), x [i(2)] m(2), . . . , x [i(M)] m(M) ⌘ X0 = ⇣ x [i(1)] m(2), x [i(2)] m(1), . . . , x [i(M)] m(M) ⌘ w (X ! X 0 ) = exp ⇥ mU (x j ) nU (x i ) ⇤ exp [ mU (x i ) nU (x j )] = exp ⇥ ( n m)( U (x j ) U (x i )) ⇤ 

Replicas T1 T2 T3 T4 Simulation Time Background: Replica Exchange Algorithm 

Background: Replica Exchange Algorithm Energetic barriers are reduced at high temperature. 300K 600K 

Background: Replica Exchange Algorithm Energetic barriers are reduced at high temperature. w/o REX w/ REX 

Motivation: Solute tempering Optimal number of replica required M ⇡ p f ln T max /T min Fukunishi H, Watanabe O, Takada S. J Chem Phys (2002) 116 Tmax=450 Tmin=300 

Motivation: Solute tempering Decompose system into 2 components, solute and solvent Solute-solute Solute-solvent Solvent-solvent U0( x ) = Uss( x ) + Usv( x ) + Uvv( x ) Liu P, Kim B, Friesner RA, Berne BJ. PNAS (2005) 102 Wang L, Friesner RA, Berne BJ. PNAS (2011) 115 

Motivation: Solute tempering Decompose system into 2 components, solute and solvent U0( x ) = Uss( x ) + Usv( x ) + Uvv( x ) REST Um( x ) = Uss( x ) + 0 + m 2 m Usv( x ) + 0 m Uvv( x ) Liu P, Kim B, Friesner RA, Berne BJ. PNAS (2005) 102 Wang L, Friesner RA, Berne BJ. PNAS (2011) 115 exp[ mUm(x)] = exp  mUss(x) + 0 + m 2 Usv(x) + 0Uvv(x) 

Motivation: Solute tempering REST Um( x ) = Uss( x ) + 0 + m 2 m Usv( x ) + 0 m Uvv( x ) Liu P, Kim B, Friesner RA, Berne BJ. PNAS (2005) 102 Wang L, Friesner RA, Berne BJ. PNAS (2011) 115 Such treatment cancels the effect from solvent degree of freedom in the exchange attempt. w (X ! X 0 ) = exp ⇥ mU (x j ) nU (x i ) ⇤ exp [ mU (x i ) nU (x j )] = exp ⇥ mUj ss ( 0 + m) / 2 Uj sv 0Uj vv nUi ss ( 0 + n) / 2 Ui sv 0Ui vv ⇤ exp h mUi ss ( 0 + m) / 2 Ui sv 0Ui vv nUj ss ( 0 + n) / 2 Uj sv 0Uj vv i = exp  ( n m) ✓ Uj ss + 1 2 Uj sv Ui ss 1 2 Ui sv ◆ 

Motivation: Solute tempering REST Um( x ) = Uss( x ) + 0 + m 2 m Usv( x ) + 0 m Uvv( x ) Liu P, Kim B, Friesner RA, Berne BJ. PNAS (2005) 102 Wang L, Friesner RA, Berne BJ. JCTC (2011) 115 Such treatment cancels the effect from solvent degree of freedom in the exchange attempt. w (X ! X 0 ) = exp ⇥ mU (x j ) nU (x i ) ⇤ exp [ mU (x i ) nU (x j )] = exp ⇥ mUj ss ( 0 + m) / 2 Uj sv 0Uj vv nUi ss ( 0 + n) / 2 Ui sv 0Ui vv ⇤ exp h mUi ss ( 0 + m) / 2 Ui sv 0Ui vv nUj ss ( 0 + n) / 2 Uj sv 0Uj vv i = exp  ( n m) ✓ Uj ss + 1 2 Uj sv Ui ss 1 2 Ui sv ◆ 

Motivation: Solute tempering Such treatment cancels the effect from solvent degree of freedom in the exchange attempt. REST2 Um( x ) = m 0 Uss( x ) + s m 0 Usv( x ) + Uvv( x ) exp[ 0Um(x)] = exp h mUss(x) + p 0 mUsv(x) + 0Uvv(x) i Liu P, Kim B, Friesner RA, Berne BJ. PNAS (2005) 102 Wang L, Friesner RA, Berne BJ. JCTC (2011) 115 

Motivation: Solute tempering Such treatment cancels the effect from solvent degree of freedom in the exchange attempt. REST2 Um( x ) = m 0 Uss( x ) + s m 0 Usv( x ) + Uvv( x ) w (X ! X 0 ) = exp ⇥ 0Um(x j ) 0Un(x i ) ⇤ exp [ 0Um(x i ) 0Un(x j )] = exp ⇥ mUj ss p 0 mUj sv 0Uj vv nUi ss p 0 mUi sv 0Ui vv ⇤ exp h mUi ss p 0 mUi sv 0Ui vv nUj ss p 0 mUj sv 0Uj vv i = exp  ( m n) ✓ Uj ss Ui ss + p 0 p m + p n ( Uj sv Ui sv) ◆ Liu P, Kim B, Friesner RA, Berne BJ. PNAS (2005) 102 Wang L, Friesner RA, Berne BJ. JCTC (2011) 115 

Implementation Two additional parameters are newly added to control scaling. Um( x ) = 1Uss( x ) + 2Usv( x ) + Uvv( x ) REST2 1 = m 0 2 = p 1 REST 1 = 0 m 2 = 0 + m 2 m 

Implementation: REST2 # ComputeNonbondedBase2.h if (sptOn && sptScaleFactor2 == -1) { if(node->molecule->get_spt_type(pExt_i.id) && node->molecule->get_spt_type(pExt_j->id)){ scaling_tmp = sptScaleFactor * scaling;} else if (node->molecule->get_spt_type(pExt_i.id) || node->molecule->get_spt_type(pExt_j->id)){ scaling_tmp = sqrt(sptScaleFactor) * scaling;} else {scaling_tmp = scaling;} } const BigReal A = scaling_tmp * lj_pars->A; const BigReal B = scaling_tmp * lj_pars->B; Nonbond parameters are simply scaled in REST2 

Implementation: REST2 # ComputeNonbondedBase.h BigReal scaling_qi = scaling; if (sptOn && node->molecule->get_spt_type(pExt_i.id) && sptScaleFactor2 == -1) { scaling_qi = sqrt(sptScaleFactor) * scaling;} # ComputeNonbondedBase2.h BigReal scaling_qj = 1.0; if(sptOn && node->molecule->get_spt_type(pExt_j->id) && sptScaleFactor2 == -1){ scaling_qj = sqrt(sptScaleFactor);} BigReal kqq = kq_i * p_j->charge * scaling_qj; # ComputePme.c if(sptOn && node->molecule->get_spt_type(xExt[i].id) && sptScaleFactor2 == -1) { data_ptr->cg = sqrt(sptScaleFactor) * coulomb_sqrt * x[i].charge; } else { data_ptr->cg = coulomb_sqrt * x[i].charge;} Nonbond parameters are simply scaled in REST2 

Implementation: REST2 Bonded parameters are scaled in REST2, but bond and angle may not be scaled (user option). # ComputeSelfTuple.h & ComputeHomeTuple.h int has_spt = sptOn && node->molecule->get_spt_type(t.atomID[0]); for (i=1; i < T::size; i++) { has_spt |= sptOn && node->molecule->get_spt_type(t.atomID[i]); } if (T::size < 4 && !sptScaleAll) has_spt = false; if ( samepatch ) { t.scale = (!has_spt) ? 1.0 : sptScaleFactor; 

Implementation: User Parameters • sptScaleFactor & sptScaleFactor2 are newly added to control scaling. • If sptScaleFactor2 is not set, REST2 is activated automatically. spt on sptScaleFactor 1.0 sptScaleFactor2 1.0 # REST2 if not set sptFile ../solute.pdb sptCol B sptScaleAll yes # Scale bond/angle also? • REST requires additional PME operation to remove image contribution. https://github.com/sunhwan/NAMD-REST 

Implementation: User Parameters 

Application: Enhanced Conformational Sampling of (AAQAA)3 • (AAQAA)3 peptide solvated in ~8000 water molecules (25 K atoms) • Optimal number of replica ~ 200 to span 300K and 600K. Only 16 was used to achieve ~40% acceptance ratio. 

Application: Enhanced Conformational Sampling of (AAQAA)3 • 200 ns of simulation with explicit water was performed. • We observed folding and unfolding transitions during the simulation. 

At the limit of T →∞, REST solute-solvent term survive. REST2 exp[ 0Um(x)] = exp h mUss(x) + p 0 mUsv(x) + 0Uvv(x) i REST exp[ mUm(x)] = exp  mUss(x) + 0 + m 2 Usv(x) + 0Uvv(x) Solute-solvent interaction favors unfolded structure at high energy. Thus, REST2 is more preferable for folding simulation. Application: Enhanced Conformational Sampling of (AAQAA)3 

DOPC-CHOL lipids are selected as a solute Application: Lipid Diffusion Huang K, Garcia AE. JCTC (2014) 10 

DOPC-CHOL are selected as solute molecules. Application: Lipid Diffusion Huang K, Garcia AE. JCTC (2014) 10 

Application: Alchemical/REST2 Binding Free Energy • REST2 is compatible with FEP as long as the perturbed atoms are not selected as solute. • Sampling side-chain around the binding pocket is difficult, so we can apply REST2 to the pocket residues while perturb the ligand to enhance convergence. 

Application: Alchemical/REST2 Binding Free Energy • Because only the end-state matters, we can raise and lower the effective temperature during the alchemical transformation. =0 =1 Teff=T0 Teff=Tmax 

Application: Alchemical/REST2 Binding Free Energy • Because only the end-state matters, we can raise and lower the effective temperature during the alchemical transformation. 

Application: Umbrella Sampling/REST2 Binding Free Energy • Val 111 favors gauche conformation at holo state where as trans conformation is favored at apo state. • If a simulation starts at the holo state (crystal structure), it does not interconvert to trans readily during the FEP simualtion, which lead to overestimation of G. FEP FEP/REST2 

Application: Umbrella Sampling/REST2 Binding Free Energy • Val 111 favors gauche conformation at holo state where as trans conformation is favored at apo state. • If a simulation starts at the holo state (crystal structure), it does not interconvert to trans readily during the FEP simualtion, which lead to overestimation of G. FEP FEP/REST2 

Application: Umbrella Sampling/REST2 Binding Free Energy • Protein-protein binding requires many sidechain-sidechain interaction. • Separation of bound protein complex using umbrella sampling may results in many non-native contacts. • REST2 can be applied to the interface residues to “loosen” these non-native interactions. Abl kinase/ p41 peptide complex 

Application: Umbrella Sampling/REST2 Binding Free Energy Bias-exchange Bias-exchange/ REST2 

Application: Umbrella Sampling/REST2 Binding Free Energy Bias-exchange Bias-exchange/ REST2 

Conclusion • We have implemented a general version solute tempering, which support both original REST and REST2. • REST2 can be used to enhance conformational sampling in explicit simulations, where regular T-REMD requires too many computational resource. • REST2 can be used in free energy method and we hope this could be used to enhance convergence of free energy calculations. 

Roux Lab Benoit Roux Yilin Meng Matthew Pond Donghyuk Suh Ziwei He Other Roux lab members Chris Chipot (UIUC, CNRS) Wei Jiang (ALCF) Acknowledgement