Computational methods used to explore transport events in biological

Subscriber access provided by UNIVERSITY OF TOLEDO LIBRARIES

Review

Computational methods used to explore transport events in biological systems Tomasz Pie#ko, and Joanna Trylska J. Chem. Inf. Model., Just Accepted Manuscript • DOI: 10.1021/acs.jcim.8b00974 • Publication Date (Web): 17 Mar 2019 Downloaded from http://pubs.acs.org on March 17, 2019

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 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 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.

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 37 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 Information and Modeling

Computational Methods Used to Explore Transport Events in Biological Systems Tomasz Pieńko∗,†,‡ and Joanna Trylska∗,† Centre of New Technologies, University of Warsaw, S. Banacha 2c, 02-097 Warsaw, Poland, and Department of Drug Chemistry, Faculty of Pharmacy with the Laboratory Medicine Division, Medical University of Warsaw, S. Banacha 1a, 02-097 Warsaw, Poland E-mail: [email protected]; [email protected] keywords: molecular dynamics simulations; Brownian dynamics; enhanced sampling methods; alchemical free energy calculations; membrane channels, transporters and pumps; transport through membranes; substrate binding and permeation; conformational dynamics Abstract Transport of various molecules facilitated with membrane proteins is necessary for maintaining homeostasis in living cells. In humans, dysfunction of these proteins leads to many diseases. Thus, understanding how the membrane proteins function may help using them as therapeutic targets. To successfully investigate the mechanistic aspects of transport, the choice of appropriate methods is crucial. We review the computational methods that have proven most effective in investigating transport events, specifically, deterministic time-dependent classical molecular dynamics and its enhanced sampling variants, as well as methods based on Brownian dynamics. We describe technical aspects of these methods and examples of their novel variants or combinations that have ∗

To whom correspondence should be addressed Centre of New Technologies, University of Warsaw, S. Banacha 2c, 02-097 Warsaw, Poland ‡ Department of Drug Chemistry, Faculty of Pharmacy with the Laboratory Medicine Division, Medical University of Warsaw, S. Banacha 1a, 02-097 Warsaw, Poland †

1

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 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

been recently and successfully applied in the transport studies. We also discuss the difficulties related to these methods and provide possible solutions to avoid them.

2


Page 2 of 37

Page 3 of 37 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


Introduction The Nobel Prize in Chemistry awarded to Agre and MacKinnon in 2003 for the discovery of channels in the cell membrane confirmed the importance of research on transport processes in biological systems. Transport of various molecules between cellular compartments is essential for maintaining homeostasis in organisms of all domains of life. Understanding the biophysical aspects of transport is fundamental because of its biomedical significance. Especially, it could help in developing strategies for the treatment of diseases associated with abnormalities in the activity of some channels and transporters. For example, dysfunction of calcium-conducting channels may potentially lead to neurological disorders, such as epilepsy, Parkinson’s and Alzheimer’s diseases or cardiovascular disorders, and of potassium channels – to cardiac arrhythmias and atrial fibrillation. Malfunction of neurotransmitter transporters causes numerous neurological and psychiatric disorders such as schizophrenia, epilepsy, bipolar and post-traumatic stress disorder. In addition, understanding the functioning of efflux pumps, that displace drugs from cells, could help in the rational design of new agents against multi-drug resistant bacteria or cancers. Unfortunately, studying the mechanistic details of transport phenomena on the microto nano-scale in biochemical and biophysical experiments is either unfeasible or very challenging. Although the X-ray, NMR or cryo-EM studies offer valuable information on the atomic structures of membrane channels and transporters, they often lack to provide deep insight into their dynamics. Therefore, a variety of computational methods based mainly on molecular dynamics (MD) and Brownian dynamics (BD) have been applied to understand the coupling of the structure, dynamics and function of proteins involved in transport. Most recent reviews on simulating transport concentrate on a particular group of membrane proteins 1–3 or discuss system preparation, force fields and tools. 4 In our review, we primarily focus on methods that are currently most commonly used to investigate the permeation events in membrane channels or large-scale conformational transitions in transporter proteins. To the best of our knowledge, so far this topic has been most extensively reviewed 3



by Vermaas et al. 5 with a focus on MD-based methods. Comparing to this work, we expand our review with BD and describe in more detail enhanced sampling methods such as steered MD, metadynamics, adaptive biasing force, and accelerated MD. We classified the methods according to their main purpose to help the readers rapidly find solutions to the problems they are interested in (Table 1). We illustrate many technical difficulties related to these methods and discuss how they can be circumvented. We also provide examples of novel variants or combinations of the methods that have been recently and successfully applied in the transport studies. In general, this work attempts to serve as a guide especially for those who are novice in the field of simulations of transport phenomena.

System resolution and simulation engine There is no single simulation technique that could cover all spatial and temporal scales characterizing transport, so these phenomena have to be explored at different resolutions – from nuclei and electrons, through atoms to whole residues. Because electrons have dualistic corpuscular-wave nature, their evolution is governed by time-dependent Schroedinger equation. Since quantum effects decay with increasing mass of a particle, dynamics of atoms and larger entities may be characterized using classical Newtonian laws. The type of mechanics used depends on whether the process involves chemical reactions. If so, the part of the system involving the breaking or forming of chemical bonds is treated quantum mechanically (QM). Because of the enormous cost of solving QM equations, the remainder of the system is treated using classical mechanics. This computational approach is known as QM/MM 34 and was found helpful in studying the transfer of electrons, 35 protons, 27 ammonium, 20 formate, 33 and sodium 36 in membrane proteins.

4


Page 4 of 37

Page 5 of 37

Substrate binding and permeation

Table 1: List of topics in the investigation of transport processes and computational methods used to achieve particular aims with example references.

Conformational dynamics of transporters

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


Aim

Method

Examples 6–8

Find initial substrate permeation pathway

SMD MetaD BRODEA US ABF MetaD SM SMD

6,8,12,13

Calculate substrate binding free energy to channel/transporter

FEP/TI

15,19

Calculate relative binding free energy of different substrates to channel/transporter

FEP/TI

19–21

8,22–24

Find initial conformational transition of transporter

TMD SMD aMD MDeNM

Calculate PMF along conformational transition of transporter

US SM

22,25

Calculate relative free energy of mutation of transporter residue

FEP/TI

19,30,33

Calculate PMF along substrate transport pathway

5


9,10 11

7,14,15 9,10,16 10,12 17,18

25–27 15,18,28 29

30–32


Molecular dynamics (MD) When the simulated process does not involve any chemical reactions, its dynamics may be described using Newtonian mechanics in a method called molecular dynamics (MD). 37 Individual particles may represent atoms (in full atomistic models) or even groups of atoms in a coarse-grained (CG) representation. The potential which governs the motion of these particles is provided by an empirical force field, and their positions and velocities as a function of time are generated by integrating Newton’s equations of motion. To mimic biological conditions, the simulated system is coupled with a thermostat (in the NVT and NPT ensemble) and barostat (in the NPT ensemble). MD simulations are typically carried out using periodic boundary conditions to circumvent the artificial boundary effects of the finite system. The typical time step between the calculation of forces is limited to 1 fs. Extending the time step is possible by constraining the bonds involving hydrogen atoms (to 2 fs), using hydrogen mass repartitioning (to 4 fs) or virtual hydrogen sites (to 5 fs). 38 In practice, such time step limitations restrict the routine simulation scales to microseconds. Unfortunately, biological events, such as binding a substrate and its permeation through a channel or large conformational changes within transporter protein, often take place on a longer time scale. Thus, to capture these events, the system’s conformations should be extensively sampled. Nowadays, millisecond scale with classical MD simulations is only reachable if using specialized hardware such as the Anton supercomputer. 39 To further lengthen the time step, the CG representation of the system may be used. 32,40 The increased mass of single particles reduces the frequency of motions, thus the time step may reach 20–40 fs, depending on the CG resolution. Although the GC models allow extending the simulation time by at least an order of magnitude comparing to full atomistic ones, some important features, specific for the transport process, could be missing due to oversimplifying the system. Also, the CG force fields are often not readily available and their parameterization is required. 41,42 Therefore, the CG modeling is most commonly used for studying permeation of peptides through lipid 6


Page 6 of 37

Page 7 of 37 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


membranes where the specificity of intermolecular interactions is not crucial. 43,44

Brownian dynamics (BD) To further reduce the number of degrees of freedom, fast-moving solvent molecules can be replaced by implicit representation in a form of frictional and random forces that act on explicit solute particles. Friction forces are proportional to particle velocity and account for the viscosity of the medium. Random forces mimic the background noise from the thermal motions of the fluid. The relation between these forces is dictated by the fluctuation-dissipation theorem. The explicit particles experience also a force arising from their interactions. This force is calculated from the many-body potential of mean force (PMF) that is felt by a particle at a certain position, while other degrees of freedom are averaged. All these features describe dynamics of the system governed by a generalized Langevin equation. If the friction (or damping) coefficient, that determines the friction forces and scales the random forces, is sufficiently large (overdamping), the effect of particle acceleration (inertia) due to instantaneous interactions with other particles decays to zero. In this form, the overdamped Langevin equation gives the basis of Brownian dynamics (BD). 45 So far, BD simulations have been an attractive alternative to study the transport of ions across membrane channels. 46,47 BD allowed calculations of pore properties, such as ion conductance and selectivity, that are directly comparable with experiments. To include realistic boundary conditions (the concentration and transmembrane potential) in the ion channel simulations, BD was conjugated with grand canonical Monte Carlo scheme (GCMC/BD). 48 GCMC performs a series of random walks of configurations during which particles can be added or removed from the system according to the Metropolis criterion. BD was also extended to an atomistic model using high-resolution 3D maps of PMF determined from all-atom MD simulations. 49 Although GCMC/BD can routinely produce tens of microseconds of simulation time and describes ion permeation through nanopores more accurately, 50 it still has many limita7



tions. Continuum dielectric representation of solvent cannot reflect the differences between the properties of water in bulk and in extremely narrow pore confinement. Similarly, the effect of protein dielectric constant is small in wide channels, but large in narrow pores. In addition, the hydrogen bonding effects between the solute and solvent, as well as between the protein and solvent, are neglected. The assumption of a rigid protein forming a channel may be inaccurate to describe substrate permeation, especially in the presence of charged and polar amino acids that interact strongly with the substrate. To include the degrees of freedom relevant to transport, such as of residues forming the permeation pathway, the hybrid GCMC/BD-MD method was formulated. Flexibility of these residues was also implemented through modifications in the many-body PMF and Langevin equations of motion. 51 Despite the improved accuracy in describing ion permeation, this approach was unable to characterize translocation pathways of larger molecules. Recently, an approach termed BRODEA (Brownian Dynamics including Explicit Atoms), extended with temperature accelerated method, was used for the first GCMC/BD-MD simulation of antibiotic permeation through porins. 11 The BRODEA obtained transport of ciprofloxacin through the OmpC porin agreed qualitatively with all-atom MD simulations. 10 Therefore, BRODEA could be helpful as the starting point because it provides reliable data at low computational cost. However, to obtain a more detailed picture, advanced methods based on MD should follow.

Substrate binding and permeation Classical microsecond MD was shown sufficient to estimate PMF, diffusion and osmotic permeability of water through aquaporins (AQP) 52 and other channels 53 since there are no high energy barriers on the permeation path. MD simulations also successfully identified putative binding sites and delivery pathways of gaseous molecules, such as molecular oxygen, using explicit (flooding) and implicit ligand sampling. 54,55 However, classical MD is typically

8


Page 8 of 37

Page 9 of 37 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


insufficient to observe binding and transport of larger molecules, such as enzymatic co-factors, drugs, nucleic acids and proteins because they occur on time scales that range from hundreds of microseconds to even seconds. Interestingly, a recent study on the semi-SWEET protein of E. coli, which passively shuttles glucose to the bacterial cell, demonstrated that several microsecond MD simulations captured the structural transitions of the transporter and led to the passage of glucose. 56 Nevertheless, usually, to examine the binding and translocation of larger substrates or conformational transitions of the transporters on the time scale of full atomistic MD, external forces (also known as biasing forces) are introduced. The methods based on this idea are termed enhanced sampling or biased.

Steered Molecular Dynamics (SMD) One of the most common biased methods, used to identify pathways of binding, translocation, and release of a substrate from channels and pumps, is steered MD (SMD). 57 In the most frequently used SMD variant, the constant-velocity SMD (cv-SMD) (Figure 1), the force is tuned to pull the solute along the reaction coordinate with constant speed. A common problem associated with solute pulling is drifting of the whole transporter. This can be eliminated by imposing harmonic restraints on protein residues located far from the binding and permeation pathway. The cv-SMD technique has been successfully applied to find translocation pathways of a variety of substrates including sodium ions, 7,36 small molecules such as urea, 6 formate, 33 aspartate, 24 GABA, 58 glucose, 18,40 sucrose, 59 maltose derivatives, 60 cyclodextrins, 61 and drugs such as tiagabine 58 and doxorubicin. 8 For all these substrates, the coordinate defining the permeation pathway was a normal to the membrane (usually z-axis), therefore the force steering the solute was parallel to this direction. However, when the solute binding or permeation pathway deviates from linearity, the usage of cv-SMD becomes problematic. To make SMD more flexible, a collective-variable SMD (CVSMD) was introduced. 62 In CVSMD, the steering forces are calculated based on collective variables (CVs) that may represent a distance, 9



angle, dihedral, orientation, number of hydrogen bonds, radius of gyration, principal component, etc. For example, Park et al 18 used this method to dock glucose to the binding site of human GLUT1 transporter via steering the distances of several hydrogen bonds formed between the ligand and binding pocket residues. CVSMD was also applied to control large conformational changes within the transporters, such as rigid-body rotations, 25,26 and association/dissociation of protein domains. 25,27

Figure 1: The steering force Fz imposed on a molecule along the z-axis during permeation through the membrane channel with constant velocity of 1 Å/ns. In standard SMD, the pulling force is usually imposed on the center of mass of heavy atoms of the solute, which works well for pulling ions or small and rigid molecules. However, this approach may drastically distort flexible and long linear polymers, such as peptides or nucleic acids, during permeation through the pores. To tackle this problem, Grid-SMD (G-SMD) 63 was designed, in which the 3D electrostatic potential is computed from an MD simulation and applied only to the transported solute in order to accelerate its permeation. The other approach to translocate longer solutes is to use hybrid SMD (hSMD), 64 in which 10


Page 10 of 37

Page 11 of 37 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


the molecule is pulled via two mutually exclusive centers. Since SMD enhances sampling along the chosen reaction coordinate, it is possible to use SMD to find how the free energy (also known as the potential of mean force – PMF) changes upon pulling. 18 The PMF can be calculated employing Jarzynski’s equality that couples the work performed through the nonequilibrium process (SMD) and free energy (PMF), which is an equilibrium property. 65 In practice, many SMD simulations are carried out using large force constant for the steering potential to reduce fluctuations of the reaction coordinate among the trajectories (also known as a stiff spring approximation). It was demonstrated that the slower the pulling, the better the estimation of PMF. 66 Since using Jarzynski’s equality to obtain reliable PMF demands performing many simulations, a computational technique based on parallel stepwise pulling protocol and employing Jarzynski’s equality was proposed to compute free energy at a lower computational cost. 17,67,68 Considering all SMD variants, this method can enhance translocation of any substrate. Currently, SMD serves mainly as a good starting point to obtain potential pathways of substrate permeation. However, when used to assess the free energy change associated with substrate binding or translocation, SMD does not offer sufficient accuracy at a reasonable computational cost when compared with other biased methods (described further).

Umbrella Sampling (US) An enhanced sampling method that is more reliable for free energy calculations than SMD, is umbrella sampling (US). 69 In US, the reaction coordinate is partitioned into a series of smaller regions, called windows or images. Each window represents a replica of the system, whose position is restrained to a certain coordinate with a single harmonic potential (Figure 2A). Thus, a restricted range of coordinates, but along the entire reaction coordinate, is sampled. The most commonly considered coordinate used to evaluate the free energy of permeation through membrane channels or transporters is the distance between a reference point (e.g., 11



Figure 2: Distribution of harmonic potentials along the reaction coordinate (A). In classical US, windows and biasing potentials are constantly coupled (B) whereas in REUS they can occasionally exchange between neighboring pairs (C). the center of mass of a protein) and the center of mass of the permeating solute projected onto a normal to the membrane (usually z-axis). Exceptionally, additional coordinates such as orientation of the permeating molecule are considered. 8 Typically, passage of only one solute molecule at a time is simulated. However, for urea 6 and potassium ions, 12 a concomitant transport of multiple species was investigated. In such cases, the problem of projection of high-dimensional PMF to a lower dimension (e.g. 4D to 2D) may arise, so in practice the maximal number of reaction coordinates used in US is limited to three. 12 To recover an unbiased PMF from the biased simulations, the weighted histogram analysis method (WHAM) 70,71 is typically used. Sufficient overlap of the normalized distributions of 12


Page 12 of 37

Page 13 of 37 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


the coordinate values collected in the neighboring windows is required to obtain a converged free energy profile. Usually, numerous windows are distributed uniformly along the reaction coordinate, often generated by SMD, with a force constant of several to tens of kcal/(mol·Å2 ). For the distance in z, the windows are commonly created at 0.1-1.0 Å intervals. Using such an approach, one can avoid sampling discontinuities along the reaction coordinate but insufficient overlap between adjacent windows may occur (if the force constant is too large). A narrow range of the reaction coordinate sampled in the window allows decreasing the simulation time necessary for convergence. However, the more images, the more MD steps for window equilibration. Note, that in some cases, shorter simulation time of the window may lead to insufficient equilibration of the slowly evolving degrees of freedom orthogonal to the reaction coordinate, what influences the PMF. To enhance the quality of sampling, replica-exchange umbrella sampling 72 (REUS), also known as Hamiltonian replica-exchange molecular dynamics (H-REMD), can be applied. 25,27,30,60 While in classical US windows are simulated independently (Figure 2B), in REUS they are computed simultaneously within one run, and all their neighboring pairs may exchange at a regular time interval (Figure 2C). For efficient sampling, the acceptance ratio, i.e., the occurrence of successful exchanges, should span the range of 0.2–0.5. While the total computational cost of REUS is the same as of equivalent US, the use of REUS is often limited due to the need of large computing resources at a time. To reduce the computational effort of US simulations, a new approach called self-learning adaptive umbrella was designed. 13,73 It allows to iteratively create simulation windows by exploring local free energy landscape. This approach is especially advantageous for high dimensionality cases, because it limits sampling of the subspace of the reaction coordinates to regions where the PMF is lower than a certain threshold. US is often considered the gold standard in simulating transport phenomena. It was used to obtain PMF associated with transport of numerous chemical species, therefore we refer only to selected applications. US helped explore the mechanisms of transport of am-

13



monia through the aquaporin-4 channel 74 and urea through the urea transporter dvUT. 6 Multi-dimensional US revealed that simultaneous binding of multiple potassium ions is critical for the selectivity of the voltage-gated potassium channel KscA. 12 In another study on KscA, 13 the free energy calculations suggested that the key for ion conduction is the conformational change in the selectivity filter and not in the intracellular gate, as previously thought. Furthermore, US was used to investigate the impact of external electric field on the permeation of calcium ions through the voltage-gated calcium channel CaV Ab. 75 The mechanism of inhibition of the voltage-gated potassium channel Kv1.5 by anti-arrhythmic agents 76 was also studied with US. US was used to characterize the energetics of rotational mechanism in the AcrB efflux pump 8 and to identify, so far unexplored, possible sites of binding of the P-glycoprotein substrates. 77 Additionally, US was applied to determine the free energy change associated with large conformational changes of transporters, such as glycerol-3-phosphate:phosphate antiporter GlpT 30 and bacterial ABC exporter MsbA. 25

Adaptive Biasing Force (ABF) An alternative enhanced sampling strategy is provided by a method called Adaptive Biasing Force (ABF). 78 Similarly to US, ABF relies on the predefined reaction coordinate. While in US the biasing potential is predetermined, in ABF the bias evolves to flatten the PMF "on-the-fly", in the sampled region of the reaction coordinate. Therefore, in contrast to US, ABF does not require any prior knowledge of the free energy landscape. As the simulation progresses, the adaptive biasing potential is calculated from the accumulated (in bins) average force experienced by the system at any position along the reaction coordinate. Applying the bias cancels the existing free energy barriers along the reaction coordinate, facilitating its exploration. The biasing potential should be updated until the reaction coordinate is uniformly sampled. Because accomplishing a sufficiently homogeneous distribution for the entire reaction coordinate would be too time-consuming, the reaction coordinate is stratified into windows. For the distance in z, the windows are usually dis14


Page 14 of 37

Page 15 of 37 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


tributed at 2.0–6.0 Å intervals. It was demonstrated that increasing the number of windows effectively reduces the time required for PMF to converge. 79 The choice of the number of samples and width of the bin, in which forces are collected, are essential for obtaining a reliable PMF. As the number of samples increases, estimation of the average force is more accurate, but the sampling progress along the reaction coordinate is slowed down. In turn, a small number of samples provides a too coarse estimation of the average force, which may lead to a nonequilibrium state. To circumvent both of these difficulties, the reasonable number of samples to use is typically between 200 and 500. If the bin size is too large, PMF calculations are less accurate. However, small bin size requires more simulation time to yield satisfactory statistics of average force in every bin. Contrary to US, the ABF method does not require overlapping of adjacent windows. The PMF along the entire reaction coordinate is reconstructed by joining the average force gradients from individual images at their boundaries. Discontinuities of the average force gradient between adjacent windows indicate the problem with ergodic sampling. ABF is not as well-established method for calculating the free energy change during transport as US, so it has been applied less frequently than US. Recently, ABF was used to assess the energetics of dopamine release from the human DA transporter (hDAT). 15 ABF was also applied to calculate the free energy profile of such processes as sodium ion permeation through the voltage-gated sodium channel, 7 diffusion of ammonia through the bacterial ammonium transporter AmtB, 20 and permeation of formate through the FocA transporter. 33 Furthermore, ABF helped explain the mechanisms underlying selectivity of the cation-selective 5-HT3 A receptor channel 80 and voltage-dependent anion channel. 14

Metadynamics (MetaD) Similarly to ABF, metadynamics (MetaD) relies on introducing an adaptive history-dependent biasing potential imposed on a selected subset of CVs, 81 thus, no prior knowledge on the free energy surface is required. The potential is constructed from a sum of Gaussian functions 15



Page 16 of 37

periodically deposited in space of d selected CVs (Si (R)) (Eq. 1) to encourage the system to escape from the already sampled regions and visit new ones. The widths (σ) and heights (W ) of Gaussians and the frequency of updating the biasing potential (τ ) have to be wisely chosen (Eq. 1,2), which often requires experience. In the long-time limit, the biasing potential converges to the negative of the underlying free energy (Eq. 3). Since the history-dependent potential alters the probability distribution, after achieving reasonable convergence, the free energy surface as a function of CVs can be reconstructed using reweighting techniques. 82 t

Z

dt0 ω exp(−

V (S, t) = 0

d X (Si (R) − Si (R(t0 )))2 i=1

ω=

2σi2

W τ

lim V (S, t) = −F (S) + C

t→∞

)

(1)

(2)

(3)

In standard MetaD simulation, the height of added Gaussian functions is constant so the biasing potential may overfill the underlying free energy surface, and, eventually, drive the system towards high-energy regions of the reaction coordinate space (Figure 3B). To avoid this problem, well-tempered metadynamics (WTMetaD) was introduced 83 (Figure 3C), in which the biasing potential is tuned to decrease the bias deposition rate (ω) over simulation time. This method is often coupled with multiple-walker strategy 84 to enhance exploration of the CV space and speed the convergence. In this technique, several MetaD simulations (walkers) are carried out in parallel, whereas each walker contributes to the overall historydependent biasing potential. A common practice is to start with WTMetaD to achieve full permeation of the substrate in order to use different substrate locations in the CV space as an input for multiple walker WTMetaD. 9,10 So far, MetaD has been especially successful in accelerating the permeation through bacterial porins of various antibiotics, such as ampicillin and moxifloxacin, 85 meropenem, 9,16 16


Page 17 of 37 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


imipenem, 9 norfloxacin, 86 and ciprofloxacin. 10 In these studies, two CVs: the distance from the protein in z axis and substrate orientation, were used to describe translocation of the substrate. While in most cases, the 2D multiple-walker WTMetaD simulations have converged, Prajapati et al 10 reported convergence issues for high energy areas, corresponding to the constriction of the pore, even after several microseconds. To circumvent this problem, they performed 1D multiple-walker WTMetaD (biasing only the distance in z) and further recomputed the obtained free energy surface as a 2D function (including the orientation). They used a reweighting procedure of Tiwary and Parrinello. 82 D’Agostino et al. 16 showed that this procedure can be also applied to build a kinetic model of antibiotic diffusion through porins and the calculated transition rates may be directly comparable with kinetic data extracted from electrophysiology experiments.

Alchemical free energy calculations The biased methods we have described so far have one fundamental common feature, namely they all use the concept of geometrical reaction coordinate for free energy calculations and are often termed geometrical transformations. They provide information on how the free energy varies with respect to the system’s conformational changes, thus they have become so popular in the field of simulating transport events. However, these methods cannot assess the free energy change upon introducing structural modifications to the system. This problem can be solved using the, so called, alchemical transformations performed with Free Energy Perturbation (FEP) 87 and Thermodynamic Integration (TI) 88 methods. The difference between FEP and TI lies in the formalism of free energy calculation but, in practice, they are almost equivalent. Both methods rely on the changes of a parameter (typically termed λ or coupling parameter) in the Hamiltonian that modifies nonbonded interactions. There are two types of alchemical transformations, namely creation/annihilation (selected atoms are gradually formed or vanished) and mutation (one group of atoms is progressively replaced by a different group). Transformation is controlled by λ, which serves as 17



a ”reaction coordinate” connecting two distinct thermodynamic states of the system. The initial state corresponds to λ=0 and the targeted state to λ=1, with a path between these states discretized into a series of windows with λ ranging from 0 to 1. Since alchemical transformations are non-physical, difficulties may appear. For example, as the van der Waals radii of atoms are reduced during annihilation or mutation, these atoms can approach too closely to the rest of the system, so that the electrostatic interactions artificially increase and cause numerical instabilities. To circumvent this problem, decoupling of van der Waals interactions should be delayed from decoupling of electrostatic interactions. Another issue is known as ”end-point catastrophes”. It may arise if λ is close to 0 or 1 and growing atoms strongly overlap with already existing atoms, leading to a virtually infinite potential. To solve this problem, a soft-core potential is applied that gradually scales the short-range nonbonded interactions of the created atoms. Alternatively, the number of windows at the end points may be increased. After performing transformations, the free energy difference can be calculated using the Bennett acceptance ratio method. The convergence of free energy calculation may be assessed by monitoring the free energy in each window and the overlap between the densities of states in adjacent windows. To estimate the statistical error of the free energy calculation, both forward (λ 0→1) and backward (λ 1→0) transformations are considered. If the free energy change predicted from both transformations is not consistent, this suggests either inadequate restraining potentials 15 or insufficient sampling of the binding site residues that directly interact with the transformed residue. To alleviate this problem, FEP has been recently coupled with H-REMD. 19,89 In the transport simulation field, methods based on alchemical transformations are valuable for estimating binding free energies of substrates to their respective transporters, which often occupy different conformational states. For example, FEP was used to calculate the binding affinity of dopamine for human dopamine transporter in several conformations corresponding to different stages of the transport cycle. 15 FEP and TI have also been widely

18


Page 18 of 37

Page 19 of 37 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


used to assess relative binding free energies to explore the mechanisms of transporters’ selectivity, such as ammonium 20,21 and formate 33 transporters, glycerol-3-phosphate:phosphate antiporter GlpT 30 and Na+ /K+ -pump. 19

Conformational dynamics of transporters X-ray crystallography delivers high-resolution images of stable conformations of transporter proteins, however, it often cannot provide the intermediate states that occur only temporarily during the transport cycle. Thus, crucial information about the processes underlying the transport mechanism is missing. Since large-scale conformational transitions occur on time scales usually inaccessible to classical MD, some specialized biasing protocols have been developed. If intermediate conformational states are unknown and identifying CVs related to transitions is problematic, applying accelerated molecular dynamics (aMD) 90 may be advantageous. An alternative method could be MD with excited normal modes (MDeNM), 91 in which slow collective motions of the protein are kinetically enhanced so extended conformational space is sampled. 29 However, if the intermediate conformational states of transporters are experimentally resolved, and the question is to find transitions between these states, targeted MD (TMD), 92 collective variable-based SMD 62 (described earlier), and string method (SM) 93,94 could be employed.

Accelerated Molecular Dynamics (aMD) In accelerated MD, a boost potential (∆V (r)) is added to system’s potential energy (V (r)), if it is below a certain energy threshold E (Eqs. 4, 5). The energy barriers separating local minima on such modified potential energy (V ∗ (r)) surface are reduced, which accelerates sampling (Figure 3D). As a result, events occurring on the millisecond time scale can be already detected in a nanosecond aMD simulation. The boost potential is routinely applied

19



Page 20 of 37

to a dihedral energy because conformational dynamics of biomolecules arises primarily from fluctuations in torsions. However, the dual-boost mode, in which total potential energy is also affected, becomes increasingly more common.

V ∗ (r) = V (r) + ∆V (r)

∆V (r) =

   0

if V (r) ≥ E

(E − V (r))2    α + E − V (r)

if V (r) < E

(4)

(5)

The rate of acceleration is controlled by the energy threshold E and α parameter. Defining these parameters requires testing. In principle, E should be greater than the average potential energy V (r) of the system. Then, the larger the E and the lower the α, the more aggressive acceleration of the system’s dynamics. The biggest advantage of aMD is that it does not require any prior knowledge on the free energy surface and predefined reaction coordinates. Indeed, for several transporters, 15,18,28 it was demonstrated that nanosecond aMD captured large conformational changes without any a priori definition of targeted transitions. The aMD method can be used for free energy calculations but since the boost potential modifies the potential energy surface, the biased ensemble averages have to be properly reweighted to recover the unbiased ones. 95 To automatize the procedure of creating the boost potential and improving accuracy of reweighting, Gaussian accelerated MD (GaMD) has been proposed. 96

Targeted Molecular Dynamics (TMD) In targeted MD, intermediate conformations between two known structures are generated by applying, to a subset of atoms, a biasing harmonic potential (Eq. 6). During a TMD simulation, root-mean-square (RMS) distance (RM S ∗ (t)) between these structures evolves linearly from the initial to target value. At each time step, the instantaneous best-fit RMS

20


Page 21 of 37 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


Figure 3: The problem of crossing the energy barrier in classical MD (A). In MetaD, deposition of Gaussians over simulation time along the collective variable enables overcoming the energy barrier, but may also lead to overfilling the free energy surface (B). In well-tempered MetaD, deposition rate of Gaussians decreases over simulation time preventing overfilling the free energy surface (C). In accelerated MD, the boost potential flattening the free energy surface increases with a decrease in α (D). distance of the current coordinates from the target coordinates (RM S(t)) is calculated. The harmonic potential is scaled by a force constant (k) and number of atoms (N ).

U (t) =

1k [RM S(t) − RM S ∗ (t)]2 2N

(6)

In order not to hinder the conformational freedom of protein side chains, it is advisable to impose forces only on backbone Cα carbons. Low force constant of the biasing potential (1–2 kcal/(mol·Å2 ) per atom) and slow TMD simulation also allow the protein to relax during the transition. Further, good practice is to perform transitions in forward and backward directions to check if the pathways agree. Recently, an approach based on TMD conjugated with simulated annealing was proposed to optimize the intermediate structures more efficiently. 23 Simulated annealing is a global optimization algorithm which applies, consecutively, heating and controlled slow cooling of

21



the system in order to find its globally optimal conformation. To circumvent the difficulty in translocation of larger solutes in response to conformational changes of the transporter, TMD was also coupled with SMD that was used to initiate the substrate motion. 8 Targeted MD has been used to study transport cycle for a wide range of transporters including bacterial ABC transporter MsbA, 25 Na+ /K+ -pump, 19 bacterial efflux pump AcrB, 8 human P-glycoprotein, 97 human glucose transporter GLUT1, 18 glutamate-gated ionotropic receptor, 23 glutamate transporter, 24 leucine transporter LeuT, 28 pentameric ligand-gated ion channel, 31 and SecDF protein. 22 Since TMD is a nonequilibrium method, the results are only qualitative. Therefore, to obtain more quantitative data TMD is usually followed by more sophisticated biasing methods, such as US 22,25 or SM. 31

String method (SM) The transition pathway between two conformational states of a transporter 30–32 or a substrate permeation path 10,12 obtained initially with TMD, SMD, US or MetaD may be further refined using the SM method. 93,94 SM allows finding the closest minimum free energy path between two states and calculating the free energy along that path as a function of a series of defined CVs. A path (known as a string) is represented by a set of evenly distributed intermediate conformations called images. Each image in the string is evolved using a harmonicallyrestrained MD (near the initial position in the CV space). Next, a new image center is calculated from an estimate of local mean force in the CV space. To avoid clustering of images near the free energy minima, the images are held equally distant along the string. 1D free energy profile along the path is calculated via integration of the mean force along the string. The optimization of images is an iterative process that finishes when the string converges to the nearest minimum free energy path. In the most common SM variant, called SM with swarms-of-trajectories (SMwST), 98 multiple (N ) independent copies of (M ) images (a swarm of trajectories) are simulated. Contrary to original SM, the image positions in the CV space are only initially restrained 22


Page 22 of 37

Page 23 of 37 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


for several picoseconds (tres ). Subsequently, they are released into short (several picoseconds tunres ) unbiased MD simulations. Next, the image location for each swarm of trajectories is updated based on the average drift in the CV space during simulations. In both SM variants, to reduce the number of iterations (Niter ) required for string convergence, initial low free energy pathway is required. As a solution, Das et al 32 proposed a multi-step approach in which a CG pathway was first calculated using the ANMPathway method, 99 based on the anisotropic network model, 100 and then was used to generate an initial full-atomic pathway using SMD simulations, which was further refined with SMwST. A different procedure to approximate minimum free energy path was proposed by Moradi et al. 30 They developed an analysis technique, post-hoc string method (PHSM), which was applied to extract an initial minimum free energy path directly from prior REUS simulations. Then, the obtained path was relaxed using SMwST. Note, that SMwST is computationally demanding and requires massively parallel machines. For example, Das et al. 32 used 32 images with 32 copies, each equilibrated for 10 ps, simulated for 10 ps, and converged after 375 iterations, so their total computational effort equaled to N · M · (tres + tunres ) · Niter = 32 · 32 · 20 ps · 375 = 7.68 µs. Some SM variants, such as zero-temperature SM, are not so computationally expensive. In this method, the initial minimum free energy path in the CV space is minimized to the closest path using steepest descent algorithm instead of MD. Prajapati et al 10 combined MetaD, to yield free energy surface, and zero-temperature SM to find minimum free energy path of the transport of ciprofloxacin through OmpC porin. Another interesting SM variant, is based on Monte Carlo coupled with simulated annealing and was used to find the minimum free energy path of permeation of potassium ions through KscA from data generated with 3D US. 12

23



Summary Computational methods, especially based on MD, have become extremely valuable for studying rare transport events across biological membranes. At current performance, MD simulations cannot replace experiments, but they provide atomistic view and physical descriptors of the processes that are hardly detectable in experiments. On the other hand, the experimental methods, such as X-ray crystallography and more recently cryo-EM, that provide high-resolution atomic structures of transporter membrane proteins, are necessary as the simulation input. In turn, the efficacy of calculations depends on the progress in methods, as well as software and hardware development. We have described the MD and BD-based methods, together with their enhanced sampling variants, most often used to investigate transport events including substrate binding, permeation and conformational changes of the transporters. Recent examples of methods and their applications to these problems are summarized in Table 1. Even though some of these methods were formulated years ago, their novel variants are still being designed, such as self-learning adaptive US, BD with explicit atoms and accounting for protein flexibility, collective variable-based SMD, and SM with swarms-of-trajectories. Since software packages are becoming more scalable on CPU clusters, and GPU usage assures good efficiency to cost ratio, simulations of molecular transport across membranes will certainly trigger more methodology development.

Acknowledgements The authors acknowledge funding from the National Science Centre, Poland (grants UMO2014/12/W/ST5/00589 and UMO-2017/27/N/NZ1/00986). J.T. acknowledges support from the Polish-U.S. Fulbright Commission and the Kosciuszko Foundation.

24


Page 24 of 37

Page 25 of 37 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


References (1) Pothula, K. R.; Solano, C. J.; Kleinekathofer, U. Simulations of Outer Membrane Channels and their Permeability. Biochim. Biophys. Acta 2016, 1858, 1760–1771. (2) Furini, S.; Domene, C. Computational Studies of Transport in Ion Channels Using Metadynamics. Biochim. Biophys. Acta 2016, 1858, 1733–1740. (3) Miranda, W. E.; Ngo, V. A.; Perissinotti, L. L.; Noskov, S. Y. Computational Membrane Biophysics: From Ion Channel Interactions with Drugs to Cellular Function. Biochim. Biophys. Acta Proteins Proteom. 2017, 1865, 1643–1653. (4) Goossens, K.; De Winter, H. Molecular Dynamics Simulations of Membrane Proteins: An Overview. J. Chem. Inf. Model. 2018, 58, 2193–2202. (5) Vermaas, J. V.; Trebesch, N.; Mayne, C. G.; Thangapandian, S.; Shekhar, M.; Mahinthichaichan, P.; Baylon, J. L.; Jiang, T.; Wang, Y.; Muller, M. P.; Shinn, E.; Zhao, Z.; Wen, P. C.; Tajkhorshid, E. Microscopic Characterization of Membrane Transporter Function by In Silico Modeling and Simulation. Meth. Enzymol. 2016, 578, 373–428. (6) Padhi, S.; Priyakumar, U. D. Urea-Aromatic Stacking and Concerted Urea Transport: Conserved Mechanisms in Urea Transporters Revealed by Molecular Dynamics. J. Chem. Theory Comput. 2016, 12, 5190–5200. (7) Li, Y.; Sun, R.; Liu, H.; Gong, H. Molecular Dynamics Study of Ion Transport Through an Open Model of Voltage-gated Sodium Channel. Biochim. Biophys. Acta Biomembr. 2017, 1859, 879–887. (8) Vargiu, A. V.; Ramaswamy, V. K.; Malvacio, I.; Malloci, G.; Kleinekathofer, U.; Ruggerone, P. Water-mediated Interactions Enable Smooth Substrate Transport in a Bacterial Efflux Pump. Biochim. Biophys. Acta Gen. Subj. 2018, 1862, 836–845. 25



(9) Bajaj, H.; Scorciapino, M. A.; Moynie, L.; Page, M. G.; Naismith, J. H.; Ceccarelli, M.; Winterhalter, M. Molecular Basis of Filtering Carbapenems by Porins from ÃŐÂšLactam-resistant Clinical Strains of Escherichia coli. J. Biol. Chem. 2016, 291, 2837– 2847. (10) Prajapati, J. D.; Fernandez Solano, C. J.; Winterhalter, M.; Kleinekathofer, U. Characterization of Ciprofloxacin Permeation Pathways across the Porin OmpC Using Metadynamics and a String Method. J. Chem. Theory Comput. 2017, 13, 4553–4566. (11) Solano, C. J. F.; Prajapati, J. D.; Pothula, K. R.; Kleinekathofer, U. Brownian Dynamics Approach Including Explicit Atoms for Studying Ion Permeation and Substrate Translocation across Nanopores. J. Chem. Theory Comput. 2018, 14, 6701–6713. (12) Medovoy, D.; Perozo, E.; Roux, B. Multi-ion Free Energy Landscapes Underscore the Microscopic Mechanism of Ion Selectivity in the KcsA Channel. Biochim. Biophys. Acta 2016, 1858, 1722–1732. (13) Heer, F. T.; Posson, D. J.; Wojtas-Niziurski, W.; Nimigean, C. M.; Berneche, S. Mechanism of Activation at the Selectivity Filter of the KcsA K+ Channel. Elife 2017, 6, e25844. (14) Krammer, E. M.; Vu, G. T.; Homble, F.; Prevost, M. Dual Mechanism of Ion Permeation Through VDAC Revealed with Inorganic Phosphate Ions and Phosphate Metabolites. PLoS ONE 2015, 10, e0121746. (15) Cheng, M. H.; Kaya, C.; Bahar, I. Quantitative Assessment of the Energetics of Dopamine Translocation by Human Dopamine Transporter. J. Phys. Chem. B 2018, 122, 5336–5346. (16) D’Agostino, T.; Salis, S.; Ceccarelli, M. A Kinetic Model for Molecular Diffusion Through Pores. Biochim. Biophys. Acta 2016, 1858, 1772–1777.

26


Page 26 of 37

Page 27 of 37 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


(17) Ngo, V.; Stefanovski, D.; Haas, S.; Farley, R. A. Non-equilibrium Dynamics Contribute to Ion Selectivity in the KcsA Channel. PLoS ONE 2014, 9, e86079. (18) Park, M. S. Molecular Dynamics Simulations of the Human Glucose Transporter GLUT1. PLoS ONE 2015, 10, e0125361. (19) Rui, H.; Artigas, P.; Roux, B. The Selectivity of the Na(+)/K(+)-pump is Controlled by Binding Site Protonation and Self-correcting Occlusion. Elife 2016, 5, e16616. (20) Wang, S.; Orabi, E. A.; Baday, S.; Berneche, S.; Lamoureux, G. Ammonium Transporters Achieve Charge Transfer by Fragmenting Their Substrate. J. Am. Chem. Soc. 2012, 134, 10419–10427. (21) Baday, S.; Wang, S.; Lamoureux, G.; Berneche, S. Different Hydration Patterns in the Pores of AmtB and RhCG Could Determine their Transport Mechanisms. Biochemistry 2013, 52, 7091–7098. (22) Ficici, E.; Jeong, D.; Andricioaei, I. Electric-Field-Induced Protein Translocation via a Conformational Transition in SecDF: An MD Study. Biophys. J. 2017, 112, 2520– 2528. (23) Mesbahi-Vasey, S.; Veras, L.; Yonkunas, M.; Johnson, J. W.; Kurnikova, M. G. Allatom NMDA Receptor Transmembrane Domain Model Development and Simulations in Lipid Bilayers and Water. PLoS ONE 2017, 12, e0177686. (24) Wang, J.; Albers, T.; Grewer, C. Energy Landscape of the Substrate Translocation Equilibrium of Plasma-Membrane Glutamate Transporters. J. Phys. Chem. B 2018, 122, 28–39. (25) Moradi, M.; Tajkhorshid, E. Mechanistic Picture for Conformational Transition of a Membrane Transporter at Atomic Resolution. Proc. Natl. Acad. Sci. U.S.A. 2013, 110, 18916–18921. 27



(26) Lee, J.; Ren, Z.; Zhou, M.; Im, W. Molecular Simulation and Biochemical Studies Support an Elevator-type Transport Mechanism in EIIC. Biophys. J. 2017, 112, 2249– 2252. (27) Di Luca, A.; Gamiz-Hernandez, A. P.; Kaila, V. R. I. Symmetry-related Proton Transfer Pathways in Respiratory Complex I. Proc. Natl. Acad. Sci. U.S.A. 2017, 114, E6314–E6321. (28) Cheng, M. H.; Bahar, I. Complete Mapping of Substrate Translocation Highlights the Role of LeuT N-terminal Segment in Regulating Transport Cycle. PLoS Comput. Biol. 2014, 10, e1003879. (29) Araujo, G. C.; Silva, R. H.; Scott, L. P.; Araujo, A. S.; Souza, F. P.; de Oliveira, R. J. Structure and Functional Dynamics Characterization of the Ion Channel of the Human Respiratory Syncytial Virus (hRSV) Small Hydrophobic Protein (SH) Transmembrane Domain by Combining Molecular Dynamics with Excited Normal Modes. J. Mol. Model. 2016, 22, 286. (30) Moradi, M.; Enkavi, G.; Tajkhorshid, E. Atomic-level Characterization of Transport Cycle Thermodynamics in the Glycerol-3-phosphate:phosphate Antiporter. Nat. Commun. 2015, 6, 8393. (31) Lev, B.; Murail, S.; Poitevin, F.; Cromer, B. A.; Baaden, M.; Delarue, M.; Allen, T. W. String Method Solution of the Gating Pathways for a Pentameric Ligand-gated Ion Channel. Proc. Natl. Acad. Sci. U.S.A. 2017, 114, E4158–E4167. (32) Das, A.; Rui, H.; Nakamoto, R.; Roux, B. Conformational Transitions and AlternatingAccess Mechanism in the Sarcoplasmic Reticulum Calcium Pump. J. Mol. Biol. 2017, 429, 647–666. (33) Lv, X.; Liu, H.; Ke, M.; Gong, H. Exploring the pH-dependent Substrate Transport

28


Page 28 of 37

Page 29 of 37 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


Mechanism of FocA using Molecular Dynamics Simulation. Biophys. J. 2013, 105, 2714–2723. (34) Warshel, A.; Levitt, M. Theoretical Studies of Enzymic Reactions: Dielectric, Electrostatic and Steric Stabilization of the Carbonium Ion in the Reaction of Lysozyme. J. Mol. Biol. 1976, 103, 227–249. (35) Watanabe, H. C.; Yamashita, Y.; Ishikita, H. Electron Transfer Pathways in a Multiheme Cytochrome MtrF. Proc. Natl. Acad. Sci. U.S.A. 2017, 114, 2916–2921. (36) Suomivuori, C. M.; Gamiz-Hernandez, A. P.; Sundholm, D.; Kaila, V. R. I. Energetics and Dynamics of a Light-driven Sodium-pumping Rhodopsin. Proc. Natl. Acad. Sci. U.S.A. 2017, 114, 7043–7048. (37) Leach, A. R. Molecular Modelling: Principles and Applications, 2nd ed.; Pearson, 2001. (38) Loubet, B.; Kopec, W.; Khandelia, H. Accelerating All-Atom MD Simulations of Lipids Using a Modified Virtual-Sites Technique. J. Chem. Theory Comput. 2014, 10, 5690–5695. (39) Shaw, D. E.; Grossman, J. P.; Bank, J. A.; Batson, B.; Butts, J. A.; Chao, J. C.; Deneroff, M. M.; Dror, R. O.; Even, A.; Fenton, C. H.; Forte, A.; Gagliardo, J.; Gill, G.; Greskamp, B.; Ho, C. R.; Ierardi, D. J.; Iserovich, L.; Kuskin, J. S.; Larson, R. H.; Layman, T.; Lee, L.-S.; Lerer, A. K.; Li, C.; Killebrew, D.; Mackenzie, K. M.; Mok, S. Y.-H.; Moraes, M. A.; Mueller, R.; Nociolo, L. J.; Peticolas, J. L.; Quan, T.; Ramot, D.; Salmon, J. K.; Scarpazza, D. P.; Ben Schafer, U.; Siddique, N.; Snyder, C. W.; Spengler, J.; Tang, P. T. P.; Theobald, M.; Toma, H.; Towles, B.; Vitale, B.; Wang, S. C.; Young, C. Anton 2: Raising the Bar for Performance and Programmability in a Special-Purpose Molecular Dynamics Supercomputer. SC ’14:

29



Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis. 2014; pp 41–53. (40) Irudayanathan, F. J.; Wang, N.; Wang, X.; Nangia, S. Architecture of the Paracellular Channels Formed by Claudins of the Blood-brain Barrier Tight Junctions. Ann. N. Y. Acad. Sci. 2017, 1405, 131–146. (41) Leonarski, F.; Trovato, F.; Tozzini, V.; Leś, A.; Trylska, J. Evolutionary Algorithm in the Optimization of a Coarse-Grained Force Field. J. Chem. Theory Comput. 2013, 9, 4874–4889. (42) Leonarski, F.; Trylska, J. RedMDStream: Parameterization and Simulation Toolbox for Coarse-grained Molecular Dynamics Models. Biophys. J. 2015, 108, 1843–1847. (43) Zhao, J.; Zhao, C.; Liang, G.; Zhang, M.; Zheng, J. Engineering Antimicrobial Peptides with Improved Antimicrobial and Hemolytic Activities. J. Chem. Inf. Model. 2013, 53, 3280–3296. (44) Farrotti, A.; Bocchinfuso, G.; Palleschi, A.; Rosato, N.; Salnikov, E. S.; Voievoda, N.; Bechinger, B.; Stella, L. Molecular Dynamics Methods to Predict Peptide Locations in Membranes: LAH4 as a Stringent Test Case. Biochim. Biophys. Acta 2015, 1848, 581–592. (45) Ermak, D. L.; McCammon, J. A. Brownian Dynamics with Hydrodynamic Interactions. J. Chem. Phys. 1978, 69, 1352–1360. (46) Schirmer, T.; Phale, P. S. Brownian Dynamics Simulation of Ion Flow Through Porin Channels. J. Mol. Biol. 1999, 294, 1159–1167. (47) Kopfer, D. A.; Song, C.; Gruene, T.; Sheldrick, G. M.; Zachariae, U.; de Groot, B. L. Ion Permeation in K+ Channels Occurs by Direct Coulomb Knock-on. Science 2014, 346, 352–355. 30


Page 30 of 37

Page 31 of 37 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


(48) Im, W.; Seefeld, S.; Roux, B. A Grand Canonical Monte Carlo-Brownian Dynamics Algorithm for Simulating Ion Channels. Biophys. J. 2000, 79, 788–801. (49) Comer, J.; Aksimentiev, A. Predicting the DNA Sequence Dependence of Nanopore Ion Current Using Atomic-resolution Brownian Dynamics. J. Phys. Chem. C Nanomater. Interfaces 2012, 116, 3376–3393. (50) Lee, K. I.; Rui, H.; Pastor, R. W.; Im, W. Brownian Dynamics Simulations of Ion Transport Through the VDAC. Biophys. J. 2011, 100, 611–619. (51) Solano, C. J.; Pothula, K. R.; Prajapati, J. D.; De Biase, P. M.; Noskov, S. Y.; Kleinekathofer, U. BROMOCEA Code:

An Improved Grand Canonical Monte

Carlo/Brownian Dynamics Algorithm Including Explicit Atoms. J. Chem. Theory Comput. 2016, 12, 2401–2417. (52) Padhi, S.; Priyakumar, U. D. Microsecond Simulation of Human Aquaporin 2 Reveals Structural Determinants of Water Permeability and Selectivity. Biochim. Biophys. Acta Biomembr. 2017, 1859, 10–16. (53) Wong-Ekkabut, J.; Karttunen, M. Molecular Dynamics Simulation of Water Permeation Through the α-hemolysin Channel. J. Biol. Phys. 2016, 42, 133–146. (54) Oliveira, A. S.; Damas, J. M.; Baptista, A. M.; Soares, C. M. Exploring O2 Diffusion in A-type Cytochrome c Oxidases: Molecular Dynamics Simulations Uncover Two Alternative Channels Towards the Binuclear Site. PLoS Comput. Biol. 2014, 10, e1004010. (55) Mahinthichaichan, P.; Gennis, R. B.; Tajkhorshid, E. All the O2 Consumed by Thermus thermophilus Cytochrome ba3 Is Delivered to the Active Site through a Long, Open Hydrophobic Tunnel with Entrances within the Lipid Bilayer. Biochemistry 2016, 55, 1265–1278.

31



(56) Bera, I.; Klauda, J. B. Structural Events in a Bacterial Uniporter Leading to Translocation of Glucose to the Cytosol. J. Mol. Biol. 2018, 430, 3337–3352. (57) Izrailev, S.; Stepaniants, S.; Balsera, M.; Oono, Y.; Schulten, K. Molecular Dynamics Study of Unbinding of the Avidin-biotin Complex. Biophys. J. 1997, 72, 1568–1581. (58) Skovstrup, S.; David, L.; Taboureau, O.; Jorgensen, F. S. A Steered Molecular Dynamics Study of Binding and Translocation Processes in the GABA Transporter. PLoS ONE 2012, 7, e39360. (59) Sun, L.; Bertelshofer, F.; Greiner, G.; Bockmann, R. A. Characteristics of Sucrose Transport through the Sucrose-Specific Porin ScrY Studied by Molecular Dynamics Simulations. Front. Bioeng. Biotechnol. 2016, 4, 9. (60) Luo, Y.; Rossi, A. R.; Harris, A. L. Computational Studies of Molecular Permeation through Connexin26 Channels. Biophys. J. 2016, 110, 584–599. (61) Bhamidimarri, S. P.; Prajapati, J. D.; van den Berg, B.; Winterhalter, M.; Kleinekathofer, U. Role of Electroosmosis in the Permeation of Neutral Molecules: CymA and Cyclodextrin as an Example. Biophys. J. 2016, 110, 600–611. (62) Fiorin, G.; Klein, M. L.; HÃľnin, J. Using Collective Variables to Drive Molecular Dynamics Simulations. Molecular Physics 2013, 111, 3345–3362. (63) Wells, D. B.; Abramkina, V.; Aksimentiev, A. Exploring Transmembrane Transport through α-hemolysin with Grid-steered Molecular Dynamics. J. Chem. Phys. 2007, 127, 125101. (64) Chen, L. Y. Computing Membrane-AQP5-phosphatidylserine Binding Affinities with Hybrid Steered Molecular Dynamics Approach. Mol. Membr. Biol. 2015, 32, 19–25. (65) Jarzynski, C. Nonequilibrium Equality for Free Energy Differences. Phys. Rev. Lett. 1997, 78, 2690–2693. 32


Page 32 of 37

Page 33 of 37 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


(66) Park, S.; Schulten, K. Calculating Potentials of Mean Force from Steered Molecular Dynamics Simulations. J. Chem. Phys. 2004, 120, 5946–5961. (67) Ngo, V. A. Parallel-pulling Protocol for Free-energy Evaluation. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 2012, 85, 036702. (68) Ngo, V. A.; Haas, S. Demonstration of Jarzynski’s Equality in Open Quantum Systems using a Stepwise Pulling Protocol. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 2012, 86, 031127. (69) Torrie, G. M.; Valleau, J. P. Nonphysical Sampling Distributions in Monte Carlo FreeEnergy Estimation: Umbrella Sampling. J. Comput. Phys. 1977, 23, 187–199. (70) Kumar, S.; Rosenberg, J.; Bouzida, D.; Swendsen, R.; Kollman, P. Multidimensional Free-energy Calculations using the Weighted Histogram Analysis Method. J. Comput. Chem. 1995, 16, 1339–1350. (71) Roux, B. The Calculation of the Potential of Mean Force using Computer Simulations. Comput. Phys. Commun. 1995, 91, 275–282. (72) Sugita, Y.; Kitao, A.; Y., O. Multidimensional Replica-exchange Method for Freeenergy Calculations. J. Chem. Phys. 2000, 113, 6042. (73) Wojtas-Niziurski, W.; Meng, Y.; Roux, B.; Berneche, S. Self-Learning Adaptive Umbrella Sampling Method for the Determination of Free Energy Landscapes in Multiple Dimensions. J. Chem. Theory Comput. 2013, 9, 1885–1895. (74) Assentoft, M.; Kaptan, S.; Schneider, H. P.; Deitmer, J. W.; de Groot, B. L.; MacAulay, N. Aquaporin 4 as a NH3 Channel. J. Biol. Chem. 2016, 291, 19184– 19195. (75) Adiban, J.; Jamali, Y.; Rafii-Tabar, H. Simulation of the Effect of an External GHz

33



Electric Field on the Potential Energy Profile of Ca2+ ions in the Selectivity Filter of the CaV Ab Channel. Proteins 2018, 86, 414–422. (76) Chen, R.; Chung, S. H. Inhibition of Voltage-Gated K+ Channel Kv1.5 by Antiarrhythmic Drugs. Biochemistry 2018, 57, 2704–2710. (77) Subramanian, N.; Condic-Jurkic, K.; Mark, A. E.; O’Mara, M. L. Identification of Possible Binding Sites for Morphine and Nicardipine on the Multidrug Transporter P-Glycoprotein Using Umbrella Sampling Techniques. J. Chem. Inf. Model. 2015, 55, 1202–1217. (78) Darve, E.; Pohorille, A. Calculating Free Energies Using Average Force. J. Chem. Phys. 2001, 115, 9169–9183. (79) Comer, J.; Gumbart, J. C.; Henin, J.; Lelievre, T.; Pohorille, A.; Chipot, C. The Adaptive Biasing Force Method: Everything You Always Wanted to Know but Were Afraid to Ask. J. Phys. Chem. B 2015, 119, 1129–1151. (80) Di Maio, D.; Chandramouli, B.; Brancato, G. Pathways and Barriers for Ion Translocation through the 5-HT3A Receptor Channel. PLoS ONE 2015, 10, e0140258. (81) Laio, A.; Parrinello, M. Escaping Free-energy Minima. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 12562–12566. (82) Tiwary, P.; Parrinello, M. A Time-independent Free Energy Estimator for Metadynamics. J. Phys. Chem. B 2015, 119, 736–742. (83) Barducci, A.; Bussi, G.; Parrinello, M. Well-tempered Metadynamics: a Smoothly Converging and Tunable Free-energy Method. Phys. Rev. Lett. 2008, 100, 020603. (84) Raiteri, P.; Laio, A.; Gervasio, F. L.; Micheletti, C.; Parrinello, M. Efficient Reconstruction of Complex Free Energy Landscapes by Multiple Walkers Metadynamics. J. Phys. Chem. B 2006, 110, 3533–3539. 34


Page 34 of 37

Page 35 of 37 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


(85) Hajjar, E.; Kumar, A.; Ruggerone, P.; Ceccarelli, M. Investigating Reaction Pathways in Rare Events Simulations of Antibiotics Diffusion Through Protein Channels. J. Mol. Model. 2010, 16, 1701–1708. (86) Bajaj, H.; Acosta Gutierrez, S.; Bodrenko, I.; Malloci, G.; Scorciapino, M. A.; Winterhalter, M.; Ceccarelli, M. Bacterial Outer Membrane Porins as Electrostatic Nanosieves: Exploring Transport Rules of Small Polar Molecules. ACS Nano 2017, 11, 5465–5473. (87) Zwanzig, R. W. High-Temperature Equation of State by a Perturbation Method.I.Nonpolar Gases. J. Chem. Phys. 1954, 22, 1420–1426. (88) Kirkwood, J. G. Statistical Mechanics of Fluid Mixtures. J. Chem. Phys. 1935, 3, 300–313. (89) Jiang, W.; Roux, B. Free Energy Perturbation Hamiltonian Replica-Exchange Molecular Dynamics (FEP/H-REMD) for Absolute Ligand Binding Free Energy Calculations. J Chem Theory Comput 2010, 6, 2559–2565. (90) Hamelberg, D.; Mongan, J.; McCammon, J. A. Accelerated Molecular Dynamics: a Promising and Efficient Simulation Method for Biomolecules. J. Chem. Phys. 2004, 120, 11919–11929. (91) Costa, M. G.; Batista, P. R.; Bisch, P. M.; Perahia, D. Exploring Free Energy Landscapes of Large Conformational Changes: Molecular Dynamics with Excited Normal Modes. J Chem Theory Comput 2015, 11, 2755–2767. (92) Schlitter, J.; Engels, M.; Krtiger, P.; Jacoby, E.; Wollmer, A. Targeted Molecular Dynamics Simulation of Conformational Change Application to the TâĘŤR Transition in Insulin. Mol. Sim. 1993, 10, 291–308.

35



(93) Weinan, E.; Ren, W.; Eijnden, E. String Method for the Study of Rare Events. Phys. Rev. B 2002, 66, 052301. (94) Maragliano, L.; Fischer, A.; Vanden-Eijnden, E.; Ciccotti, G. String Method in Collective Variables: Minimum Free Energy Paths and Isocommittor Surfaces. J. Chem. Phys. 2006, 125, 24106. (95) Miao, Y.; Sinko, W.; Pierce, L.; Bucher, D.; Walker, R. C.; McCammon, J. A. Improved Reweighting of Accelerated Molecular Dynamics Simulations for Free Energy Calculation. J. Chem. Theory Comput. 2014, 10, 2677–2689. (96) Miao, Y.; Feher, V. A.; McCammon, J. A. Gaussian Accelerated Molecular Dynamics: Unconstrained Enhanced Sampling and Free Energy Calculation. J. Chem. Theory Comput. 2015, 11, 3584–3595. (97) McCormick, J. W.; Vogel, P. D.; Wise, J. G. Multiple Drug Transport Pathways through Human P-Glycoprotein. Biochemistry 2015, 54, 4374–4390. (98) Pan, A. C.; Sezer, D.; Roux, B. Finding Transition Pathways Using the String Method with Swarms of Trajectories. J. Phys. Chem. B 2008, 112, 3432–3440. (99) Das, A.; Gur, M.; Cheng, M. H.; Jo, S.; Bahar, I.; Roux, B. Exploring the Conformational Transitions of Biomolecular Systems Using a Simple Two-state Anisotropic Network Model. PLoS Comput. Biol. 2014, 10, e1003521. (100) Bahar, I.; Rader, A. J. Coarse-grained Normal Mode Analysis in Structural Biology. Curr. Opin. Struct. Biol. 2005, 15, 586–592.

36


Page 36 of 37

Page 37 of 37 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


Table of Contents

37


Computational methods used to explore transport events in biological

Recommend Documents