Saturation mutagenesis by efficient free-energy calculation - Journal of

Dec 20, 2017 - Single point mutations in proteins can greatly influence protein stability, binding affinity, protein function or its expression per se...
0 downloads 10 Views 9MB Size
Subscriber access provided by AUBURN UNIV AUBURN

Article

Saturation mutagenesis by efficient free-energy calculation Zuzana Jandova, Daniel Fast, Martina Setz, Maria Pechlaner, and Chris Oostenbrink J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b01099 • Publication Date (Web): 20 Dec 2017 Downloaded from http://pubs.acs.org on December 27, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1

Journal of Chemical Theory and Computation

Saturation mutagenesis by efficient free-energy calculation

2 3

Zuzana Jandova1, Daniel Fast1, Martina Setz1, Maria Pechlaner1,‡, Chris Oostenbrink1*

4 5

1

6

Vienna

7



8

8093 Zürich, Switzerland

Institute of Molecular Modeling and Simulation, University of Natural Resources and Life Sciences,

Current address: Laboratory of Physical Chemistry, Swiss Federal Institute of Technology, ETH,

9 10

*[email protected]

11 12

Abstract

13

Single point mutations in proteins can greatly influence protein stability, binding affinity,

14

protein function or its expression per se. Here, we present accurate and efficient predictions of

15

the free energy of mutation of amino acids. We divided the complete mutational free energy

16

into an uncharging step, which we approximate by a third-power fitting (TPF) approach, and

17

an annihilation step, which we approximate using the one-step perturbation (OSP) method. As

18

a diverse set of test systems, we computed the solvation free energy of all amino acid side-

19

chain analogues and obtained an excellent agreement with thermodynamic integration (TI)

20

data. Moreover, we calculated mutational free energies in model tripeptides and established

21

an efficient protocol involving a single reference state. Again, the approximate methods

22

agreed excellently with the TI references, with a root-mean-square error of only 3.6 kJ/mol

23

over 17 mutations. Our combined TPF+OSP approach does not only show a very good

24

agreement but also a two-fold higher efficiency than full blown TI calculations.

25

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

26

Introduction

27

One point mutation in a carefully chosen position in a protein may have a huge impact on a

28

number of various properties, such as protein stability1, protein secondary structure2, catalytic

29

function,3 oligomerization,4 binding of small ligands,5 DNA6 or protein-protein interactions.7

30

It is of a great interest to understand and be able to predict these effects, for which it is highly

31

relevant to compute the associated free energy of mutation.

32

Plentiful methods to calculate the change of the free energy in different systems have been

33

described. By far the fastest are the empirical scoring functions. They gain their speed by

34

working with rigid molecules, thereby largely neglecting the entropic term.8 A class of

35

relatively fast free-energy methods is formed by the so-called end-point methods. As the

36

name suggests, end-point methods save time by omitting intermediates and simulating solely

37

two end states (e.g. bound and unbound, folded and unfolded, or charged and neutral states).

38

Nonetheless, proper sampling in the end states is necessary for convergence of the free

39

energies and accurate results. Sharir-Ivry et al.9 were able to correctly predict the trend of the

40

free energy of mutation for variants of haloalkane dehalogenase via the linear response

41

approximation (LRA) and Almlöf et al.10 accurately predicted the relative binding free energy

42

of two proteins upon mutations using a similar approach. Other methods, such as molecular

43

mechanics/Poisson-Boltzmann surface area (MM-PBSA)11,

44

generalized Born surface area (MM-GBSA)13 decompose the free energy term into various

45

thermodynamic contributions and have been applied for a quick alanine scan with a good

46

agreement with alchemical methods or experiments.14,

47

theoretically robust are alchemical free-energy calculation methods. They work, in contrast to

48

the end-point methods, with physical or unphysical intermediates, often employing non-

49

interacting dummy atoms. Thus computationally more accessible processes can be simulated

15

2 ACS Paragon Plus Environment

12

and molecular mechanics/

Somewhat less efficient, yet

Page 2 of 41

Page 3 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

50

and – because free energy is a state function – relevant free-energy differences can be

51

calculated from thermodynamic cycles. As the value around the cycle is equal to zero, the

52

nature of intermediate does not influence the final free energy value, but can significantly

53

affect the efficiency of the calculation.16-19 In thermodynamic integration (TI)20 the free energy

54

is computed along a continuous path connecting state A and state B. TI has been proven to be

55

accurate at calculating pKa values,21 relative22, 23 or absolute binding free energies24 of protein-

56

ligand or protein-protein complexes and the free energy of mutation in DNA.25

57

Various methods that are still derived from robust statistical mechanics, but can lead to more

58

efficient free-energy calculations have been reported in the past. One prominent example is an

59

one-step perturbation (OSP) approach in which a judiciously chosen reference state is

60

designed, such that the free-energy difference to multiple relevant end-states can be computed

61

in a single step.26-28 The reference state itself does not have to represent a physical molecule,

62

but the molecular sampling can be enhanced by using soft-core potential energy functions18, 29

63

or floppy molecules.30 With OSP one can compare the free energy difference of two end states

64

A and B from the differences to the reference state. OSP has found applications in design of

65

BACE-1 inhibitors,31 binding free energy calculations of α-thrombin and p38-MAP kinase

66

ligands,32 in combination with quantum mechanics33 and in a number of automated

67

programs.34, 35 Chiang and Wang used OSP to predict the free energy changes of benzene to its

68

derivatives with different substituents on the ring.36 Without using any soft-core potential,

69

however, only predictions for small neutral substituents meet accuracy requirements (< 2.5

70

kJ/mol), while larger changes may be accommodated by using more elaborate reference

71

states.28, 37 Still, OSP was repeatedly shown to perform best for nonpolar changes, while large

72

changes in polarity lead to poor overlap in the conformational ensembles of the reference state

73

and the end-states.18, 37-39



3 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

74

The free energy of charging a neutral compound may be more appropriately approximated by

75

the end-state methods, e.g. based on linear response or the linear interaction energies

76

method.40-42 Some years ago, we introduced the third-power fitting method, which

77

approximates the thermodynamic integration profile for charging or uncharging from

78

simulations at the end-states by using a cumulant expansion to determine second derivatives

79

of the free energy.26, 43 De Ruiter accurately calculated charging free energies of benzamidines

80

in water and trypsin using third power fitting (TPF) method.43 As TPF only takes the two end

81

points of a TI trajectory it still significantly decreases the necessary simulation time while

82

maintaining a comparable accuracy.

83

84

The aim of the current work is twofold. On the one hand, we want to systematically

85

investigate the use of OSP and TPF on a chemically diverse set of compounds, calculating the

86

solvation free energy of the amino acid side chain analogues. On the other hand, we aim to

87

efficiently compute the free energy of mutation in a model peptide, with an eye to

88

computational saturation mutagenesis in future work. We will use the fast and inexpensive

89

OSP approach to calculate the Van der Waals contributions to the free energy (nonpolar free

90

!"# energy, ∆𝐺!!! ) and the third power fitting method (TPF) to calculate the Coulombic

91

!"# contributions to the free energy (uncharging free energy, ∆𝐺!!! ).

92

Our work can be divided into two parts. In the first part, we calculate the free energy of

93

solvation of amino acid sidechains. The free energy of desolvation (negative solvation free

94

energy) can be computed by transforming an amino acid sidechain in solution into a non-

95

interacting dummy particle. While this can be done using an alchemical method using 10-20

96

intermediate states, we here propose to use the scheme in Figure 1. First, the free energy of

97

!" annihilation of the reference state into dummy atoms (∆𝐺!!! ) is calculated via TI, so that the

98

entire desolvation free energy of sidechain analogues is calculated as the sum of the

99

!" , the Van der Waals contribution to the annihilation free energy of the reference state ∆𝐺!!!



4 ACS Paragon Plus Environment

Page 4 of 41

Page 5 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

100

!"# !"# free energy, ∆𝐺!!! , and the Coulombic contributions to the free energy, ∆𝐺!!! . The

101

reference states for OSP calculation of the solvation free energies are formed by one or two

102

soft spheres, noted as R1 and R2, respectively (see Figure 2). As indicated in Figure 1, the

103

!" uncharging and annihilation of all compounds is also calculated using TI, yielding ∆𝐺!!! and

104

!" , respectively. ∆𝐺!!!

105 106

In the second part of this work, we calculate the free energy of mutation of the tripeptide Ala

107

– Ref – Ala, where Ref represents the reference state. In order to compare calculations from

108

the tripeptides with the ones of sidechains, the first two reference states consist of one (R3)

109

and two soft spheres (R4) connected to the peptide backbone. However, the conformational

110

sampling of the side chains, relative to the backbone appears to play a large role, for which

111

we have designed an additional reference state with one non-interacting dummy atom and one

112

soft atom (R5) and one with the same construct but with three additional dummy atoms, R6.

113

The additional dummy atoms in R6 are remnants from previous simulations but should not

114

influence the free energies, such that these simulations can be considered independent repeats.

115

Furthermore, a more elaborate thermodynamic cycle (Figure 3) is devised, employing a

116

conformational library of amino acid sidechains, which was developed in the GROMOS

117

forcefield and can be further used for future applications.

118 119



5 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

120

Methods

121

MD simulations of reference states

122

All MD simulations were carried out using the GROMOS11 software simulation package,44

123

employing the GROMOS 54a8 forcefield.45 The molecular topology building blocks for the

124

reference states (Figure 2) are provided in the Supporting Information. The Van der Waals

125

parameters for soft reference atoms, noted as A, were modified according to Schafer18 to

126

(C6)1/2 = 0.27322 (kJ mol−1 nm6)1/2 and (C12)1/2 = 0.056143 (kJ mol−1 nm12)1/2. This amounts

127

to an interactions between water and the soft reference atoms with a radius, σ = 0.43 nm and a

128

well depth, ε of 0.536 kJ/mol. In the reference state R3 a bond of 0.252 nm was used between

129

the Cα and A, which corresponds to the average distance between Cα and A in R5, in which

130

the sidechain bonds have a length of 0.153 nm and the Cα – D – A angle has an optimal value

131

of 111°. For the reference states R2 and R4 a bond of length of 0.351 nm was used between

132

the soft atoms. This distance was calculated as the regular CH-CH bond (0.153 nm)

133

subtracted from twice the distance of Cα and A in R3 (2 x 0.252 nm). To enhance the

134

sampling of the reference a Lennard-Jones soft core parameter of 1.51 was used for the

135

reference atoms.18, 29

136

Initial structures of the reference states were modelled in MOE.46 Reference states were

137

energy-minimized in vacuum using the steepest-descent algorithm and subsequently solvated

138

in a rectangular, periodic and pre-equilibrated box of simple point charge (SPC) water.47

139

Minimum solute to wall distances were set to 1.8 nm and the minimum solute-solvent

140

distance to 0.23 nm. This led to systems containing around 6000 atoms.

141

During equilibration, initial velocities were randomly assigned according to a Maxwell–

142

Boltzmann distribution at 60 K. All solute atoms were positionally restrained with a harmonic

6 ACS Paragon Plus Environment

Page 6 of 41

Page 7 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

143

potential using a force constant of 2.5 × 104 kJ mol−1 nm−2. In each of the four subsequent

144

20 ps MD simulations, the force constant of the positional restraints was reduced by one order

145

of magnitude and the temperature was increased by 60 K. Subsequently, the positional

146

restraints were removed and rototranslational constraints were introduced on all solute

147

atoms.48 In order to ensure that the systems are equilibrated, the simulations at 300K were

148

prolonged for another 10 ns while keeping the temperature (300 K) as well as the pressure

149

constant (1 atm). To sustain a constant temperature, we used the weak-coupling thermostat49

150

with a coupling time of 0.1 ps. The pressure was maintained using a weak coupling barostat

151

with a coupling time of 0.5 ps and an isothermal compressibility of 7.627 x 10-4 kJ-1 mol nm3

152

for reference states R1 and R2 and 4.575 × 10–4 kJ−1·mol·nm3 for reference states R3-R6.

153

Solute and solvent were coupled to separate temperature baths. Implementation of the

154

SHAKE algorithm50 to constrain bond lengths of solute and solvent to their optimal values

155

allowed for a 2-fs time-step. Nonbonded interactions were calculated using a triple range

156

scheme. Interactions within a short-range cutoff of 0.8 nm were calculated at every time step

157

from a pair list that was updated every fifth step. At these points, interactions between 0.8 and

158

1.4 nm were also calculated explicitly and kept constant between updates. A reaction field51

159

contribution was added to the electrostatic interactions and forces to account for a

160

homogenous medium outside the long-range cutoff using a relative dielectric constant of 61

161

as appropriate for the SPC water model52. Coordinate and energy trajectories were stored

162

every 0.5 ps for subsequent analysis. Production runs of reference states R1 and R2 in water

163

were subsequently performed for 10 ns, while the production runs of tripeptides (reference

164

states R3 – R6) were performed for 50 ns.

165 166

MD simulations of tripeptides for conformational library

7 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 41

167

To obtain a conformational library of the real amino acids, simulations of all Ala – X – Ala

168

tripeptides except for glycine and proline were performed. We followed the simulation

169

protocol of the reference states, with production simulations of 100 ns. These trajectories

170

were subsequently clustered using the algorithm by Daura53 according to the conformation of

171

the sidechain of the second residue after a rotational fit on the backbone atoms. The cut-off

172

for clustering varied depending on the size of the sidechain, so that the first five clusters cover

173

~90 % of trajectories. This clustering led to five central member structures (CMSs) for each

174

amino acid which were subsequently used for fitting in OSP.

175 176

Thermodynamic Integration

177

Solvation free energies of amino acid sidechain analogues and mutation free energies for

178

tripeptides were calculated in two steps. In the first step the uncharging of the molecules was

179

!" performed (∆𝐺!!! ) which was followed by disappearing of the van der Waals radii, i. e.

180

!" !" turning the neutral sidechains into dummy atoms (∆𝐺!!! ) or to the Alanine residue (∆𝐺!!! ).

181

The side chain analogues were prepared by breaking the Cα - Cβ bond and increasing the

182

number of (united atom) hydrogen atoms on Cβ. Thermodynamic integration allows us to

183

calculate the free energy difference between state A and B via multiple discrete intermediate

184

steps using a coupling parameter 𝜆. 1 𝜕𝐻 𝜆 𝑑𝜆 𝜕𝜆 𝜆

!" ∆𝐺!!! = 0

185

(1)

186

In state A, represented by Hamiltonian 𝐻! , 𝜆 is equal to 0 and in state B, represented by 𝐻! ,

187

λ is equal to 1. These two states are linked via 𝐻(𝜆) and the free energy difference between

188

them is calculated from the derivative of the free energy with respect to 𝜆.

189

The uncharging simulations were performed at 11 evenly spaced 𝜆 points, using a linear

190

coupling of the states, without soft-core potentials,

8 ACS Paragon Plus Environment

Page 9 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

𝐻!!! 𝜆 = 1 − 𝜆 𝐻! + 𝜆𝐻!

191 192

(2)

such that !"

193

!" !

=

!" !" !

= 𝐻! − 𝐻!

!

= − 𝑉!"!" ! .

(3)

194

𝐻! and 𝐻! indicate the Hamiltonians from the charged and neutral states,

195

ensemble average obtained from a simulation at λ, and 𝑉!"!" is the electrostatic energy of the

196

molecule with its surroundings. The intramolecular interactions are not included in the

197

equation, as they are not expected to change significantly with the surrounding and cancel in

198

the thermodynamic cycle to compute e.g. the solvation free energy or the relative free energy

199

of mutation in the folded and unfolded state.

200

20 ps of equilibration at each 𝜆 value were followed by 1 ns of production run. To reach an

201

estimated error smaller than 1.0 kJ/mol in the TI-calculations for every amino acid analogue,

202

two more 𝜆 points (0.05 and 0.95) were added to the systems where the solutes carry either a

203

positive or negative net charge. No counter charge was introduced to neutralize these systems

204

upon creation or annihilation of the full charge. The simulations at the end states were

205

!"# prolonged to 10 ns and used to estimate the uncharging free energies with LRA (∆𝐺!!! ) and

206

!"# TPF (∆𝐺!!! ) as well (see below). The calculations of the free energies of annihilation were

207

performed using a similar scheme. Here, soft-core parameters of 0.5 for the van der Waals29

208

were used for the perturbed atoms, such that a nonlinear path in λ is taken.

!

denotes an

209 210

Linear response approximation

211

To approximate the electrostatic contribution to the free energy, ∆𝐺!!! , we applied the linear

212

response approximation (LRA).54,55 !"# ∆𝐺!!! = 𝛽 𝐻! − 𝐻!

213



!

+ 𝐻! − 𝐻!

!

= 𝛽( −𝑉!"!"

9 ACS Paragon Plus Environment

!

+ −𝑉!"!"

!)

(4)

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 41

214

The neutral and charged sidechains are simulated explicitly to obtain the appropriate

215

ensemble averages in equation (4). According to the linear response theory, the parameter β

216

takes a theoretical value of !!.54,55 The term −𝑉!"!"

217

preorganization energy.56 Note that equation (4) with β = !! gives exactly the same result as

218

equation (1) with the parameterization of H(λ) in eqs. (2) and (3) if −𝑉!"!"

219

The LRA estimates were based on 10 ns of simulations in the charged and neutral states,

220

which were the prolonged end states of the TI calculations.

!

is also referred to as the electrostatic

!

is linear in 𝜆.

221 222

Third power fitting

223

In the work of de Ruiter43 we observed that the derivative of the electrostatic contribution to

224

the free energy with respect to 𝜆 does not always correlate linearly with 𝜆. Rather, it shows a

225

slight curvature, which can be described with a third order polynomial. In the third-power

226

fitting approach, we approximate the integration in equation (1) by !"# ∆𝐺!!! =

227

228

! !" 𝑑𝜆 ! !"

=

! 𝑓 !

𝜆 𝑑𝜆

(5)

with

229

𝑓 𝜆 = 𝑎𝜆! + 𝑏𝜆! + 𝑐𝜆 + 𝑑

230

and the parameters are fitted to the first and second derivative of the free energy with respect

231

to 𝜆 in the end states (𝜆 = 0 and 𝜆 = 1). Using the Hamiltonian of equation (2), the first

232

derivative is given by equation (3) and through a cumulant expansion,57 the second derivative

233

of the free energy, is given by !!!

234

!!! !



=

! !! !

( 𝑉!"!"

! !

− (𝑉!"!" )! ! )

10 ACS Paragon Plus Environment

(6)

(7)

Page 11 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

235

Here, k is Bolzmann’s constant, and T is the temperature in K. The second derivative of the

236

free energy with respect to 𝜆 is equal to the negative of the fluctuations of the derivative of

237

the Hamiltonian with respect to 𝜆.

238

The TPF estimates were based on 10 ns of simulations in the charged (𝜆 = 0) and neutral

239

states (𝜆 = 1), which were the prolonged end states of the TI calculations.

B

240 241

One-step perturbation

242

The one-step perturbation is based on Zwanzig’s free-energy perturbation equation:58 !"# = 𝐺! − 𝐺! = −𝑘! 𝑇ln 𝑒 !(!! !!! )/(!! !) ∆𝐺!!!

243

!

(8)

244

where 𝐻! refers to the Hamiltonian of an end state and 𝐻! to the Hamiltonian of a reference

245

state and

246

relevant end state is the neutral amino acid sidechain N. The efficiency of this method results

247

from the fact that there is only one reference state simulation needed to compute the free-

248

energy differences to multiple end states. However, the appropriate reference state has to be

249

chosen very carefully. One precondition for accurate free-energy estimates is that all of the

250

reference state Hamiltonian singularities should be shared by the Hamiltonians of all end

251

states,59 i.e that all the relevant conformations of the end states should be covered by the

252

reference state. To avoid intermolecular singularities between the end state atoms and their

253

surroundings, we use the soft-core interaction.29 Liu et al.26, 60 and Schäfer18 showed that

254

atoms with a soft-core potential in reference state simulation help to accurately predict the

255

solvation free energies of nonpolar molecules in water.

256

In the current work 10 ns MD trajectories of the R1 or R2 reference state was used for the

257

OSP calculation of the amino acid side chain analogues. In case of R1, the centre of geometry

258

of the sidechain was fitted onto the centre of geometry of the reference state atom. In case of

!

to the ensemble average from a simulation of the reference state. Here, the

11 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 41

259

R2, the line between the two atoms with the longest intramolecular distance in the sidechain

260

analogue was aligned with the bond between the two atoms of the reference state. This way,

261

we ensured that the amino acid sidechain would fit into both atoms of the reference state.

262

Additionally, to improve statistics, this fit was repeated every 36 degrees, leading to 10

263

rotations for every sidechain.38

264

For the tripeptides, 50 ns of simulations of the reference states R3 – R6 were used. Because of

265

the common backbone, the fitting was not that straightforward. In a first step, the Cα atom of

266

the fitted amino acid was aligned with the Cα atom of the reference state. In the second step,

267

the centre of geometry of the sidechain was aligned to the center of geometry of the soft

268

atoms in the reference state. Next, the interaction energies were calculated for the sidechain

269

conformation attached to the backbone of the reference state.

270

In the tripeptides, the fitting procedure was repeated for each of the central member structures

271

obtained from separate plain MD simulations, to include multiple conformations of the side

272

!"# chains. This procedure leads to five separate ∆𝐺!!!" values per sidechain. The relative free

273

energies between the conformations 𝑖, was computed using !"#$

∆𝐺!

274

= 𝑘! 𝑇ln𝑃!

(9)

275

where Pi is the relative occurrence of this conformation, as obtained from the conformational

276

clustering. ∆𝐺!

277

!"# obtain a single estimate of ∆𝐺!!!

!"#$

!"# is added to ∆𝐺!!!" and the resulting sum was exponentially averaged to

!"# ∆𝐺!!!

278

= −𝑘! 𝑇ln

!

!"#$ ! ∆!!"# !!!" !∆!!

!

!! !

(10)

279

This equation represents the proper averaging of the multiple paths from the reference state to

280

the real, neutral state of the amino acid in Figure 3.

281



12 ACS Paragon Plus Environment

Page 13 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

282

Results and discussion

283

Solvation free energies

284

Uncharging free energies

285

In the first step, we calculated the solvation free energies of amino acid sidechain analogues

286

in water. As outlined in Figure 1, this process is split up into an uncharging and a cavity-

287

formation, or decoupling step. For the uncharging process, we compare two end-point

288

methods, LRA and TPF, with TI (Table 1). Table S1 shows the individual ligand surrounding

289

energies and their fluctuations in the neutral and charged states, which are used to compute

290

the TPF free-energy differences. The convergence of the fluctuations is represented

291

graphically in figure S1. All fluctuations seem well converged within about 5 ns. For all

292

molecules studied, the agreement between TPF and TI is significantly better than between

293

LRA and TI. As a representative example, Figure 5 shows the free-energy profiles for the

294

uncharging of methionine. It can be seen that the curve of dG/dλ vs. 𝜆 (in black) is not strictly

295

linear, but shows some curvature. While LRA assumes linearity, with TPF we approximate

296

the TI profile better, leading to a better free-energy estimate from the end-state simulations.

297

!"# The total root mean square error (RMSE) over all compounds in Table 1 between ∆𝐺!!! and

298

!" !"# !" ∆𝐺!!! amounts to 11.8 kJ/mol, while between ∆𝐺!!! and ∆𝐺!!! this is reduced to 3.3

299

kJ/mol. The only two cases, where TPF deviates from TI by more than 4 kJ/mol are the

300

compounds with a full positive charge (Arg and Lys sidechain analogues). TPF outperforms

301

LRA especially in the case of compounds bearing a full negative charge (Asp and Glu

302

sidechain analogues), where LRA deviates from TI by around 20 kJ/mol.

303

Decoupling free energies

304

The second step towards the solvation free energy involves the cavity formation step, here

305

computed as the decoupling of the neutral sidechain analogues by turning off the Van der

306

Waals interactions of all atoms with their surroundings. Following Figure 1, this process was

307

!" and a combination of TI and OSP, which is performed by the TI approach, yielding ∆𝐺!!!



13 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

308

!"# !"# !" !" computed as ∆𝐺!!! = ∆𝐺!!! + ∆𝐺!!! . The value of ∆𝐺!!! was -11.5 ± 0.1 kJ/mol for

309

R1, and -16.5 ± 0.2 kJ/mol for R2. Note that these values only need to be computed once

310

to obtain the complete solvation free energy, while it is not needed to represent relative

311

solvation free energies. The annihilation free energies are collected in Table 2 and

312

compared between TI and OSP with reference states R1 and R2. For reference state R1 it

313

!"# is clear that ∆𝐺!!! amounts to extremely high values for larger molecules, such as the

314

sidechain analogues of Arg, Gln, Met, Tyr and Trp. This is a clear indication that no relevant

315

conformations for these compounds are observed in the simulation of R1. Table S3 in the

316

Supporting Information shows the percentage of snapshots in the reference state simulations

317

that contribute significantly to the free-energy estimates. This value is obtained by counting

318

!"# + 𝑘! 𝑇. the number of configurations for which (𝐻! − 𝐻! ) in equation (8) is less than ∆𝐺!!!

319

For the indicated compounds, no significant number of contributing snapshots is obtained,

320

suggesting that the extremely high free-energy estimates.

321

In order to accommodate larger compounds better, we introduced reference state R2 and

322

included 10 rotational states of the molecules in the calculations, as outlined in the Methods

323

!" !"# section. Table 2 shows that the agreement between ∆𝐺!!! and ∆𝐺!!! is much better for the

324

large compounds, giving a RMS Error over the indicated 5 compounds of 4.4 kJ/mol. On the

325

other hand, some of the smaller compounds (e.g. the sidechain anlogues of Cys and Ser) are

326

not so well accommodated in the larger reference state, yielding deviations of roughly 20

327

kJ/mol. A viable strategy to choose the optimal reference state for a compound seems to be to

328

pick the larger reference state, only if the percentage of contributing frames from Table S3 is

329

less than 0.1 % for simulation R1. The RMSE between the TI and OSP data is reduced to 3.1

330

kJ/mol when following this rule of thumb.

331

332

Combining the uncharging and the decoupling simulations into the full transfer of the

333

sidechain analogues from the hydrated to the vacuum state, we can compare the solvation free

334

energies from experiment61 to the value obtained with TI and using the OSP-TPF approach.

14 ACS Paragon Plus Environment

Page 14 of 41

Page 15 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

335

Table 3 collects this data, using the reference states R1 and R2 in the OSP approach. Note that

336

for the sidechains of Ala, Leu, Ile and Val no TPF contribution was added, as these sidechains

337

are completely neutral in the united atom GROMOS 54a8 force field. No experimental data is

338

given for the charged compounds (Arg, Asp, Glu, and Lys), as both the experimental as the

339

computational ones require extensive processing before they can be compared directly.62, 63 In

340

this work, we did not attempt to compensate for the artifacts that arise in free-energy

341

calculations involving a full charge change. These artifacts are identical between the various

342

approaches and are hence irrelevant for a comparison between the TI and TPF+OSP data. The

343

large deviations between the TI data and the TPF+OSP approach are due to the corresponding

344

deviations in the decoupling states. Figure 6 compares the solvation free energies calculated

345

via TI and TPF+OSP to the experimental values. Except for the bulky outliers, which could

346

not be shown in the figure, the calculations using the R1 reference state agree better with the

347

experimental as well as with the TI results. Using the same rule of thumb as above, applying

348

reference state R2 only when the percentage of contributing frames drops below 0.1 % for R1,

349

the RMSE between TI and TPF+OSP is reduced to 2.6 kJ/mol, due to a fortuitous cancellation

350

of the small errors in OSP (RMSE 3.1 kJ/mol) and TPF (RMSE 3.3 kJ/mol). The largest

351

remaining deviations are for Trp (using R2; 6.4 kJ/mol), Lys (using R1; 4.4 kJ/mol) and Glu

352

(using R1; 4.3 kJ/mol). For Trp the deviation still largely comes from the OSP calculation, as

353

this compound is still relatively large for the (larger) reference state as reflected by the lowest

354

percentage of contributing frames for R2 in Table S3. For Lys, however, the deviation largely

355

comes from the TPF calculation (8.9 kJ/mol), one of the two compounds with a full positive

356

charge, for which TPF was seen to perform worst (even though still better than LIE or LRA).

357

For Glu, finally, the deviation of 4.3 kJ/mol is the result of errors of 2.8 kJ/mol for the

358

charging free energy, which is quite reasonable for generating a full charge, and 1.5 kJ/mol

359

for the decoupling free energy.

360 361

Mutation free energies



15 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

362

Eventually, we aim to calculate the free energies of mutation in a protein, therefore, we here

363

validate the TPF+OSP method using a tripeptide of Ala – Ref – Ala. All of the amino acids

364

except for proline and glycine were mutated into Ala via TI or the combination of TPF+OSP,

365

following the scheme in Figure 3.

366

TI was performed in the tripeptide Ala – X – Ala, X being the mutated amino acid, in two

367

separate steps. Firstly, the partial charges of amino acid sidechains were removed followed by

368

the changing of Cß atom of the sidechain into a united CH3 atom and annihilation of all

369

remaining sidechain atoms into dummy atoms. End points of the uncharging TI, i. e. λ=0 and

370

λ=1 were simulated for 10 ns and used for the LRA and TPF calculation.

371

Simulations of the tripeptide with four reference states, R3-R6 were carried out for 50 ns in

372

water. Seeing that in a tripeptide the conformation of a sidechain can, through interactions

373

with the backbone, significantly influence the mutation free energies, we created and used a

374

library of the most relevant sidechain conformations. A simulation of 100 ns of each

375

tripeptide Ala – X – Ala was performed and subsequently clustered to obtain the five most

376

relevant sidechain conformations for each amino acid. CMSs from the five clusters were then

377

fitted into the reference states to apply the OSP. The contributions for the individual

378

conformations and a contribution due to the size of the cluster were exponentially averaged

379

using equation (10) to obtain an overall estimate. An example of the conformational diversity

380

of the five CMSs of Arg is shown in Figure 4.

381

Uncharging free energies

382

Table 4 shows the comparison of uncharging free energies using TI, LRA and TPF and Table

383

S2 shows the individual ligand surrounding energies and their fluctuations in the neutral and

384

charged states, which are used to compute the TPF free-energy differences. The results are

385

very comparable to the values obtained for the individual sidechain analogues. As before, the

386

LRA values agree poorly, with an overall RMSE value of 12 kJ/mol. Similarly, we observe a

16 ACS Paragon Plus Environment

Page 16 of 41

Page 17 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

387

good agreement between TPF and TI uncharging free energies with an overall RMSE of 3.4

388

kJ/mol. The biggest differences, of more than 5 kJ/mol are for Arg, Gln and Lys. Amino acids

389

with a full charge show large uncharging free energies of more than 100 kJ/mol. Surprisingly,

390

we found comparable values for the neutral residues Asn and Gln, which uncharging free

391

energies reaching more than 250 kJ/mol. This could be traced to very strong intramolecular

392

interactions of the acetamide group.

393

Free energies of annihilation

394

In the calculation of the free energies of annihilation, we compare the value obtained by TI,

395

!" !"# !"# !"# ∆𝐺!!! to the corresponding free energy obtained from OSP as ∆∆𝐺!!! = ∆𝐺!!! − ∆𝐺!!! in

396

Table 5. Following up on the OSP calculation for the solvation free energies in the previous

397

section, we here included reference states R3 and R4. Considering that the smaller R1 gave

398

reasonable results for most of the sidechains, we furthermore introduced two additional

399

reference states, R5 and R6. Both of them contain one dummy atom between the soft sphere

400

and Cα, for the sake of a broader sampling of the soft sphere with respect to the backbone. As

401

noted before, R6 is identical to R5, except for some additional dummy particles that may

402

influence the dynamics, but not the conformational ensemble. All five CMSs were fitted into

403

the configurations collected for the four reference states in the Ala – Ref – Ala tripeptide and

404

interaction energies were recalculated. Subsequently, the results from the individual CMS

405

were averaged into one final value, using equation (10). This approach worked reasonably

406

well for most of the sidechains, but in Ile we could see a discrepancy between TI and OSP

407

approach by more than 4 kJ/mol in all reference states. This could be avoided by including

408

eight CMSs instead of five, better representing the conformational ensemble of Ile. Also the

409

largest sidechain of Trp seems most difficult to accommodate in all of the reference states,

410

with deviations from the TI data of 5.7 to 16 kJ/mol.

411

Overall, Table 5 shows that reference state R3 shows the worst agreement between TI and

412

OSP with RMSE of 5.6 kJ/mol, followed by R4 with RMSE of 4.6 kJ/mol. R5 and R6, on the

17 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

413

other hand, show very good agreement between the OSP and TPF values. R5 shows excellent

414

agreement with TI, with RMSE of 2.2 kJ/mol, which is below kBT at 300 K.

415

Free energies of mutation

416

Combing the data from the uncharging and the annihilation of nonpolar sidechains, the total

417

free energy of mutation to Ala was calculated and correlated in Figure 7. The difference to the

418

𝑻𝑷𝑭!𝑶𝑺𝑷 !𝑻𝑰 TI data ( ∆∆𝑮𝑸!𝑨 ) is given in Table 6. Following the trend observed for the

419

annihilation of nonpolar particles with OSP, the reference states showing the best agreement

420

with TI are R5 and R6, with RMSE of 3.2 kJ/mol and 4.1 kJ/mol, respectively. R3 and R4

421

show RMSE of 6.7 kJ/mol and 6.4 kJ/mol, respectively. In all cases the biggest outlier is Lys,

422

with a deviation of more than 9 kJ/mol, Gln, with a deviation of at least 6 kJ/mol and Trp,

423

with a difference of at least 5 kJ/mol. The reasons for these deviations are different. As for the

424

amino acid side chain analogues, the full positive charge, hampers the accuracy of the TPF

425

approach, while for Trp the large size limits the accuracy of the OSP approach. Not

426

surprisingly, the large TPF estimates for Gln described above, lead to a large remaining

427

deviation in the mutation free energy. An average overall deviation that is on the order of 4

428

kJ/mol (~ 1 kcal/mol) is quite acceptable considering the wide range of the mutational free

429

energies in Figure 7.

430 431 432

Efficiency of the method

433

The combined TPF+OSP approach has the potential to significantly reduce the amount of

434

computational time needed for calculation of the free energy of mutation. In TI, 10-20

435

simulations at different l-values are needed for every pair of mutants that is studied. Most of

436

the simulation time is spent on unphysical intermediate states. In the TPF+OSP approach, the

437

charging free energy is collected from only two simulations in the end-states, of which the

438

charged state is a simulation of the actual mutants one is interested in, giving access to

18 ACS Paragon Plus Environment

Page 18 of 41

Page 19 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

439

structural properties of the mutant. Furthermore the OSP contribution is calculated from a

440

single (longer) simulation of the reference state, which is efficiently re-used for all mutants.

441

Taken together, the TPF+OSP approach requires less overall simulation time, and a larger

442

fraction of the simulation time is spent on simulations that have physical relevance and can be

443

further used to study e.g. the interactions of the amino acids.

444

Without having optimized the simulation time for either of the approaches, Table S4

445

summarizes the overall simulation time used in the TI calculations and the TPF+OSP

446

calculations of the tripeptides. The latter approach was based on approximately half the

447

simulation length as the TI data.

448

It is important to mention that the computational resources spent on creating of the sidechain

449

conformation library amounted to 100 ns for each amino acid. This library is however,

450

already created and can be further used in future calculations. Alternatively, if necessary, the

451

end point trajectories of TPF could be used to create additional libraries in future calculations

452

of the free energy of mutation in very different systems.

453



19 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

454

Conclusion

455

In the current work, we systematically studied the performance of an approach combining

456

third-power fitting (TPF) and one-step perturbation (OSP) approaches, as compared to more

457

robust thermodynamic integration (TI) data. The solvation free energies of a large range of

458

compounds (charged, polar, nonpolar, small aliphatic, aromatic, …) were computed and

459

compare excellently to the TI data. When two reference states are considered, depending on

460

the size of the solute, a root-mean-square deviation of only 2.6 kJ/mol was obtained.

461

Furthermore, we extended this approach to the calculation of mutational free energies in

462

model tripeptides and could establish an efficient protocol involving a single reference state.

463

The overall deviation between the TI data and the TPF+OSP approach amounts to a very

464

good 3-4 kJ/mol, which is still very acceptable considering the large range of free energies

465

considered. Without explicit optimization of the method, the TPF+OSP approach can be

466

about two times more efficient than full TI calculations.

467 468



20 ACS Paragon Plus Environment

Page 20 of 41

Page 21 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

469

Acknowledgements

470

Financial support of the Austrian Science Fund, project I 1999 is gratefully acknowledged.

471 472 473

Supporting information

474

Tables S1 and S2 contain the individual energy terms and their fluctuations as used in the TPF

475

method. Table S3 contains the percentage of contributing configurations in the OSP

476

calculations for all amino acid sidechain analogues (reference states R1 and R2) and for all

477

amino acids (reference states R3, R4, R5 and R6). Table S4 compares the computational

478

efforts of TI and the TPF+OSP approach. Figure S1 shows the convergence of the energy

479

fluctuations as a function of time for selected amino acid side chain analogues. Building

480

blocks of all reference states for use in the GROMOS 54a8 force field.

481

This information is available free of charge at http://pubs.acs.org

482



21 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

483

Figures

484

485 486

Figure 1: The solvation free energy is calculated by turning the electrostatic and Van der

487

Waals interactions of the amino acid side chain analogues off. This is done in two steps: 1)

488

the uncharging free energy (Q>N) is computed by TI and TPF and 2) the decoupling free

489

energy (N>D), is computed by TI and OSP, using R1 or R2 as reference states.

490



22 ACS Paragon Plus Environment

Page 22 of 41

Page 23 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

491 492

Figure 2: Chemical structure of the reference states used for the OSP approach. Atom A

493

represents a soft-core particle as defined in.18

494



23 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

495 496

Figure 3: The free energy of mutation is computed using TI and the TPF and OSP approaches.

497

In the OSP approach, multiple sidechain conformations are considered, which are

498

subsequently averaged appropriately (see methods section).

499



24 ACS Paragon Plus Environment

Page 24 of 41

Page 25 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

CMS1

CMS2

CMS3

CMS4

CMS5

All

CMSs

500

Figure 4: Overlay of central member structures of arginine sidechains in the tripeptide. The

501

backbone of the first CMS of the tripeptide Ala – Arg – Ala shown in orange and five central

502

member structures of the sidechain of Arg in differently colored sticks.

503 504



25 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

505 506

Figure 5: Example of the uncharging free energy computed by TI (black curve), LRA (blue

507

curve) and TPF (red curve) for the sidechain analogue of Met.

508



26 ACS Paragon Plus Environment

Page 26 of 41

Page 27 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

509 !"#

510

Figure 6: Comparison of experimental (∆𝐺!"#$ ) and calculated (∆𝐺!!! ) solvation free

511

energies of amino acid sidechain analogues that do not carry a net charge. Calculations using

512

TI are compared to using the TPF+OSP approach with reference states R1 and R2. Outliers

513

Gln, Met, Trp and Tyr are not shown for the calculations using R1 (see Table 3).

514 515



27 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

516

517

Figure 7: Comparison of the free energies of mutation of amino acids into alanine in the

518

!" ) and TPF+OSP approach with all reference states tripeptide calculated via TI (∆G!!!

519

!"#!!"# (∆∆G!!! ). The right panel is a zoom of the indicated square in the left panel.

520



28 ACS Paragon Plus Environment

Page 28 of 41

Page 29 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

521

Journal of Chemical Theory and Computation

Tables

522 ∆𝑮𝑻𝑰 𝑸!𝑵

a.a.

∆𝑮𝑳𝑹𝑨 𝑸!𝑵

∆∆𝑮𝑳𝑹𝑨!𝑻𝑰 𝑸!𝑵

∆∆𝑮𝑻𝑷𝑭!𝑻𝑰 𝑸!𝑵

arg

138.9

±

0.3

152.6

±

0.1

13.8

144.0

±

0.2

5.1

asn

52.1

±

0.2

64.8

±

0.04

12.7

54.0

±

0.1

1.9

asp

338.6

±

0.4

359.2

±

0.1

20.6

339.5

±

0.3

0.9

cys

10.2

±

0.1

12.7

±

0.02

2.5

10.4

±

0.03

0.2

gln

51.8

±

0.2

64.9

±

0.04

13.1

54.3

±

0.1

2.5

glu

337.7

±

0.4

356.9

±

0.1

19.2

339.2

±

0.3

1.5

hisa

51.5

±

0.2

61.9

±

0.1

10.4

54.2

±

0.3

2.8

hisb

66.3

±

0.2

76.0

±

0.04

9.7

68.6

±

0.1

2.3

lys

178.4

±

0.3

192.0

±

0.1

13.7

187.3

±

0.3

8.9

met

14.7

±

0.1

17.4

±

0.02

2.8

14.9

±

0.03

0.2

phe

9.0

±

0.1

10.7

±

0.02

1.8

9.3

±

0.02

0.3

ser

29.5

±

0.1

40.4

±

0.03

10.8

32.5

±

0.1

3.0

thr

28.6

±

0.1

39.3

±

0.03

10.7

31.7

±

0.1

3.1

trp

32.1

±

0.1

35.9

±

0.03

3.8

32.7

±

0.1

0.6

tyr

32.5

±

0.1

43.0

±

0.03

10.5

35.9

±

0.1

3.4

RMSE

523 524 525 526 527

∆𝑮𝑻𝑷𝑭 𝑸!𝑵

11.8

3.3

Table 1: Calculated uncharging free energies of sidechain analogues in kJ/mol via TI, LRA !"#!!" !"#!!" and TPF and their absolute differences ( ∆∆𝐺!!! and ∆∆𝐺!!! ). Statistical error estimates are obtained from 1000 bootstrap replicates on the original data. The root mean square error (RMSE) is computed between TI and LRA and between TI and TPF.



29 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

R1 ∆𝑮𝑻𝑰 𝑵!𝑫

a.a. ala

-9.1

arg

-15.0

asn

± 0.3

Page 30 of 41

R2 ∆𝑮𝑶𝑺𝑷 𝑵!𝑫 -9.3

∆∆𝑮𝑶𝑺𝑷!𝑻𝑰 𝑵!𝑫

± 0.1

± 0.8 -1408

± 3

-9.5

± 0.5

-9.8

± 0.4

asp

-11.0

± 0.6

-10.9

cys

-2.9

± 0.4

-5.6

gln

-11.1

± 0.6 -1325

± 3

glu

-12.2

± 0.6

-9.4

± 0.7

hisA

-7.5

± 0.7

-11.8

hisB

-7.2

± 0.6

ile

-9.2

leu

∆𝑮𝑶𝑺𝑷 𝑵!𝑫

0.2

-

-21.2 ± 0.3

6.2

0.3

-20.7 ± 0.2

11.2

± 0.4

0.1

-19.0 ± 0.3

8.1

± 0.4

2.7

-21.3 ± 0.2

18.4

-16.2 ± 0.2

5.1

2.8

-13.2 ± 0.2

1.1

± 1.2

4.3

-13.1 ± 0.2

5.6

-9.2

± 1.3

2.0

-14.0 ± 0.2

6.8

± 0.6

-6.5

± 1.4

2.7

-9.4 ± 0.2

0.2

-10.3

± 0.7

-10.7

± 1.1

0.4

-12.5 ± 0.3

2.2

lys

-12.4

± 0.6

-16.9

± 2.0

4.5

-16.8 ± 0.3

4.4

met

-6.6

± 0.6

-7.2 ± 0.2

0.5

phe

-7.6

± 0.8

-12.7

± 2.4

5.1

-10.0 ± 0.3

2.4

ser

-5.4

± 0.3

-7.1

± 0.4

1.7

-25.4 ± 0.25

20.0

thr

-6.3

± 0.4

-8.6

± 0.3

2.3

-10.5 ± 0.57

4.2

trp

-7.0

± 1.0

-149

± 3

142

-1.3 ± 0.51

5.7

tyr

-6.9

± 0.9 -8583

± 3

8576

-6.5 ± 0.62

0.5

val

-8.5

± 0.5

± 0.4

-15.0 ± 0.24

6.5

RMSE 528 529 530 531

-438

-9.3

± 3

1393

-

∆∆𝑮𝑶𝑺𝑷!𝑻𝑰 𝑵!𝑫

1314

431

0.8 2018

8.1

Table 2: Decoupling free energies (in kJ/mol) of the neutral amino acid sidechain analogues !"# !"# !" via TI and OSP (∆𝐺!!! = −∆𝐺!!! + ∆𝐺!!! ) using reference states R1 and R2, together !"#!!" with their absolute differences (∆∆𝐺!!! ).



30 ACS Paragon Plus Environment

Page 31 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

R1 𝒆𝒙𝒑

R2

∆𝑮𝑻𝑰 𝑸!𝑫

∆𝑮𝑻𝑷𝑭!𝑶𝑺𝑷 𝑸!𝑫

∆𝑮𝑻𝑷𝑭!𝑶𝑺𝑷 𝑸!𝑫

∆∆𝑮𝑻!𝑶!𝑻𝑰 𝑸!𝑫

8.1

9.1

9.3

0.2

-

-

Ile

9.0

9.2

6.5

2.7

9.3

0.2

Leu

9.5

10.3

10.7

0.4

12.5

2.2

Val

8.3

8.5

9.3

0.8

14.9

6.4

Arg

-

-123.9

-122.8

1.1

Asn

-40.5

-42.6

-44.2

1.6

-33.3

9.2

Asp

-

-327.6

-328.6

1.0

-320.5

7.1

Cys

-5.2

-7.3

-4.8

2.5

10.8

18.1

Gln

-39.3

-40.6

-38.1

2.5

Glu

-

-325.5

-329.8

4.3

-326.0

0.5

HisA

-43.0

-44.0

-42.4

1.6

-41.2

2.8

HisB

-

-59.1

-59.5

0.4

-54.6

4.5

Lys

-

-166.0

-170.4

4.4

-170.5

4.5

Met

-6.2

-8.0

-7.8

0.3

Phe

-3.2

-1.4

3.4

4.8

0.7

2.1

Ser

-21.2

-24.1

-25.4

1.3

-7.1

17.0

Thr

-20.4

-22.3

-23.2

0.9

-21.2

1.1

Trp

-24.6

-25.1

116

142

-31.4

6.4

Tyr

-25.6

-25.5

8547

8572

-29.5

3.9

a.a.

∆𝑮𝒔𝒐𝒍𝒗

Ala

∆∆𝑮𝑻!𝑶!𝑻𝑰 𝑸!𝑫

1264

1388

1271

1312

423

432

RMSE

532 533 534 535

2018

7.2 !"#

Table 3: Solvation free energies of amino acid sidechains in kJ/mol. ∆𝐺!"#$ are !" experimentally measured solvation energies,61 ∆𝐺!!! obtained from thermodynamic !"#!!"# !!!!!" integration and ∆𝐺!!! obtained with the combined TPF+OSP approach. ∆∆𝐺!!! are !"#!!"# !" the absolute differences between ∆𝐺!!! and ∆𝐺!!! .

536 537



31 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

∆𝑮𝑻𝑰 𝑸!𝑵

a.a.

∆𝑮𝑳𝑹𝑨 𝑸!𝑵

∆∆𝑮𝑳𝑹𝑨!𝑻𝑰 𝑸!𝑵

∆𝑮𝑻𝑷𝑭 𝑸!𝑵

∆∆𝑮𝑻𝑷𝑭!𝑻𝑰 𝑸!𝑵

arg

217.0

±

0.1

230.8

±

0.1

13.8

222.5

±

0.3

5.5

asn

252.0

±

0.1

263.8

±

0.1

11.9

253.9

±

0.1

1.9

asp

415.9

±

0.2

437.3

±

0.1

21.5

416.9

±

0.5

1.1

cys

18.3

±

0.03

21.0

±

0.03

2.7

18.7

±

0.04

0.4

gln

268.2

±

0.1

284.6

±

0.1

16.5

274.3

±

0.1

6.2

glu

406.1

±

0.2

425.8

±

0.1

19.7

407.1

±

0.5

1.0

hisa

29.3

±

0.1

40.0

±

0.1

10.7

32.2

±

0.1

2.8

hisb

60.9

±

0.1

71.8

±

0.1

10.9

64.0

±

0.1

3.1

lysh

250.4

±

0.1

264.4

±

0.1

14.0

258.5

±

0.4

8.1

met

12.6

±

0.04

15.1

±

0.03

2.6

12.7

±

0.04

0.1

phe

0.6

±

0.03

2.4

±

0.02

1.7

1.0

±

0.03

0.3

ser

44.4

±

0.1

53.5

±

0.04

9.1

46.2

±

0.1

1.7

thr

44.0

±

0.1

51.9

±

0.04

7.9

45.4

±

0.1

1.4

trp

53.0

±

0.1

57.6

±

0.04

4.6

54.1

±

0.1

1.1

tyr

97.9

±

0.1

108.0

±

0.05

10.1

100.8

±

0.1

2.9

RMSE

538 539 540 541 542

Page 32 of 41

12.0

3.4

Table 4: Calculated uncharging free energies of sidechains in the tripeptides in kJ/mol via TI, !"#!!" !"#!!" LRA and TPF and their absolute differences (∆∆𝐺!!! and ∆∆𝐺!!! ). Statistical error estimates are obtained from 1000 bootstrap replicates on the original data. The root mean square error (RMSE) is computed between TI and LRA and between TI and TPF.



32 ACS Paragon Plus Environment

Page 33 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

∆∆𝑮𝑶𝑺𝑷!𝑻𝑰 [kJ/mol] 𝑵!𝑨 tripeptides

543 544 545 546

R3

R4

R5

R6

arg

2.2

4.3

1.8

2.7

asn

2.8

2.3

1.2

1.1

asp

6.2

3.6

0.8

1.3

cys

1.7

1.5

2.6

2.3

gln

1.6

4.8

0.3

0.1

glu

1.8

2.8

1.5

0.7

hisa

4.3

4.4

0.0

0.6

hisb

1.0

0.3

2.3

1.3

ile*

6.8

8.6

3.1

6.6

leu

3.5

7.4

0.3

0.3

lysh

4.3

3.8

1.7

3.6

met

7.0

6.0

0.6

1.3

phe

6.3

4.7

1.2

1.1

ser

2.6

3.3

2.8

2.5

thr

4.8

3.2

2.1

1.8

trp

16.0

6.3

5.7

8.0

tyr

4.5

2.8

0.3

1.6

val

2.7

4.5

2.2

2.4

RMSE

5.6

4.6

2.2

3.0

Table 5: Absolute difference between the free energy of annihilation calculated via TI and !"#!!" !"# !" OSP, ∆∆𝐺!!! = ∆∆𝐺!!! − ∆𝐺!!! , for the different reference states. *) Eight CMSs were used instead of five.



33 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

𝑻𝑷𝑭!𝑶𝑺𝑷 !𝑻𝑰

∆∆𝑮𝑸!𝑨

547 548 549 550

tripeptides

R3

Ile*

6.8

Leu

R4

[kJ/mol] R5

R6

8.6

3.1

6.6

3.5

7.4

0.3

0.3

Val

2.8

4.5

2.2

2.4

Arg

3.2

9.7

3.6

2.7

Asn

4.7

4.2

0.7

0.8

Asp

7.1

4.5

1.7

2.2

Cys

1.3

1.1

2.2

1.9

Gln

7.9

11.2

6.0

6.4

Glu

2.8

3.8

0.5

0.3

Hisa

4.6

4.7

0.3

0.9

Hisb

5.0

4.3

1.6

2.7

Lys

12.5

12.0

9.9

11.8

Met

7.1

6.1

0.7

1.4

Phe

6.7

5.1

0.8

1.4

Ser

0.8

1.5

1.0

0.7

Thr

3.4

1.9

0.8

0.4

Trp

14.9

5.2

6.8

6.9

Tyr

7.9

6.2

3.7

1.8

RMSE

6.7

6.4

3.6

4.1

Table 6: Absolute difference between the mutational free energy calculated via TI and OSP, (!"#!!"#)!!" !"#!!"# !" ∆∆𝐺!!! = ∆∆𝐺!!! − ∆𝐺!!! , for the different reference states. *) 8 CMSs were used instead of 5.

551



34 ACS Paragon Plus Environment

Page 34 of 41

Page 35 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

552

References

553 554

1.

Zhang, X. F.; Yang, G. Y.; Zhang, Y.; Xie, Y.; Withers, S. G.; Feng, Y., A general and efficient strategy for generating the stable enzymes. Sci. Rep. 2016, 6, 12.

555 556 557

2.

Won, A.; Pripotnev, S.; Ruscito, A.; Ianoul, A., Effect of Point Mutations on the Secondary Structure and Membrane Interaction of Antimicrobial Peptide Anoplin. J. Phys. Chem. B 2011, 115 (10), 2371-2379.

558 559 560

3.

Aceved, J. P.; Reetz, M. T.; Asenjo, J. A.; Parra, L. P., One-step combined focused epPCR and saturation mutagenesis for thermostability evolution of a new cold-active xylanase. Enzyme Microb. Technol. 2017, 100, 60-70.

561 562 563 564

4.

Moutaoufik, M. T.; Morrow, G.; Maaroufi, H.; Ferard, C.; Finet, S.; Tanguay, R. M., Oligomerization and chaperone-like activity of Drosophila melanogaster small heat shock protein DmHsp27 and three arginine mutants in the alpha-crystallin domain. Cell Stress Chaperones 2017, 22 (4), 455-466.

565 566 567

5.

Cortina, G. A.; Kasson, P. M., Excess positional mutual information predicts both local and allosteric mutations affecting beta lactamase drug resistance. Bioinformatics 2016, 32 (22), 3420-3427.

568 569 570

6.

Budzko, L.; Jackowiak, P.; Kamel, K.; Sarzynska, J.; Bujnicki, J. M.; Figlerowicz, M., Mutations in human AID differentially affect its ability to deaminate cytidine and 5methylcytidine in ssDNA substrates in vitro. Sci. Rep. 2017, 7 (1), 3873.

571 572

7. Bogan, A. A.; Thorn, K. S., Anatomy of hot spots in protein interfaces. J. Mol. Biol. 1998, 280 (1), 1-9.

573 574 575

8.

Du, X.; Li, Y.; Xia, Y. L.; Ai, S. M.; Liang, J.; Sang, P.; Ji, X. L.; Liu, S. Q., Insights into Protein-Ligand Interactions: Mechanisms, Models, and Methods. Int. J. Mol. Sci. 2016, 17 (2).

576 577 578

9.

Sharir-Ivry, A.; Varatharaj, R.; Shurki, A., Challenges within the Linear Response Approximation When Studying Enzyme Catalysis and Effects of Mutations. J. Chem. Theory. Comput. 2015, 11 (1), 293-302.

579 580 581

10. Almlof, M.; Aqvist, J.; Smalas, A. O.; Brandsdal, B. O., Probing the effect of point mutations at protein-protein interfaces with free energy calculations. Biophys. J. 2006, 90 (2), 433-442.

582 583 584 585

11. Kollman, P. A.; Massova, I.; Reyes, C.; Kuhn, B.; Huo, S. H.; Chong, L.; Lee, M.; Lee, T.; Duan, Y.; Wang, W.; Donini, O.; Cieplak, P.; Srinivasan, J.; Case, D. A.; Cheatham, T. E., Calculating structures and free energies of complex molecules: Combining molecular mechanics and continuum models. Acc. Chem. Res. 2000, 33 (12), 889-897.

586 587

12. Honig, B.; Nicholls, A., Classical electrostatic in biology and chemistry. Science 1995, 268 (5214), 1144-1149.

588 589 590

13. Onufriev, A.; Bashford, D.; Case, D. A., Exploring protein native states and large-scale conformational changes with a modified generalized born model. Proteins. 2004, 55 (2), 383-394.

35 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

591 592 593

14. Martins, S. A.; Perez, M. A. S.; Moreira, I. S.; Sousa, S. F.; Ramos, M. J.; Fernandes, P. A., Computational Alanine Scanning Mutagenesis: MM-PBSA vs TI. J. Chem. Theory. Comput. 2013, 9 (3), 1311-1319.

594 595 596 597

15. Simoes, I. C. M.; Costa, I. P. D.; Coimbra, J. T. S.; Ramos, M. J.; Fernandes, P. A., New Parameters for Higher Accuracy in the Computation of Binding Free Energy Differences upon Alanine Scanning Mutagenesis on Protein-Protein Interfaces. J. Chem. Inf. Model. 2017, 57 (1), 60-72.

598 599 600

16. Chodera, J. D.; Mobley, D. L.; Shirts, M. R.; Dixon, R. W.; Branson, K.; Pande, V. S., Alchemical free energy methods for drug discovery: progress and challenges. Curr. Opin. Struct. Biol. 2011, 21 (2), 150-160.

601 602

17. Steinbrecher, T.; Mobley, D. L.; Case, D. A., Nonlinear scaling schemes for LennardJones interactions in free energy calculations. J. Chem. Phys. 2007, 127 (21).

603 604

18. Schafer, H.; Van Gunsteren, W. F.; Mark, A. E., Estimating relative free energies from a single ensemble: Hydration free energies. J.Comput. Chem. 1999, 20 (15), 1604-1617.

605 606 607

19. Steinbrecher, T.; Joung, I.; Case, D. A., Soft-Core Potentials in Thermodynamic Integration: Comparing One- and Two-Step Transformations. J. Comput. Chem. 2011, 32 (15), 3253-3263.

608 609

20. Kirkwood, J. G., Statistical mechanics of fluid mixtures. Journal of Chemical Physics 1935, 3 (5), 300-313.

610 611

21. Liu, Y. L.; Patel, A. H. G.; Burger, S. K.; Ayers, P. W., Benchmarking pK(a) prediction methods for Lys115 in acetoacetate decarboxylase. J. Mol. Model. 2017, 23 (5), 8.

612 613 614

22. Maurer, M.; de Beer, S. B. A.; Oostenbrink, C., Calculation of Relative Binding Free Energy in the Water-Filled Active Site of Oligopeptide-Binding Protein. Molecules 2016, 21 (4).

615 616 617

23. Bekker, G. J.; Kamiya, N.; Araki, M.; Fukuda, I.; Okuno, Y.; Nakamura, H., Accurate Prediction of Complex Structure and Affinity for a Flexible Protein Receptor and Its Inhibitor. J. Chem. Theory Comput. 2017, 13 (6), 2389-2399.

618 619 620

24. Tanida, Y.; Ito, M. S.; Fujitani, H., Calculation of absolute free energy of binding for theophylline and its analogs to RNA aptamer using nonequilibrium work values. Chem. Phys. 2007, 337 (1-3), 135-143.

621 622 623

25. Beierlein, F. R.; Kneale, G. G.; Clark, T., Predicting the Effects of Basepair Mutations in DNA-Protein Complexes by Thermodynamic Integration. Biophys. J. 2011, 101 (5), 1130-1138.

624 625 626

26. Liu, H. Y.; Mark, A. E.; van Gunsteren, W. F., Estimating the relative free energy of different molecular states with respect to a single reference state. J. Phys. Chem. 1996, 100 (22), 9485-9494.

627 628 629

27. Oostenbrink, C.; van Gunsteren, W. F., Free energies of binding of polychlorinated biphenyls to the estrogen receptor from a single simulation. Proteins. 2004, 54 (2), 237246.



36 ACS Paragon Plus Environment

Page 36 of 41

Page 37 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

630 631

28. Oostenbrink, C.; van Gunsteren, W. F., Free energies of ligand binding for structurally diverse compounds. Proc. Natl. Acad. Sci. U.S.A. 2005, 102 (19), 6750-4.

632 633 634

29. Beutler, T. C.; Mark, A. E.; van Schaik, R. C.; Gerber, P. R.; van Gunsteren, W. F., Avoiding singularities and numerical instabilities in free energy calculations based on molecular simulations. Chem. Phys. Lett. 1994, 222 (6), 529-539.

635 636 637 638

30. Graf, M. M. H.; Lin, Z. X.; Bren, U.; Haltrich, D.; van Gunsteren, W. F.; Oostenbrink, C., Pyranose Dehydrogenase Ligand Promiscuity: A Generalized Approach to Simulate Monosaccharide Solvation, Binding, and Product Formation. PLoS Comput. Biol. 2014, 10 (12).

639 640 641

31. Ciordia, M.; Perez-Benito, L.; Delgado, F.; Trabanco, A. A.; Tresadern, G., Application of Free Energy Perturbation for the Design of BACE1 Inhibitors. J. Chem. Inf. Model. 2016, 56 (9), 1856-1871.

642 643 644

32. Raman, E. P.; Vanommeslaeghe, K.; MacKerell, A. D., Site-Specific Fragment Identification Guided by Single-Step Free Energy Perturbation Calculations. J. Chem.Theory Comput. 2012, 8 (10), 3513-3525.

645 646 647

33. Cave-Ayland, C.; Skylaris, C. K.; Essex, J. W., Direct Validation of the Single Step Classical to Quantum Free Energy Perturbation. J. Phys. Chem. B 2015, 119 (3), 10171025.

648 649

34. Loeffler, H. H.; Michel, J.; Woods, C., FESetup: Automating Setup for Alchemical Free Energy Simulations. J. Chem. Inf. Model. 2015, 55 (12), 2485-2490.

650 651 652

35. Gapsys, V.; Michielssens, S.; Seeliger, D.; de Groot, B. L., pmx: Automated Protein Structure and Topology Generation for Alchemical Perturbations. J. Comput. Chem. 2015, 36 (5), 348-354.

653 654

36. Chiang, Y. C.; Wang, Y., Virtual Substitution Scan Via Single-Step Free Energy Perturbation. Biopolymers 2016, 105 (6), 324-336.

655 656 657

37. Oostenbrink, C.; Van Gunsteren, W. F., Single-step perturbations to calculate free energy differences from unphysical reference states: Limits on size, flexibility, and character. J. Comput. Chem. 2003, 24 (14), 1730-1739.

658 659

38. Pitera, J. W.; van Gunsteren, W. F., One-step perturbation methods for solvation free energies of polar solutes. J. Phys. Chem. B 2001, 105 (45), 11264-11274.

660 661 662

39. Oostenbrink, C., Efficient free energy calculations on small molecule host-guest systems - a combined linear interaction energy/one-step perturbation approach. J. Comput. Chem. 2009, 30 (2), 212-21.

663 664 665

40. Sham, Y. Y.; Chu, Z. T.; Tao, H.; Warshel, A., Examining methods for calculations of binding free energies: LRA, LIE, PDLD-LRA, and PDLD/S-LRA calculations of ligands binding to an HIV protease. Proteins. 2000, 39 (4), 393-407.

666 667

41. Aqvist, J.; Medina, C.; Samuelsson, J. E., A new method for predicting binding affinity in computer-aided drug design. Protein Eng. 1994, 7 (3), 385-91.

668 669

42. Aqvist, J.; Hansson, T., On the validity of electrostatic linear response in polar solvents. J. Phys. Chem. 1996, 100 (22), 9512-9521.

37 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

670 671

43. de Ruiter, A.; Oostenbrink, C., Efficient and Accurate Free Energy Calculations on Trypsin Inhibitors. J. Chem. Theory Comput. 2012, 8 (10), 3686-3695.

672 673 674

44. Schmid, N.; Christ, C. D.; Christen, M.; Eichenberger, A. P.; van Gunsteren, W. F., Architecture, implementation and parallelisation of the GROMOS software for biomolecular simulation. Comput. Phys. Commun. 2012, 183 (4), 890-903.

675 676 677

45. Reif, M. M.; Winger, M.; Oostenbrink, C., Testing of the GROMOS Force-Field Parameter Set 54A8: Structural Properties of Electrolyte Solutions, Lipid Bilayers, and Proteins. J. Chem. Theory Comput. 2013, 9 (2), 1247-1264.

678 679

46. Chemical Computing Group ULC Molecular Operating Environment (MOE), 2013.08; 1010 Sherbooke St. West, Suite #910, Montreal, QC, Canada, H3A 2R7,, 2017.

680 681 682 683 684

47. Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J., Interaction Models for Water in Relation to Protein Hydration. In Intermolecular Forces: Proceedings of the Fourteenth Jerusalem Symposium on Quantum Chemistry and Biochemistry Held in Jerusalem, Israel, April 13–16, 1981, Pullman, B., Ed. Springer Netherlands: Dordrecht, 1981; pp 331-342.

685 686 687

48. Amadei, A.; Chillemi, G.; Ceruso, M. A.; Grottesi, A.; Di Nola, A., Molecular dynamics simulations with constrained roto-translational motions: Theoretical basis and statistical mechanical consistency. J. Chem. Phys. 2000, 112 (1), 9-23.

688 689 690

49. Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Dinola, A.; Haak, J. R., Molecular dynamics with coupling to an external bath. Journal of Chemical Physics 1984, 81 (8), 3684-3690.

691 692 693

50. Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C., Numerical integration of cartesian equations of motion of a system with constraints - molecular dynamics of N-alkanes. J. Comput. Phys. 1977, 23 (3), 327-341.

694 695

51. Tironi, I. G.; Sperb, R.; Smith, P. E.; van Gunsteren, W. F., A generalized reaction field method for molecular dynamics simulations. J. Chem. Phys. 1995, 102 (13), 5451-5459.

696 697 698

52. Heinz, T. N.; van Gunsteren, W. F.; Hünenberger, P. H., Comparison of four methods to compute the dielectric permittivity of liquids from molecular dynamics simulations. J. Chem. Phys. 2001, 115 (3), 1125-1136.

699 700

53. Daura, X.; van Gunsteren, W. F.; Mark, A. E., Folding-unfolding thermodynamics of a beta-heptapeptide from equilibrium simulations. Proteins 1999, 34 (3), 269-80.

701 702 703

54. Lee, F. S.; Chu, Z. T.; Bolger, M. B.; Warshel, A., Calculations of antibody-antigen interactions: microscopic and semi-microscopic evaluation of the free energies of binding of phosphorylcholine analogs to McPC603. Protein Eng 1992, 5 (3), 215-28.

704 705 706

55. Sham, Y. Y.; Chu, Z. T.; Tao, H.; Warshel, A., Examining methods for calculations of binding free energies: LRA, LIE, PDLD-LRA, and PDLD/S-LRA calculations of ligands binding to an HIV protease. Proteins 2000, 39 (4), 393-407.

707 708

56. Warshel, A.; Russell, S. T., Calculations of electrostatic interactions in biological systems and in solutions. Q Rev Biophys 1984, 17 (3), 283-422.



38 ACS Paragon Plus Environment

Page 38 of 41

Page 39 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

709 710

57. Chipot, C., Pohorille, A., Free energy calculations: Theory and applications in chemistry and biology . Springer: Berlin, 2007.

711 712

58. Zwanzig, R. W., High‐Temperature Equation of State by a Perturbation Method. I. Nonpolar Gases. J. Chem. Phys. 1954, 22 (8), 1420-1426.

713 714

59. Pitera, J. W.; van Gunsteren, W. F., One-step perturbation methods for solvation free energies of polar solutes. Journal of Physical Chemistry B 2001, 105 (45), 11264-11274.

715 716 717

60. Mark, A. E.; Xu, Y. W.; Liu, H. Y.; van Gunsteren, W. F., Rapid non-empirical approaches for estimating relative binding free energies. Acta Biochim. Pol. 1995, 42 (4), 525-535.

718 719

61. Wolfenden, R.; Andersson, L.; Cullis, P. M.; Southgate, C. C. B., Affinities of aminoacid side-chains for solvent water. Biochemistry 1981, 20 (4), 849-855.

720 721 722 723

62. Reif, M. M.; Hünenberger, P. H., Computation of methodology-independent single-ion solvation properties from molecular simulations. IV. Optimized Lennard-Jones interaction parameter sets for the alkali and halide ions in water. J. Chem. Phys. 2011, 134 (14).

724 725 726

63. Reif, M. M.; Hünenberger, P. H.; Oostenbrink, C., New Interaction Parameters for Charged Amino Acid Side Chains in the GROMOS Force Field. J. Chem. Theory Comput. 2012, 8 (10), 3705-3723.

727 728



39 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

729

Table of contents graphic

730



40 ACS Paragon Plus Environment

Page 40 of 41

Page 41 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

table of contents graphic 145x64mm (300 x 300 DPI)

ACS Paragon Plus Environment