Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Open Access

Peer-reviewed

Research Article

Assigning crystallographic electron densities with free energy calculations—The case of the fluoride channel Fluc

Abstract

Approximately 90% of the structures in the Protein Data Bank (PDB) were obtained by X-ray crystallography or electron microscopy. Whereas the overall quality of structure is considered high, thanks to a wide range of tools for structure validation, uncertainties may arise from density maps of small molecules, such as organic ligands, ions or water, which are non-covalently bound to the biomolecules. Even with some experience and chemical intuition, the assignment of such disconnected electron densities is often far from obvious. In this study, we suggest the use of molecular dynamics (MD) simulations and free energy calculations, which are well-established computational methods, to aid in the assignment of ambiguous disconnected electron densities. Specifically, estimates of (i) relative binding affinities, for instance between an ion and water, (ii) absolute binding free energies, i.e., free energies for transferring a solute from bulk solvent to a binding site, and (iii) stability assessments during equilibrium simulations may reveal the most plausible assignments. We illustrate this strategy using the crystal structure of the fluoride specific channel (Fluc), which contains five disconnected electron densities previously interpreted as four fluoride and one sodium ion. The simulations support the assignment of the sodium ion. In contrast, calculations of relative and absolute binding free energies as well as stability assessments during free MD simulations suggest that four of the densities represent water molecules instead of fluoride. The assignment of water is compatible with the loss of these densities in the non-conductive F82I/F85I mutant of Fluc. We critically discuss the role of the ion force fields for the calculations presented here. Overall, these findings indicate that MD simulations and free energy calculations are helpful tools for modeling water and ions into crystallographic density maps.

Citation: Ariz-Extreme I, Hub JS (2018) Assigning crystallographic electron densities with free energy calculations—The case of the fluoride channel Fluc. PLoS ONE 13(5): e0196751. https://doi.org/10.1371/journal.pone.0196751

Editor: Freddie Salsbury Jr, Wake Forest University, UNITED STATES

Received: January 12, 2018; Accepted: April 18, 2018; Published: May 17, 2018

Copyright: © 2018 Ariz-Extreme, Hub. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability: All relevant data are within the paper and its Supporting Information files.

Funding: This study was supported by the Deutsche Forschungsgemeinschaft (grant numbers HU 1971-1/1, HU 1971-3/1) to JSH (URL: www.dfg.de). IA was supported by the Göttingen Graduate School for Neurosciences, Biophysics, and Molecular Biosciences (GGNB). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

Introduction

Roughly 130000 macromolecular structures have been deposited in the Protein Data Bank (PDB) by 2017 [1, 2], of which more than 120000 have been determined using X-ray crystallography or cryo-electron microscopy. Overall, the quality of macromolecular structures in the PDB is very high, in part thanks to protein modeling and validation routines such as COOT, MolProbity, Procheck, NQ-Flipper, WHAT_CHECK, and several others [3–7]. Besides, the project PDB_REDO provides a protein re-refinement procedure in which the entire content of the PDB is re-evaluated [8]. However, whereas crystal structure analysis has been partially automated, the interpretation of certain electron density maps still needs a significant amount of human supervision and expertise [9], specially for non-covalently bound molecules such as organic ligands, water, and ions, which are usually more disordered than nearby residues [10, 11]. For certain ions, the anomalous scattering signal may guide the correct assignment. However, because the anomalous signal is often weak, this approach is not always applicable.

For instance, since the structure of K+ channel KcsA from Streptomyces lividans was obtained in 1998 (PDB ID 1BL8) [12], the interpretation of the electron densities found in the selectivity filter has triggered lively debates. Whether the electron densities represent two cations intercalated with water molecules [13, 14] or four cations in a row [15] has remained controversial. As another example of ambiguous electron densities may serve the obsolete models of β-glucosidase that contained 252 and 199 sodium ions, respectively, later mostly reassigned to water molecules (PDB IDs 3FJ0 and 3FIY, superseded by 4HZ8 and 4HZ6, respectively) [16]. Further, in the model of cyclic bovine pancreatic trypsin inhibitor (PDB ID 1K6U), a single water molecule should have been modeled instead of ethylene glycol molecule [17, 18]. Some tools such as PURY [19], CheckMyMetal [20], Twilight [21], or PHENIX [22] have been developed in an attempt to resolve the nature of isolated electron densities.

Molecular dynamics (MD) simulations are an established tool for refining atomic structures against experimental data, as for instance implemented in the CNS and MDFF/xMDFF software suites [23–26]. Further, Malde and Mark (2011) [27] showed that MD simulations and free energy calculations provide a useful framework to assign ligands in X-ray crystal complexes. For example, this proved to be a succesful strategy for determining the stereochemistry of the HIV-1 protease inhibitor JG-365 (PDB ID 7HVP) [28, 29] or the orientation of the ligand L-Serine in the binding pocket of Pyrococcus abyssi threonil-tRNA synthetase (PDB ID 2HKZ) [30]. Given the constant improvement of hardware and algorithms, MD simulations and free energy calculations are increasingly resorted techniques that may help to solve ambiguous assignments of electron densities. Because ions and water play key functional roles in many proteins, the correct assignment of such densities may aid the understanding of protein function at molecular level. With this aim in mind, we here investigate electron densities of the fluoride extruding Fluc channel.

Fluc is a family of small membrane proteins found in a wide range of unicellular organisms [31–33]. The expression of Fluc channels is controlled by a fluoride-responsive riboswitch, and their function is to allow the export of fluoride anions, thereby keeping the cytoplasmic F− concentration below the enzyme inhibitory levels [34, 35]. Electrophysiological evidence states that Fluc channels are highly selective for fluoride anions over chloride or sodium ions, with a time- and voltage-independent current, reaching ∼ 107 fluoride ions/sec at −200 mV [32]. However, the molecular details about the permeation mechanism remain unknown.

Fluc channels arrange in an antiparallel homodimeric double-barreled structure that is new among membrane proteins (PDB ID 5NKQ) [36]. The crystal structure of the Bordetella pertussis Fluc (Fluc-Bpe) is blocked by L2 monobodies (Fig 1A) [37, 38], and reveals five intriguing disconnected electron densities inside the protein (Fig 1B). The first of the electron densities is coordinated by four carbonyl groups of residues G77 and T80 at the break of the two transmembrane helices 3 (TM3 site). This density was assigned to a structural sodium ion although this type of tetrahedral coordination is only seen in ∼ 3.5% of sodium binding sites, according to the MetalPDB database [39]. The remaining four electron densities are surrounded by the electropositive side chains of N43, S108, S112, F82, and F85. Because this environment is compatible with a halide coordination shell, these four electron densities were “provisionally” identified as fluoride anions, as emphasized the authors of the 5NKQ structure [36]. Interestingly, Last et al. showed that crystal structures of the non-conductive F80I and F83I mutants of the Escherichia coli Fluc (PDB ID 5KBN and 5KOM, respectively), lack the electron densities associated to F80 and F83 residues (analogues of F82 and F85 residues in Fluc-Bpe, respectively) [40]. In other words, the loss of electron densities at the F80 and F83 sites correlates with a loss of protein function. This intriguing finding was taken as an indication that the electron densities at the F82 and F85 sites represent functionally relevant fluoride binding sites.

thumbnail
Download:
Fig 1. Fluc channel simulation system and structure.

(A) Fluc channel simulation system with blocking L2 monobodies. Fluc is depicted as orange and yellow ribbons, L2 monobodies are showed as pink thick ribbons, lipids heads are shown as purple volume and lipid tails as white sticks. Water molecules not shown for the sake of clarity. (B) Close-up of Fluc channel electron densities (cyan) at F82 (blue) and F85 (green) sites. The sodium ion at the TM3 site is also highlighted in yellow. (C) Fluc channel simulation system without blocking L2 monobodies.

https://doi.org/10.1371/journal.pone.0196751.g001

In this study, taking Fluc as an example, we illustrate how free MD simulations and free-energy calculations may help with the assignment of ambiguous disconnected electron densities. For the TM3 site, we find that variations in the structural stability during free simulation indeed favor sodium over water, in line with the original assignments. For the four densities at the F82 and F85 residues, however, free-energy calculations as well as stability assessments strongly favor water over fluoride.

Materials and methods

System setup for MD simulations

Amber ff99SB*-ILDN forcefield.

Initial coordinates of Fluc-Bpe were taken from the Protein Data Bank (PDB ID 5FXB, now superseded by 5NKQ) [36]. Initial coordinates for the lipid membrane were taken from a 1 ns equilibrium simulation of a hydrated membrane of 328 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) lipids. The protein was embedded in the membrane with the g_membed software [41]. The simulation box was filled by explicit TIP3P [42] water and neutralized with counter ions. Four different simulation systems were set up either of the wild-type (WT) or of the F82I/F85I mutant, each including or excluding the blocking L2 monobodies (which were present in the crystal). The final simulation system with blocking L2 monobodies contained one homodimer Fluc channel, two blocking L2 monobodies, 188 POPC molecules, 17685 water molecules, 6 sodium ions (including the structural Na+), and 4 fluoride ions (Fig 1A). Similarly, the final simulation system without L2 monobodies contained one homodimer Fluc channel 204 POPC molecules, 18052 water molecules, 1 structural sodium ion, and 7 fluoride ions (Fig 1C). The ion parameters were taken from Joung and Cheatham (2008) [43]. The F82I/F85I mutant was generated using PyMOL [44].

Interactions of the protein atoms were described by the Amber ff99SB*-ILDN forcefield [45], and lipid parameters were taken from Berger et al. and Cordomi et al. [46, 47]. The simulations were carried out with the GROMACS simulation software (version 5.1.4) [48]. Electrostatic interactions were calculated at every step with the particle-mesh Ewald method [49]. Short-range repulsive and attractive dispersion interactions were described by a Lennard-Jones (LJ) potential, with a cut-off at 1 nm. The geometry of water molecules was constrained with the SETTLE algorithm [50], and all other bond lengths were constrained with LINCS [51]. Hydrogen atoms of the protein were constructed as virtual sites, allowing a 4 fs time step [52].

The simulation temperature was controlled at 300 K using velocity rescaling (τ = 2.5 ps) [53], and the pressure was kept at 1 bar with a semiisotropic Berendsen barostat (τ = 2 ps) [54]. The systems were equilibrated for 20 ns with position restraints on the protein heavy atoms and on the molecules placed at the electron density sites (fluoride or water molecules at F82 site and F85 site, see Fig 1B), with a force constant of 1000 kJ mol−1nm−2. After that, the systems were equilibrated for 400 ns with position restraints (k = 1000 kJ mol−1nm−2) on the protein Cα atoms and on the molecules placed at F82 and F85 sites before production.

To investigate the stability of the fluoride and water ions assigned to the densities at the F82 and F85, multiple free production simulations were started from random time points of the 400 ns equilibration simulations. Here, four different system types were defined for free simulations: (i) four fluoride ions at F82 and F85 sites / without monobodies; (ii) fluoride ions at F82 and F85 sites / with monobodies; (iii) four water at F82 and F85 sites / without monobodies; (iv) water at F82 and F85 sites / with monobodies. For systems with or without blocking L2 monobodies, we carried out 10 replica simulations of 80 ns each or 20 replica simulations of 40 ns each, respectively. All the production simulations were carried out without position restraints, and using the parameters described above.

CHARMM36 forcefield.

To exclude that the qualitative results depend on the applied protein force field, another set of simulation systems were setup using the CHARMM force field family [55, 56]. Setup of the system was identical to the setup for the Amber force field (see above), except that Fluc was solvated in a membrane of 1,2-dipalmitoyl-sn-phosphatidilcholine (DPPC) instead of POPC. The final CHARMM simulation system contained one homodimer of WT Fluc, 228 DPPC molecules, 13273 TIP3P water molecules, 1 structural sodium ion, 3 chloride ions, and 4 fluoride ions. The fluoride ion parameters were again taken from Joung and Cheatham (2008) [43], and the sodium and chloride parameters correspond to the CHARMM36 defaults.

Interactions of the protein atoms and lipid atoms were described by the CHARMM36 forcefield [55, 56]. The simulations were carried out with the GROMACS simulation software (version 5.1.4) [48]. Electrostatic interactions were calculated at every step with the particle-mesh Ewald method [49], and using a cut-off at 1.2 nm for LJ interactions. The geometry of water molecules was constrained with the SETTLE algorithm [50], and all bonds of hydrogen atoms were constrained with LINCS [51]. A 2 fs timestep was used. The simulation temperature was controlled at 323 K using velocity rescaling (τ = 2.5 ps) [53], and the pressure was kept at 1 bar with a semiisotropic Berendsen barostat (τ = 2 ps) [54]. The system was equilibrated for 20 ns with position restraints on the protein heavy atoms and on the fluoride ions placed at the F82 and F85 sites with a force constant of 1000 kJ mol−1nm−2. Subsequently, the system was equilibrated for 5 ns with position restraints on the protein Cα atoms and on the fluoride ions placed at F82 and F85 sites before production.

Finally, by assigning new random velocities, we generated 10 independent replicas from the last snapshot taken from the 5 ns equilibration trajectory, and we ran each replica for 40 ns.

Polarizable DRUDE forcefield.

To furter exclude that polarizability could affect the qualitative results, a simulation system was set up using the polarizable DRUDE forcefield [57–59]. Here, initial coordinates for the lipid membrane were taken from a 3 ns equilibrium simulation of a hydrated membrane of 256 DPPC lipids. The protein was again embedded in the membrane with the g_membed software [41]. The simulation box was filled by explicit SWM4-NDP [60] water and neutralized with counter ions. The final simulation system contained one homodimer WT Fluc channel, 228 DPPC molecules, 10181 water molecules, 1 sodium ion, and 7 fluoride ions. The ion parameters were taken from Yu et al. (2010) [61].

Interactions of the protein atoms were described by the polarizable DRUDE forcefield [57–59], and lipid parameters were taken from Chowdhary et al. [62]. The GROMACS DPPC topology file was generated from original DRUDE-CHARMM topology files using TopoGromacs [63]. The simulations were carried out with the GROMACS simulation software (version DRUDE) [64, 65]. Electrostatic interactions were calculated at every step with the particle-mesh Ewald method [49], and using a cut-off at 1.2 nm for LJ interactions. The geometry of water molecules was constrained with the SETTLE algorithm [50], and all other bond lengths were left unconstrained.

To avoid polarization catastrophes, the system was carefully equilibrated with multi-step protocol. First, we equilibrated the simulation system using the self-consistent field (SCF) approach, as implemented by Van Maaren and Van der Spoel [64]. The simulation temperature was controlled at 323 K using velocity rescaling (τ = 0.1 ps) [53], and the pressure was kept at 1 bar with a semiisotropic Berendsen barostat (τ = 1 ps) [54]. Using a 0.2 fs time step, the system was equilibrated for 37 ps with position restraints on the protein heavy atoms (force constant of 1000 kJ mol−1nm−2). Next, the system was equilibrated for 800 ps with position restraints on the protein Cα atoms and on the fluoride ions placed at the F82 and F85 sites. In this last equilibration, the simulation temperature was controlled at 323 K for the non-drude atoms, and at 1 K for the drude particles, using the Nose-Hoover thermostat (τ = 0.1 ps, and τ = 0.005 ps, respectively) [66, 67]. No pressure-coupling was used, and a 1 fs time step was used with the extended Lagrangian dynamics implemented in Gromacs by Lemkul et al. [65].

Finally, we generated 10 different replicas from the last snapshot taken from the 800 ps equilibration trajectory and ran each of them for 10 ns.

Free-energy calculations of absolute and relative binding

We carried out two sets of free energy calculations: First, to compute relative binding free energies, an ion was alchemically transformed into a water molecule at the binding site while simultaneously transforming a water molecule into an ion in bulk solvent. According to a thermodynamic cycle (Fig 2), the associated free energy difference equals the free energy difference ΔΔGbind of binding an ion versus a water molecule into the binding site. A negative ΔΔGbind implies that a water molecule is more likely to translocate from solvent to the binding site as compared to the ion. More quantitatively, in equilibrium, the probabilities of finding a water molecule or of finding an ion at the binding site are related via , where cwat and cion denote the water and ion concentrations in the solvent, kB is the Boltzmann constant, and T the absolute temperature. Hence, the ΔΔGbind is in principle the key quantity to detect the correct assignment of the ambiguous electron densities. However, ΔΔGbind is associated with an increased uncertainty, mainly because force fields might not represent ion-protein interactions as accurately as water-protein interactions [68–72].

thumbnail
Download:
Fig 2. Alchemical transformation thermodynamic cycle.

From the thermodynamic cycle, the binding free energy difference ΔΔGbind = ΔG1 − ΔG4 between the ion-bound Fluc channel and the water-bound Fluc channel can be calculated via ΔG3 − ΔG2.

https://doi.org/10.1371/journal.pone.0196751.g002

Second, absolute binding free energies ΔGbind for water molecules were computed, that is, free energy for translocating a water molecule from solvent into the binding site. Because the probability of finding a water molecule at the binding site is proportional to exp(−ΔGbind/kBT), the assignment of a water molecule to an ambiguous electron density is only plausible if ΔGbind is negative or at least close to zero.

All free energy calculations were conducted using discrete thermodynamic integration (TI) along an alchemical reaction coordinate λ, where λ = 0 and λ = 1 denote the initial and the final state. During TI calculations, a time step of 2 fs was applied, and the temperature was controlled at 300 K using a stochastic dynamics integration scheme (τ = 0.1 ps) [73]. Starting configurations were taken from the final snapshot of the 20 ns equilibrium simulation. All the other simulation parameters were identical to the equilibrium simulations.

Relative binding free-energy computational details.

The parameter λ was chosen such that λ = 0 described the state with the ion in the protein (fluoride at F82 or at F85 site, or sodium TM3 site) while having a water in bulk solvent. λ = 1 described the ion (fluoride or sodium) in bulk solvent while having a water molecule at the F82, F85, or TM3 site. The free energy calculations with fluoride were done in the presence or absence of blocking L2 monobodies, while the free energy calculations with sodium were done only in presence of monobodies. To keep the structure close to the crystal structure, position restraining potentials (k = 1000 kJ mol−1nm−2) were applied on the protein Cα atoms and on the ion and water molecules that were alchemically transformed. Here, the molecule in solvent was restrained at a position as far as possible from the protein and membrane. To avoid numerical instabilities close to λ = 0 and λ = 1, soft-core potentials were used for both electrostatics and Lennard-Jones interactions with α = 0.5, σ = 0.3 and a soft-core power of 1. The transformation from λ = 0 to λ = 1 was done in 41 λ-steps, where each λ-steps was simulated for 10 ns. Derivatives of the Hamiltonian with respect to λ were recorded at every step. Free energy differences were calculated from (1) where 〈⋅〉 denote the average over the trajectory, after removing the first 50 ps for equilibration.

Absolute binding free energy computational details.

Absolute binding free energies ΔGbind of water were calculated using TI following the scheme shown in Fig 3. Accordingly, ΔGbind was calculated from the difference in free energy between (i) turning off the interactions of a water molecule in a pure water box (ΔG2 in Fig 3), and (ii) turning off the interactions of a water molecule in the WT Fluc channel at F82 and F85 sites, or in the mutant channel at F82I and F85I sites (ΔG1 in Fig 3). Since there are two F82 sites and two F85 sites in Fluc, we calculated the absolute ΔGbind of water to each of the two binding sites.

thumbnail
Download:
Fig 3. Absolute binding free energy thermodynamic cycle.

From the thermodynamic cycle, the absolute binding free energy is calculated via ΔGbind = ΔG2 − ΔG1, where ΔG3 ≡ 0.

https://doi.org/10.1371/journal.pone.0196751.g003

As described by Deng and Roux [74], the sampling is more efficient during calculations of absolute binding free energies if the molecule is restraint at the binding site. Therefore, to keep the water molecules at the positions indicated by the electron densities in the crystal structure, we applied position restraints to the oxygen atom of water molecules at Fluc binding sites. Likewise, we also applied position restraints on the water molecule in the pure water box. Following the scheme shown in Fig 3, the calculation of ΔG1 was separated in five sub-steps, namely, (a) , turning on the position restraints with a force constant Kx = 200 kJ mol−1nm−2; was typically in the order of only 1 kJ mol−1, (b) , turning off electrostatic, and (c) , turning off Lennard-Jones interactions, (d) translating the non-interacting solute from the protein into a pure-solvent box, and (e) release of the position restraints in pure solvent. To obtain , 11 λ-steps were simulated for 2 ns each for all systems.

was computed from 11 λ-steps simulated for 2 ns each, and was computed in 41 λ-steps simulated for 2 ns each. To account for atomic overlaps occurring close to λ = 0 and λ = 1, soft-core potentials were used for Lennard-Jones interactions with α = 0.5, σ = 0.3 and a soft-core power of 1. As before, the first 50 ps of each trajectory were removed for equilibration, and the free energy differences were calculated according to Eq 1. The calculation of ΔG2 in Fig 3 was separated into the electrostatic part and Lennard-Jones interactions part, and the free energy due to applying position restraints on a non-interacting water particle was calculated analytically following ref. [75]. Accordingly, applying a translational restraining potential Ur,x with a force constant Kx and the form Ur,x = Kx(x − x0)2/2 confines the molecule to a small effective volume Veff = ∫V exp[−Kx(x − x0)2/2kBT]dx that can be analytically evaluated as . Hence, the corresponding free energy can be calculated via where V∘ = 0.02992 nm3 is the molecular volume of water [76]. We used Kx = 200 kJ mol−1nm−2, yielding 0.774225 kJ mol−1.

Maximum-likelihood estimates for the lifetime of fluoride and water molecules inside Fluc

To calculate the lifetime τ of the fluoride or water molecules at F82 and F85 positions inside Fluc, we used the maximum-likelihood estimate as described recently [77]. In summary, given M simulations of simulation time Ts, during which m fluoride or water leave the binding site at times ti (i = 1, …, m; m ≤ M), then the maximum-likelihood estimate τml for the lifetime of fluoride or water at the binding site is given via .

Quantum chemical calculations

The relaxed potential energy surface (PES) calculations of fluoride and benzene were carried out using density functional theory with the B3LYP functional [78], using either the ma-def2-SVP or the ma-def2-TZVP basis set [79, 80]. One additional scan was conducted using Møller-Plesset perturbation theory of second order (MP2) using the aug-cc-pVDZ basis set [81]. PES scans were computed with the ORCA 3 program system [82].

Results

Fluoride is less stable than water inside Fluc in free MD simulations

The authors of the 5NKQ structure provisionally assigned fluoride ions F− to four disconnected electron densities found in Fluc at F82 and F85 sites (Fig 1B) [36]. To test the stability of F− inside Fluc in crystallographic and physiological conditions (i.e. in the presence and absence of monobodies), we carried out MD simulations with the systems depicted in Fig 1A and 1C. To obtain sufficient statistics, we carried out 30 replicas with fluoride starting at F82 and F85, among which in 10 replicas (80 ns each) Fluc was blocked by L2 monobodies (Fig 1A), while in 20 replicas (40 ns each) Fluc was free of monobodies (Fig 1C). For comparison, we performed similar simulations with water molecules starting at F82 and F85 sites. Typical time traces of the positions of water and fluoride along the channel z-axis are shown in Fig 4. Evidently, fluoride escaped earlier and more often than water in the simulations. The time traces of all the simulations are shown in S1, S2, S3 and S4 Figs. To exclude that the choice for the forcefield may influence the qualitative results, we simulated another 10 replicas of the Fluc system without monobodies using the CHARMM36 forcefield, and another 10 replicas with the polarizable DRUDE forcefield. The time traces of these simulations are shown in S5 and S6 Figs. Overall, the time traces reveal that water is more stable in the channel than fluoride.

thumbnail
Download:
Fig 4. Fluoride and water molecules unbinding from the F82 and F85 sites of Fluc.

Position change in z-direction normal to membrane in MD simulations. (A) Fluoride position inside Fluc with blocking L2 monobodies. (B) Water molecule position in Fluc with blocking L2 monobodies. (C) Fluoride position inside Fluc without monobodies. (D) Water molecule position in Fluc without monobodies. Green curves represent the molecules starting at F85 site, and blue curves belong to molecules starting at F82 site. Different shades of green and blue indicate starting positions in the two homodimers of Fluc.

https://doi.org/10.1371/journal.pone.0196751.g004

More quantitatively, the maximum likelihood estimates for the lifetime τml of F− and water at the at F82 and F85 sites, as computed from all replicas, are summarized in Table 1. Overall, the lifetime of water is at least one order of magnitude longer than the lifetime of F−. This finding is most pronounced for the crystallographic conformation with monobodies, for which not a single escape event for water occurred within a total of 1,600 ns, whereas F− escaped from the F82 and F85 sites with lifetimes of 181 ns and 25 ns, respectively (Table 1, left columns). In the absence of the monobody, Fluc was more flexible, thereby facilitating the escape of F− and water, as evident from the reduced lifetimes (Table 1, right columns). Here, F− escaped with lifetimes of 108 ns or 3 ns, respectively, whereas water escaped with longer lifetimes of 1,600 ns and 1,595 ns, respectively. Notably, using the CHARMM36 or the CHARMM Polarizable DRUDE force field instead of the Amber force field had only a moderate effect on lifetime of F− in absence of the monobodies (Table 1, numbers in parentheses), rendering a major uncertainty in these findings due to the fluoride-protein interactions unlikely. In summary, during free simulations, fluoride ions were unstable at the F82 and F85 sites whereas water molecules were stable, irrespective of the presence or absence of the monobodies and irrespective of the applied force field. These findings suggest that the isolated electron densities at the F82 and F85 sites represent water and not fluoride.

thumbnail
Download:
Table 1. Maximum-likelihood estimates τml (nanoseconds) for the lifetime of fluoride and water molecules at Fluc F82 and F85 positions.

https://doi.org/10.1371/journal.pone.0196751.t001

Difference in binding affinity ΔΔGbind between fluoride and water to Fluc favors water over fluoride

Based on the free simulations, water molecules are more stable inside Fluc as compared with fluoride ions. To rationalize these findings in terms of free energies, we used thermodynamic integration (TI) to compute the difference in binding free energy ΔΔGbind of water versus F− to the F82 and F85 sites of Fluc.

Fig 5A shows ΔΔGbind as blue and green bars, where a negative value indicates a more favorable binding of water. If the L2 monobodies are blocking Fluc and thereby stabilizing the crystallographic conformation, the free energy differences at F82 and F85 sites are ΔΔGbind ≈ −235 kJ mol−1, indicating a much larger probability for water being bound at those sites as compared to F− (Fig 5A, solid blue/green bars). These numbers may be rationalized by the large dehydration free energy for F− of ∼500 kJ mol−1, which are only partly compensated by F−-protein interactions [83]. In the simulations carried out in the presence of pore-blocking monobodies, water cannot enter the pore, leaving the molecules at F82 and F85 sites dehydrated (Fig 6A, 6B, 6E and 6F). In contrast, in the absence of monobodies, we observed that F− ions restrained at the F82 and F85 sites rapidly dragged water into the channel, thereby partly hydrating the F− ions; subsequently, in simulations with unrestrained F−, the F− rapidly left the channel (see Figs 4C, 6C and 6D, and Table 1). Consequently, computed ΔΔGbind values were less negative, in particular for the more solvent-exposed F85 site (Fig 5A, dashed blue/green bars). Nevertheless, ΔΔGbind still reveals a higher probability for water molecules occupying the F82 and F85 sites over fluoride anions (Fig 6G and 6H). Hence, the ΔΔGbind are compatible with the results from free simulations, favoring water over fluoride at the binding sites, irrespective whether antibodies are present or absent.

thumbnail
Download:
Fig 5. Summary of free-energy calculations.

(A) Difference ΔΔGbind of the binding free energy of water versus an ion, defined such that a negative ΔΔGbind implies a more favorable binding free energy for water. Color code as follows. Blue: fluoride/water bound to F82 site; green: fluoride/water bound to F85 site; orange: sodium/water bound to TM3 site. Solid and dashed bars show results computed in the presence and absence of the L2 blocking monobodies. (B) Absolute binding free energies ΔGbind of water to F82 (green) and F85 (blue) sites in the WT channel, and to the F82I (magenta) and F85I (dark turquoise) sites in the mutant channel. The same filled/dashed bar scheme is followed. The two values for each site correspond to the two Fluc monomers.

https://doi.org/10.1371/journal.pone.0196751.g005

thumbnail
Download:
Fig 6. Water and fluoride at F82 and F85 sites, taken from snapshots of the relative binding free energy simulations.

Several representative snapshots taken from the alchemical transformation simulations are depicted. Fluoride (A,B,C,D) is shown as purple spheres and water (E,F,G,H) as red (oxygen) and white (hydrogen) spheres. Some of the amino-acid side-chains in the F82 (A,C,E,G) and F85 (B,D,F,H) sites are depicted in colored sticks. Some of the hydrogen bonds established between the protein and the water/fluoride are highlighted with black dashed triangles.

https://doi.org/10.1371/journal.pone.0196751.g006

Absolute binding free energy ΔGbind of water to Fluc

The results of ΔΔGbind between fluoride and water indicate a much higher probability for water for being located at F82 and F85 sites as compared to fluoride. However, the ΔΔGbind values mainly reflect the large dehydration cost for fluoride, but they do not reveal whether water per se exhibits a favorable binding affinity to the F82 and F85 sites, as would be required to explain the electron densities at the F82 and F85 in terms of water. Hence, we further calculated the absolute binding free energy ΔGbind of water to F82 and F85 sites, i.e. the free energy for translocating a water molecule from bulk solvent into the F82 and F85 sites, again either in presence or in absence of the monobodies.

Fig 5B (blue and light green bars) shows ΔGbind of water molecules at the F82 and F85 sites in the WT Fluc channel. The results in Fig 5B are also summarized in Table 2. Evidently, both in systems with monobodies (Fig 5B, solid bars) and without monobodies (Fig 5B, shaded bars), ΔGbind is negative, indicating a high probability for water binding to these sites. Visual inspection shows that, at the F85 sites, water may interact with the hydroxyl groups of Ser180 and Tyr104 (as hydrogen bond donor or as acceptor), and with Asn43 (as hydrogen bond acceptor), thus rationalizing the large negative ΔGbind of approx. −20 kJ mol−1 (Fig 5B, light green bars). The F82 sites offer less interaction sites; here, water mainly interacts with the hydroxyl group of Ser112 and with the backbone carbonyl of Asn43, rationalizing the less negative ΔGbind (Fig 5B, blue bars). Taken together, the coordination of the isolated densities at the F82 and F85 sites are compatible with water binding.

thumbnail
Download:
Table 2. Absolute binding free energy ΔGbind (kJ mol−1) of water at F82, F85, F82I, and F85I sites.

https://doi.org/10.1371/journal.pone.0196751.t002

The results in Fig 5B and Table 2 reveal that ΔGbind of water to F85 site is lower than ΔGbind to F82 site regardless the absence or presence of L2 monobodies.

Rationalizing the disappearence of electron densities in the F82I/F85I mutant

Last et al. showed that, in the non-conducting F80I/F83I Fluc-Ecl mutant, the disconnected electron densities associated to F80I and F83I sites (analogues of F82 and F85 sites in Fluc-Bpe) disappear [40]. The correlation between the (i) disappearance of the electron densities and (ii) loss of function was interpreted as an indication that the electron densities represent functionally relevant fluoride binding sites. Hence, to test whether the assignment of water molecules to the binding sites, as suggested by the free energy calculations reported above, is compatible with the disappearance of the electron densities in the mutant, we further computed the absolute ΔGbind of water to F82I and F85I sites in Fluc-Bpe.

Fig 5B (magenta and dark green bars) shows that most of the ΔGbind values remain negative upon mutation, suggesting that the sites may be occupied also in the mutants. An exception is the F85I site without monobody; here, as soon as the analyzed water molecule is moved to bulk, it may rapidly be replaced by bulk water, leading to ΔGbind near zero. These findings might seem to contrast the loss of electron density in the mutant crystals. Hence, we speculated that the disappearance of the electron densities in the mutants might not reflect a loss of water binding, but instead reflect increased disorder of bound water molecules in the crystal, as the Phe→Ile mutations increase the volume of the cavities. To test this hypothesis, we ran free MD simulations for 80 ns of water at the F82 and F85 sites, and we compared disorder of the water molecules between WT and mutant. Fig 7A and 7B show positions of the water oxygen atoms taken from ∼2000 MD snapshots in WT and mutant channels. To quantify the volume sampled by water molecules at each site, we computed the root mean squared fluctuations (RMSF) of water oxygen atoms at the F82 and F85 sites in each simulation, taking the mean oxygen position as reference in each case. The average RMSF of oxygen atoms at F82 and F85 sites in WT Fluc is 0.032 nm, and 0.038 nm, respectively, characterizing a highly localized water molecule. In the mutant channel, the average RMSF of oxygen atoms at F82I and F85I is increased to 0.066 nm, and 0.063 nm, respectively. Hence, water is spatially more dispersed in the mutant as compared to WT, which may rationalize that such bound water cannot be resolved in the crystal. As such, also the disappearance of the electron densities in the F82I and F85I seems compatible with water binding to these sites.

thumbnail
Download:
Fig 7. On the disappearance of electron densities at the F85 and F82 in the mutant crystals.

Water oxygen atoms (red spheres) positions in the 2000 snapshots taken from 100-nanosecond simulations. (A) WT Fluc channel is shown in grey and white ribbons. The F82 and F85 binding sites are highlighted in both monomers. (B) Mutant Fluc channel is depicted as before. Here, the F82I and F85I mutated residues are highlighted.

https://doi.org/10.1371/journal.pone.0196751.g007

Sodium binding to the TM3 site stabilizes its structure

The isolated density at the TM3 site of Fluc is coordinated by four main chain carbonyl oxygen atoms, and it was hence modeled as a sodium ion Na+ in the 5NKQ structure [36]. As a control, we tested whether MD simulations support the assignment of Na+ to this density. Hence, we computed the relative binding free energy difference ΔΔGbind between Na+ and water at the at TM3 site (Fig 1B, yellow sphere). Fig 5A shows a ΔΔGbind = −6.2 ± 1.4 kJ mol−1, which would indicate a weak preference of water over Na+. However, because (i) ΔΔGbind is close to zero, and (ii) protein-ion interactions still exhibit increased uncertainties in modern force fields, these findings are hardly significant. As such, free energy calculations alone are in this case insufficient to assign the electron densities.

Consequently, we instead investigated the stability of the coordination shell in free simulations with either water or Na+ at the TM3 site. Fig 8A and 8B show an overlay of the snapshots of G77 and T80 residues forming the TM3 site taken from free 10-nanosecond MD simulations (thin sticks), compared with the positions of G77 and T80 in the crystal structure (thick sticks). Evidently, throughout the 10 ns simulation with Na+, the G77 and both T80 residues remain close to the crystal structure positions (Fig 8A). In contrast, during the simulation with water at the TM3 site, the G77 and T80 reveal an increased flexibility, and G77 of one chain even deviates systematically from the crystallographic position (Fig 8B, circle). Visual inspection of the trajectories shows that the carbonyl group of G77 formed a hydrogen bond with the neighboring T80 residue because water can not form favorable contacts simultaneously with all four carbonyl groups, as pointed by the authors of the 5NKQ structure [36]. In conclusion, whereas free energy calculations are on the TM3 site an inconclusive, free simulations suggest that Na+ as the more likely assignment as compared to water.

thumbnail
Download:
Fig 8. G77 and T80 residue positions in MD simulations.

(A) Positions of G77 and T80 residues in a 10 ns MD simulation in which a sodium ion is at the TM3 site. The positions of G77 and T80 residues are shown as thin sticks in which oxygen is colored in red, hydrogen in white, nitrogen in blue, and carbon in cyan. The positions of the G77 and T80 residues in the crystal structure are depicted in thick sticks following the same color scheme. The protein secondary structure is showed in white and gray ribbons. (B) Positions of G77 and T80 residues in a 10 ns MD simulation in which a water molecule is at the TM3 site. The G77 residue that adopts a significantly different orientation in the simulation compared with the crystal structure is highlighted in a dashed circle.

https://doi.org/10.1371/journal.pone.0196751.g008

The role of the force field in modeling fluoride-phenylalanine interactions

To test whether the main conclusions might be affected by force field uncertainties, we have computed relaxed PESs between fluoride and benzene using quantum-chemical calculations, and we compared the scans with results from an additive and from a polarizable force field. For the force field-based PESs, benzene was modeled from atoms of the phenyl ring of phenylalanine. As shown in Fig 9, quantum-chemical calculations reveal strong interactions of fluoride with the quadrupole of benzene, leading a potential energy of binding in vacuum of ∼ −50 kJ mol−1. The non-polarizable (additive) Amber99sb force field together with ion parameters by Joung and Cheatham strongly underestimate fluoride-phenyl interactions (Fig 9, blue curve). This suggests that the relative binding free energies between fluoride and water presented in Fig 5A overestimate the preference for water by approx. 40 kJ mol−1. However, because the additive force field yielded ΔΔGbind values in the crystallographic conformation with monobodies in the order ∼−235 kJ mol−1, such correction will not change the qualitative findings of the free energy calculations. Compared to the non-polarizable force field, the polarizable CHARMM-Drude force field models stronger fluoride-benzene interactions, in better agreement with the quantum-chemical calculations (Fig 9, orange curve). However, despite stronger fluoride-benzene interactions, the life time of fluoride at the F82 and F85 sites during CHARMM-Drude simulations was even shorter as compared to Amber-based simulations (Table 1, second numbers in parentheses), suggesting that stability assessments in this case do not critically depend on the fluoride-phenyl interactions.

thumbnail
Download:
Fig 9. Relaxed potential energy scans (PESs) of fluoride and benzene.

(A) Illustration of the coordinate R as the distance of fluoride (cyan sphere) from the center of mass of the benzene ring (green/white spheres), taken in the plane of the ring. The left cyan sphere indicates the position of the potential energy minimum. (B) Relaxed PESs using quantum-chemical calculations at different quantum levels and using different basis sets, as indicated in the legend (black, red, and green curves), revealing a potential energy minimum of ∼ −50 kJ mol−1. The non-polarizable (additive) Amber99sb force field with ion parameters by Joung and Cheatham [43] strongly underestimate fluoride-benzene interactions (blue). In contrast, the polarizable CHARMM-Drude force field [57–59] provides a PES in reasonable agreement with quantum-chemical calculations.

https://doi.org/10.1371/journal.pone.0196751.g009

Discussion

Ions, water, and other small molecules play key roles in the structure and function of many proteins, taking part in activation and de-activation processes, in enzymatic reactions, and in the formation and stabilization of structural motifs. Hence, the correct assignment of disconnected electron densities in crystal structures forms the basis for deriving functional mechanisms from structural data. However, the assignment of densities is in many cases far from obvious [9–17, 84, 85].

The present study suggests that free-energy calculations and free MD simulations may aid the assignment of disconnected electron densities in the crystal structure of proteins. More specifically, three simulation methods may serve for this purpose:

  1. Relative binding free energies ΔΔGbind between two putative molecular species, such as between an ion and water, reveal which species binds with higher affinity to the site of the disconnected electron density. Hence, ΔΔGbind values in principle clarify which interpretation of the density is more plausible. However, if the computed ΔΔGbind involves the binding of ions, ΔΔGbind is dominated by strong ion-water and, possibly, strong ion-protein interactions, leading to increased uncertainties due to imperfections in current force fields (see below). Hence, purely ΔΔGbind values that strongly deviate from zero, as found for the fluoride/water alchemical transformations found here, allow for robust conclusions.
  2. Given that ΔΔGbind favor the binding of a species A over another species B, absolute binding affinities ΔGbind (or bulk to protein transfer free energies) reveal whether the binding of the species A is plausible per se. The binding of species A with high probability, such that A is visible in the crystal, requires ΔGbind to be negative or near zero. As a side remark, the fraction of binding sites occupied by species A depends, apart from ΔGbind, also on the bulk concentration of A.
  3. The stability of structures in short free MD simulations may provide additional indications on the nature of electron densities. In particular, rapid and reproducible drifts from the crystal structure hint towards incorrect assignments of densities, or towards other problems in the structure. For instance, here we found pronounced instability of F− ions at the F82 and F85, irrespective of the applied force field (Fig 4, S1–S4 Figs), and we observed a drift in the cage of the TM3 site upon assignment of a water molecule to this site (Fig 8). In turn, however, the stability of a structure in short simulations does not necessarily imply that the assignment is correct, since structural drifts may occur in a stochastic manner and on longer time scales as compared to the simulation time. Therefore, conclusions can only be drawn from comparative stability assessments and from multiple independent simulations.

Here, we illustrated these methods by analyzing the five disconnected electron densities present in the crystal structure of the fluoride channel Fluc. In summary, the free-energy calculations and stability assessments suggest that four of the disconnected densities represent water molecules but not fluoride ions, in contrast to a provisional assignment by the authors of the crystal structure. For the fifth, central density, the results from free-energy calculations were inconclusive; however, stability assessments confirm the assignment of a sodium ion, as previously suggested [36].

The accuracy of the absolute and relative binding affinities depends on the accuracy of the applied force fields. Most critically, ion-protein interactions in additive force fields are still a source of uncertainty. For instance, the potential energy surface (PES) along a sodium-carboxylate distance in vacuum computed by older force fields may deviate from quantum-chemical calculations by tens of kilojoule per mole, which has triggered extensive force field developments in recent years [70–72]. For the present study on Fluc, we compared the interactions between fluoride and the quadrupolar field of a benzene ring computed by quantum-chemical calculations with the predictions by a non-polarizable (Amber99sb) and a polarizable (CHARMM-Drude) force field (Fig 9). The PES derived with CHARMM-Drude agreed reasonably with the PESs from quantum-chemical calculations; in contrast, Amber99sb combined with the Joung/Cheatham ion model [43] strongly underestimates fluoride–phenylalanine interactions, leaving room for future force field improvements. However, the qualitative findings of this study are not affected by these force field imperfections because (i) simulations with CHARMM-Drude confirmed the instability of F− at the F82 and F85 sites, and (ii) relative binding free energy differences between F− and water are by far larger than the uncertainties imposed by the non-polarizable force fields.

In this work, simulations have been carried out in a physiologically realistic environment including a lipid bilayer, explicit solvent molecules, and counter ions. During the alchemical transformations, this setup allowed us to maintain the overall charge neutrality of the system by simultaneously morphing water to fluoride in bulk while morphing fluoride to water in the protein. Maintaining charge neutrality avoids that a so-called background charge density is built up during the alchemical transformation, which is a consequence of the particle-mesh Ewald (PME) method for computing electrostatic interactions [86].

To carry out the ΔΔGbind calculations used here in a high-throughput fashion, it would be desirable to avoid the setup of a physiologically realistic simulation system, but to simulate the crystallographic unit cell instead. Because the unit cell might not contain sufficient volume of bulk water for placing a counter charge, such setup may require alchemical transformations that do not maintain charge neutrality, leading to artifacts due the interactions of ions with the background charge. For solvation free energy calculations of ions in a uniform solvent, analytic corrections for background charge effects are available (in addition to effects due to finite size of the simulation system and finite volume of the ion) [87–95] In addition, corrections for inhomogeneous systems of simple geometry, such as hydrophobic slaps or spheres in water, have also been proposed [86]. Corrections for more complex inhomogeneous systems, such as protein crystals, are to our knowledge not yet established, but approximations might build upon previous related work.

Supporting information

S1 Fig. Fluoride ion positions in Fluc blocked by L2 monobodies.

Green and blue curves represent the molecules starting at the F85 and F82 sites, respectively. Different shades of green and blue indicate starting positions in the two homodimers of Fluc.

https://doi.org/10.1371/journal.pone.0196751.s001

(EPS)

S2 Fig. Water molecule positions in Fluc blocked by L2 monobodies.

Green and blue curves represent the molecules starting at the F85 and F82 sites, respectively. Different shades of green and blue indicate starting positions in the two homodimers of Fluc.

https://doi.org/10.1371/journal.pone.0196751.s002

(EPS)

S3 Fig. Fluoride ion positions in Fluc free of monobodies.

Green and blue curves represent the molecules starting at the F85 and F82 sites, respectively. Different shades of green and blue indicate starting positions in the two homodimers of Fluc.

https://doi.org/10.1371/journal.pone.0196751.s003

(EPS)

S4 Fig. Water molecule positions in Fluc free of monobodies.

Green and blue curves represent the molecules starting at the F85 and F82 sites, respectively. Different shades of green and blue indicate starting positions in the two homodimers of Fluc.

https://doi.org/10.1371/journal.pone.0196751.s004

(EPS)

S5 Fig. Fluoride ion positions in Fluc free of monobodies, simulated with the CHARMM36 force field for the protein and Joung/Cheatham parameters for fluoride.

Green and blue curves represent the molecules starting at the F85 and F82 sites, respectively. Different shades of green and blue indicate starting positions in the two homodimers of Fluc.

https://doi.org/10.1371/journal.pone.0196751.s005

(EPS)

S6 Fig. Fluoride ion positions in Fluc free of monobodies, simulated with the CHARMM-Drude force field for the protein and the polarizable fluoride model by Lamoureux et al.

Green and blue curves represent the molecules starting at the F85 and F82 sites, respectively. Different shades of green and blue indicate starting positions in the two homodimers of Fluc.

https://doi.org/10.1371/journal.pone.0196751.s006

(EPS)

Acknowledgments

We thank Piotr Neumann for investigating the anomalous signal of fluoride in the Fluc structure, and we thank Randy Stockbridge for insightful discussions. This study was supported by the Deutsche Forschungsgemeinschaft (grant numbers HU 1971-1/1, HU 1971-3/1) to JSH. IA was additionally supported by the Göttingen Graduate School for Neurosciences, Biophysics, and Molecular Biosciences (GGNB). The simulations were in part performed with resources provided by the North-German Supercomputing Alliance (HLRN).

References

  1. 1. Berman HM, Westbrook J, Feng Z, Gilliland G, Bhat TN, Weissig H, et al. The protein data bank. Nucleic Acids Res. 2000;28(1):235–242. pmid:10592235
  2. 2. Rose PW, Prlić A, Altunkaya A, Bi C, Bradley AR, Christie CH, et al. The RCSB protein data bank: integrative view of protein, gene and 3D structural information. Nucleic Acids Res. 2017;45(D1):D271–D281. pmid:27794042
  3. 3. Emsley P, Lohkamp B, Scott WG, Cowtan K. Features and development of Coot. Acta Crystallogr. 2010;D66(4):486–501.
  4. 4. Chen VB, Arendall WB, Headd JJ, Keedy DA, Immormino RM, Kapral GJ, et al. MolProbity: all-atom structure validation for macromolecular crystallography. Acta Crystallogr. 2010;D66(1):12–21.
  5. 5. Laskowski RA, Rullmann JAC, MacArthur MW, Kaptein R, Thornton JM. AQUA and PROCHECK-NMR: programs for checking the quality of protein structures solved by NMR. J Biomol NMR. 1996;8(4):477–486. pmid:9008363
  6. 6. Weichenberger CX, Sippl MJ. NQ-Flipper: validation and correction of asparagine/glutamine amide rotamers in protein crystal structures. Bioinformatics. 2006;22(11):1397–1398. pmid:16595557
  7. 7. Hooft RW, Vriend G, Sander C, Abola EE. Errors in protein structures. Nature. 1996;381(6580):272–272. pmid:8692262
  8. 8. Joosten RP, Long F, Murshudov GN, Perrakis A. The PDB_REDO server for macromolecular structure model optimization. IUCrJ. 2014;1(4):213–220. pmid:25075342
  9. 9. Kleywegt GJ. Validation of protein crystal structures. Acta Crystallogr. 2000;D56(3):249–265.
  10. 10. Kleywegt GJ. Crystallographic refinement of ligand complexes. Acta Crystallogr. 2007;D63(1):94–100.
  11. 11. Cooper DR, Porebski PJ, Chruszcz M, Minor W. X-ray crystallography: Assessment and validation of protein-small molecule complexes for drug discovery. Expert Opin Drug Discov. 2011;6(8):771–782. pmid:21779303
  12. 12. Doyle DA, Cabral JM, Pfuetzner RA, Kuo A, Gulbis JM, Cohen SL, et al. The structure of the potassium channel: molecular basis of K+ conduction and selectivity. Science. 1998;280(5360):69–77. pmid:9525859
  13. 13. Zhou Y, MacKinnon R. The occupancy of ions in the K+ selectivity filter: charge balance and coupling of ion binding to a protein conformational change underlie high conduction rates. J Mol Biol. 2003;333(5):965–975. pmid:14583193
  14. 14. Kratochvil HT, Carr JK, Matulef K, Annen AW, Li H, Maj M, et al. Instantaneous ion configurations in the K+ ion channel selectivity filter revealed by 2D IR spectroscopy. Science. 2016;353(6303):1040–1044. pmid:27701114
  15. 15. Köpfer DA, Song C, Gruene T, Sheldrick GM, Zachariae U, de Groot BL. Ion permeation in K+ channels occurs by direct Coulomb knock-on. Science. 2014;346(6207):352–355. pmid:25324389
  16. 16. Nam KH, Sung MW, Hwang KY. Structural insights into the substrate recognition properties of β-glucosidase. Biochem Biophys Res Comm. 2010;391(1):1131–1135. pmid:20005197
  17. 17. Botos I, Wu Z, Lu W, Wlodawer A. Crystal structure of a cyclic form of bovine pancreatic trypsin inhibitor. FEBS Lett. 2001;509(1):90–94. pmid:11734212
  18. 18. Wlodawer A, Minor W, Dauter Z, Jaskolski M. Protein crystallography for aspiring crystallographers or how to avoid pitfalls and traps in macromolecular structure determination. FEBS J. 2013;280(22):5705–5736. pmid:24034303
  19. 19. Andrejašič M, Pražnikar J, Turk D. PURY: a database of geometric restraints of hetero compounds for refinement in complexes with macromolecular structures. Acta Crystallogr. 2008;D64(11):1093–1109.
  20. 20. Zheng H, Cooper DR, Porebski PJ, Shabalin IG, Handing KB, Minor W. CheckMyMetal: a macromolecular metal-binding validation tool. Acta Crystallogr. 2017;D73(3):223–233.
  21. 21. Weichenberger CX, Pozharski E, Rupp B. Visualizing ligand molecules in twilight electron density. Acta Crystallogr. 2013;F69(2):195–200.
  22. 22. Adams PD, Afonine PV, Bunkóczi G, Chen VB, Davis IW, Echols N, et al. PHENIX: a comprehensive Python-based system for macromolecular structure solution. Acta Crystallogr. 2010;D66(2):213–221.
  23. 23. Brunger AT, Adams PD, Clore GM, DeLano WL, Gros P, Grosse-Kunstleve RW, et al. Crystallography & NMR system: A new software suite for macromolecular structure determination. Acta Crystallogr D Biol Crystallograph. 1998;54(5):905–921.
  24. 24. Brunger AT. Version 1.2 of the Crystallography and NMR system. Nature Protoc. 2007;2(11):2728.
  25. 25. McGreevy R, Teo I, Singharoy A, Schulten K. Advances in the molecular dynamics flexible fitting method for cryo-EM modeling. Methods. 2016;100:50–60. pmid:26804562
  26. 26. McGreevy R, Singharoy A, Li Q, Zhang J, Xu D, Perozo E, et al. xMDFF: molecular dynamics flexible fitting of low-resolution X-ray structures. Acta Crystallogr D Biol Crystallograph. 2014;70(9):2344–2355.
  27. 27. Malde AK, Mark AE. Challenges in the determination of the binding modes of non-standard ligands in X-ray crystal complexes. J Comput Aided Mol Des. 2011;25(1):1–12. pmid:21053051
  28. 28. Ferguson DM, Radmer RJ, Kollman PA. Determination of the relative binding free energies of peptide inhibitors to the HIV-1 protease. J Med Chem. 1991;34(8):2654–2659. pmid:1652028
  29. 29. Tropsha A, Hermans J. Application of free energy simulations to the binding of a transition-state-analogue inhibitor to HTV protease. Protein Eng. 1992;5(1):29–33. pmid:1631042
  30. 30. Malde AK, Mark AE. Binding and enantiomeric selectivity of threonyl-tRNA synthetase. J Am Chem Soc. 2009;131(11):3848–3849. pmid:19292486
  31. 31. Stockbridge RB, Lim HH, Otten R, Williams C, Shane T, Weinberg Z, et al. Fluoride resistance and transport by riboswitch-controlled CLC antiporters. Proc Natl Acad Sci. 2012;109(38):15289–15294. pmid:22949689
  32. 32. Stockbridge RB, Robertson JL, Kolmakova-Partensky L, Miller C. A family of fluoride-specific ion channels with dual-topology architecture. Elife. 2013;2:e01084. pmid:23991286
  33. 33. Macdonald CB, Stockbridge RB. A topologically diverse family of fluoride channels. Curr Opin Struct Biol. 2017;45:142–149. pmid:28514705
  34. 34. Baker JL, Sudarsan N, Weinberg Z, Roth A, Stockbridge RB, Breaker RR. Widespread genetic switches and toxicity resistance proteins for fluoride. Science. 2012;335(6065):233–5. pmid:22194412
  35. 35. Stockbridge RB, Koide A, Miller C, Koide S. Proof of dual-topology architecture of Fluc F− channels with monobody blockers. Nat Commun. 2014;5(4):5120. pmid:25290819
  36. 36. Stockbridge RB, Kolmakova-Partensky L, Shane T, Koide A, Koide S, Miller C, et al. Crystal structures of a double-barrelled fluoride ion channel. Nature. 2015;525(7570):548. pmid:26344196
  37. 37. Turman DL, Nathanson JT, Stockbridge RB, Street TO, Miller C. Two-sided block of a dual-topology F − channel. Proc Natl Acad Sci. 2015;112(18):5697–5701. pmid:25902543
  38. 38. Turman DL, Stockbridge RB. Mechanism of single- and double-sided inhibition of dual topology fluoride channels by synthetic monobodies. J Gen Physiol. 2017; p. jgp.201611747. pmid:28258203
  39. 39. Andreini C, Cavallaro G, Lorenzini S, Rosato A. MetalPDB: a database of metal sites in biological macromolecular structures. Nucleic Acids Res. 2013;41(D1):D312–D319. pmid:23155064
  40. 40. Last NB, Kolmakova-Partensky L, Shane T, Miller C. Mechanistic signs of double-barreled structure in a fluoride ion channel. Elife. 2016;5:1–10.
  41. 41. Wolf MG, Hoefling M, Aponte-Santamaría C, Grubmüller H, Groenhof G. g_membed: Efficient insertion of a membrane protein into an equilibrated lipid bilayer with minimal perturbation. J Comput Chem. 2010;31(11):2169–2174. pmid:20336801
  42. 42. Jorgensen WL, Chandrasekhar J, Madura JD, Impey RW, Klein ML. Comparison of simple potential functions for simulating liquid water. J Chem Phys. 1983;79(2):926–935.
  43. 43. Joung IS, Cheatham TE III. Determination of alkali and halide monovalent ion parameters for use in explicitly solvated biomolecular simulations. J Phys Chem B. 2008;112(30):9020–9041. pmid:18593145
  44. 44. Schrödinger, LLC. The PyMOL Molecular Graphics System, Version 1.8; 2015.
    • 45. Lindorff-Larsen K, Piana S, Palmo K, Maragakis P, Klepeis JL, Dror RO, et al. Improved side-chain torsion potentials for the Amber ff99SB protein force field. Proteins Struct Funct Bioinform. 2010;78(8):1950–1958.
    • 46. Berger O, Edholm O, Jähnig F. Molecular dynamics simulations of a fluid bilayer of dipalmitoylphosphatidylcholine at full hydration, constant pressure, and constant temperature. Biophys J. 1997;72(5):2002. pmid:9129804
    • 47. Cordomí A, Caltabiano G, Pardo L. Membrane protein simulations using AMBER force field and Berger lipid parameters. J Chem Theory Comput. 2012;8(3):948–958. pmid:26593357
    • 48. Pronk S, Páll S, Schulz R, Larsson P, Bjelkmar P, Apostolov R, et al. GROMACS 4.5: A high-throughput and highly parallel open source molecular simulation toolkit. Bioinformatics. 2013;29(7):845–854. pmid:23407358
    • 49. Essmann U, Perera L, Berkowitz ML, Darden T, Lee H, Pedersen LG. A smooth particle mesh Ewald method. J Chem Phys. 1995;103(19):8577.
    • 50. Miyamoto S, Kollman PA. SETTLE: an analytical version of the SHAKE and RATTLE algorithm for rigid water models. J Comput Chem. 1992;13(8):952–962.
    • 51. Hess B. P-LINCS: A parallel linear constraint solver for molecular simulation. J Chem Theory Comput. 2008;4(1):116–122. pmid:26619985
    • 52. Feenstra KA, Hess B, Berendsen HJ. Improving e_ciency of large timescale molecular dynamics simulations of hydrogen-rich systems. J Comput Chem. 1999;20:786–798.
    • 53. Bussi G, Donadio D, Parrinello M. Canonical sampling through velocity rescaling. J Chem Phys. 2007;126(1):014101. pmid:17212484
    • 54. Berendsen HJ, Postma Jv, van Gunsteren WF, DiNola A, Haak J. Molecular dynamics with coupling to an external bath. J Chem Phys. 1984;81(8):3684–3690.
    • 55. Best RB, Zhu X, Shim J, Lopes PE, Mittal J, Feig M, et al. Optimization of the additive CHARMM all-atom protein force field targeting improved sampling of the backbone φ, ψ and side-chain χ1 and χ2 dihedral angles. J Chem Theory Comput. 2012;8(9):3257–3273. pmid:23341755
    • 56. Klauda JB, Venable RM, Freites JA, O’Connor JW, Tobias DJ, Mondragon-Ramirez C, et al. Update of the CHARMM all-atom additive force field for lipids: validation on six lipid types. J Phys Chem B. 2010;114(23):7830–7843. pmid:20496934
    • 57. Lamoureux G, Roux B. Modeling induced polarization with classical Drude oscillators: Theory and molecular dynamics simulation algorithm. J Chem Phys. 2003;119(6):3025–3039.
    • 58. Harder E, Anisimov VM, Vorobyov IV, Lopes PEM, Noskov SY, MacKerell AD, et al. Atomic level anisotropy in the electrostatic modeling of lone pairs for a polarizable force field based on the classical drude oscillator. J Chem Theory Comput. 2006;2(6):1587–1597. pmid:26627029
    • 59. Lemkul JA, Huang J, Roux B, Mackerell AD. An Empirical Polarizable Force Field Based on the Classical Drude Oscillator Model: Development History and Recent Applications. Chem Rev. 2016;116(9):4983–5013. pmid:26815602
    • 60. Lamoureux G, Harder E, Vorobyov IV, Roux B, MacKerell AD. A polarizable model of water for molecular dynamics simulations of biomolecules. Chem Phys Lett. 2006;418(1):245–249.
    • 61. Yu H, Whitfield TW, Harder E, Lamoureux G, Vorobyov I, Anisimov VM, et al. Simulating Monovalent and Divalent Ions in Aqueous Solution Using a Drude Polarizable Force Field. J Chem Theory Comput. 2010;6(3):774–786. pmid:20300554
    • 62. Chowdhary J, Harder E, Lopes PEM, Huang L, MacKerell AD, Roux B. A Polarizable Force Field of Dipalmitoylphosphatidylcholine Based on the Classical Drude Model for Molecular Dynamics Simulations of Lipids. J Phys Chem B. 2013;117(31):9142–9160. pmid:23841725
    • 63. Vermaas JV, Hardy DJ, Stone JE, Tajkhorshid E, Kohlmeyer A. TopoGromacs: Automated Topology Conversion from CHARMM to GROMACS within VMD. J Chem Inf Model. 2016;56(6):1112–1116. pmid:27196035
    • 64. Van Maaren PJ, Van Der Spoel D. Molecular dynamics simulations of water with novel shell-model potentials. J Phys Chem B. 2001;105(13):2618–2626.
    • 65. Lemkul JA, Roux B, Van Der Spoel D, Mackerell AD. Implementation of extended Lagrangian dynamics in GROMACS for polarizable simulations using the classical Drude oscillator model. J Comput Chem. 2015;36(19):1473–1479. pmid:25962472
    • 66. Nosé S. A unified formulation of the constant temperature molecular dynamics methods. J Chem Phys. 1984;81(1):511–519.
    • 67. Hoover WG. Canonical dynamics: equilibrium phase-space distributions. Phys Rev A: At Mol Opt Phys. 1985;31(3):1695.
    • 68. Vrbka L, Vondrášek J, Jagoda-Cwiklik B, Vácha R, Jungwirth P. Quantification and rationalization of the higher affinity of sodium over potassium to protein surfaces. Proc Natl Acad Sci. 2006;103(42):15440–15444. pmid:17032760
    • 69. Jungwirth P, Winter B. Ions at aqueous interfaces: from water surface to hydrated proteins. Annu Rev Phys Chem. 2008;59:343–366. pmid:18031215
    • 70. Luo Y, Roux B. Simulation of osmotic pressure in concentrated aqueous salt solutions. J Phys Chem Lett. 2009;1(1):183–189.
    • 71. Yoo J, Aksimentiev A. Improved parametrization of Li+, Na+, K+, and Mg2+ ions for all-atom molecular dynamics simulations of nucleic acid systems. J Phys Chem Lett. 2011;3(1):45–50.
    • 72. Rembert KB, Paterová J, Heyda J, Hilty C, Jungwirth P, Cremer PS. Molecular mechanisms of ion-specific effects on proteins. J Am Chem Soc. 2012;134(24):10039–10046. pmid:22687192
    • 73. Van Gunsteren W, Berendsen H. A leap-frog algorithm for stochastic dynamics. Mol Simulat. 1988;1(3):173–185.
    • 74. Deng Y, Roux B. Computations of Standard Binding Free Energies with Molecular Dynamics Simulations. J Phys Chem B. 2009;113(8):1–30.
    • 75. Hermans J, Wang L. Inclusion of loss of translational and rotational freedom in theoretical estimates of free energies of binding. Application to a complex of benzene and mutant T4 lysozyme. J Am Chem Soc. 1997;119(11):2707–2714.
    • 76. Hodgman CD, Frankel M, et al. Handbook of chemistry and physics. Chemical Rubber Pub. Co.; 1960.
      • 77. Hub JS, Awasthi N. Probing a continuous polar defect: A reaction coordinate for pore formation in lipid membranes. J Chem Theory Comput. 2017;. pmid:28376619
      • 78. Becke A. Becke 3 term correlation functional (part of b3lyp). J Chem Phys. 1993;98:5648.
      • 79. Weigend F, Ahlrichs R. Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy. Phys Chem Chem Phys. 2005;7(18):3297–3305. pmid:16240044
      • 80. Zheng J, Xu X, Truhlar DG. Minimally augmented Karlsruhe basis sets. Theor Chem Acc. 2011;128(3):295–305.
      • 81. Dunning TH Jr. Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. J Chem Phys. 1989;90(2):1007–1023.
      • 82. Neese F. The ORCA program system. Wiley Interdiscip Rev: Comput Mol Sci. 2012;2(1):73–78.
      • 83. Lamoureux G, Roux B. Absolute hydration free energy scale for alkali and halide ions established from simulations with a polarizable force field. J Phys Chem B. 2006;110(7):3308–3322. pmid:16494345
      • 84. Morozenko A, Stuchebrukhov A. Dowser++, a new method of hydrating protein structures. Proteins: Struc Func Bioinf. 2016;84(10):1347–1357.
      • 85. Sridhar A, Ross GA, Biggin PC. Waterdock 2.0: Water placement prediction for Holo-structures with a pymol plugin. PLoS ONE. 2017;12(2):e0172743. pmid:28235019
      • 86. Hub JS, Groot BLD, Grubmueller H, Groenhof G. Quantifying artifacts in Ewald simulations of inhomogeneous systems with a net charge. J Chem Theory Comput. 2014;10(1).
      • 87. Figueirido F, Del Buono GS, Levy RM. On finite-size effects in computer simulations using the Ewald potential. J Chem Phys. 1995;103(14):6133–6142.
      • 88. Hummer G, Pratt LR, Garcia AE. Free energy of ionic hydration. J Phys Chem. 1996;100(4):1206–1215.
      • 89. Figueirido F, Del Buono GS, Levy RM. On finite-size corrections to the free energy of ionic hydration. J Phys Chem B. 1997;101(29):5622–5623.
      • 90. Hummer G, Pratt LR, García AE. Ion sizes and finite-size corrections for ionic-solvation free energies. J Chem Phys. 1997;107(21):9275–9277.
      • 91. Hummer G, Pratt LR, García AE. Molecular theories and simulation of ions and polar molecules in water. J Phys Chem A. 1998;102(41):7885–7895.
      • 92. Sakane S, Ashbaugh HS, Wood RH. Continuum corrections to the polarization and thermodynamic properties of Ewald sum simulations for ions and ion pairs at infinite dilution. J Phys Chem B. 1998;102(29):5673–5682.
      • 93. Hünenberger PH, McCammon JA. Ewald artifacts in computer simulations of ionic solvation and ion–ion interaction: a continuum electrostatics study. J Chem Phys. 1999;110(4):1856–1872.
      • 94. Hünenberger PH, McCammon JA. Effect of artificial periodicity in simulations of biomolecules under Ewald boundary conditions: a continuum electrostatics study. Biophys Chem. 1999;78(1):69–88. pmid:10343384
      • 95. Kastenholz MA, Hünenberger PH. Influence of artificial periodicity and ionic strength in molecular dynamics simulations of charged biomolecules employing lattice-sum methods. J Phys Chem B. 2004;108(2):774–788.