Brownian Dynamics Approach Including Explicit Atoms for Studying

Subscriber access provided by Kaohsiung Medical University

Biomolecular Systems

A Brownian Dynamics Approach Including Explicit Atoms For Studying Ion Permeation And Substrate Translocation Across Nanopores Carlos J.F. Solano, Jigneshkumar Dahyabhai Prajapati, Karunakar Reddy Pothula, and Ulrich Kleinekathöfer J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00917 • Publication Date (Web): 08 Nov 2018 Downloaded from http://pubs.acs.org on November 13, 2018

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 Theory and Computation

A Brownian Dynamics Approach Including Explicit Atoms For Studying Ion Permeation And Substrate Translocation Across Nanopores Carlos J. F. Solano,†,‡ Jigneshkumar D. Prajapati,†,‡ Karunakar R. Pothula,† and Ulrich Kleinekathöfer∗,† †Department of Physics and Earth Sciences, Jacobs University Bremen, Campus Ring 1, 28759 Bremen, Germany ‡Contributed equally to this work E-mail: [email protected]

Abstract A Brownian Dynamics (BD) approach including Explicit Atoms coined BRODEA is presented to model ion permeation and molecule translocation across a nanopore confinement. This approach generalizes our previous hybrid molecular dynamics-Brownian dynamics framework (J. Chem. Theory Comput. 12, 2401 (2016)) by incorporating a widespread and enhanced set of simulation schemes based on several boundary conditions and electrostatic models, as well as a temperature accelerated method for sampling free energy surfaces and determining substrate translocation pathways. As a test case, BRODEA was applied to study the ion diffusion as well as to ciprofloxacin and enrofloxacin transport through the outer membrane porin OmpC from E. coli. The

1

ACS Paragon Plus Environment

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

equivalence between the different simulation schemes was demonstrated and their computational efficiency evaluated. The BRODEA results are able to reproduce the main features of the ion currents and free energy surfaces determined by all-atom molecular dynamics simulations and validated by experiments. Furthermore, the BRODEA results are able to determine the ciprofloxacin and enrofloxacin permeation pathways showing a remarkable agreement with the results obtained from a computational protocol that combines metadynamics and a zero-temperature string method (J. Chem. Theory Comput. 13, 4553 (2017); J. Phys. Chem. B 122, 1417 (2018)). To our knowledge, this is the first time such antibiotic permeation pathways have been characterized by a technique based on Brownian dynamics.

1

Introduction

Biological channels are ubiquitous in nature and artificially designed nanopores, inspired from functions of biological channels, have potential applications, for example, in targeted (drug) delivery and DNA sequencing as well as in nanosensors. 1 Consequently, the key principles governing ion permeation and substrate translocation across nanoscale confinements have attracted considerable attention in recent years. 2–11 The computational methods for studying ion and small molecule permeation across biological and artificial pores can be subdivided mainly into three categories: 12 at the continuum level the Poisson-Nernst-Planck model; Brownian dynamics (BD), in which the environment is typically represented by a structureless dielectric medium; and all-atom molecular dynamics (MD) which is a fully microscopic description with all atoms treated explicitly. Recently, we have developed an hybrid MD-BD scheme 13 that incorporates a force field description for the dynamics of some explicit atoms into the grand canonical Monte Carlo/Brownian dynamics (GCMC/BD) algorithm. 14 Although the first results concerning the ion diffusion through a biological channel were encouraging, the simulations still required a considerable amount of CPU time when explicitly considering thermal fluctuations of some key amino acid residues. Moreover, the 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


characterization of substrate permeation pathways is out of reach. Here, we generalize our previous scheme in order to provide a framework for modeling ion permeation and substrate translocation across nanopores. Accordingly, we have implemented an enhanced sampling technique, as well as an extended set of simulation protocols by combining different boundary conditions (in addition to the GCMC algorithm), electrostatic models, constraint algorithms and numerical integration methods of the equations of motion. The result is the new approach coined BRODEA (BROwnian Dynamics including Explicit Atoms) that has been implemented in a code with the same name and that is based on the BROMOCEA code. 13 The software BRODEA allows simulations in several statistical ensembles, mimics the effect of the transmembrane potential by means of different methodologies, and incorporates a temperature accelerated method for characterizing the substrate permeation pathways. 15 Throughout this work, we have performed BRODEA simulations for a total simulation time of more than 740 µs using the outer membrane channel OmpC 16–18 from the Gramnegative bacterium E. coli as a test system. The OmpC porin has been classified as a nonspecific protein channel, since it allows the passage of ions and metabolites up to 600 Da with no clear specificity. Although the two most prominent outer membrane porins of E. coli, i.e., OmpF and OmpC, share around 60% sequence similarity, 17,19 an environment characterized by a high level of nutrients, such as in mammal intestines, favors the expression of OmpC over OmpF. 16 The OmpC porin has been well studied using electrophysiology experiments and simulation techniques at different conditions. 13,20,21 As shown in Fig. 1A, the OmpC porin is a trimer in which each monomer consists of a 16-stranded hollow βbarrel. An overall hourglass shape can be used to describe each monomer, where the narrow constriction region separates the extracellular (EC) and periplasmic (PP) vestibules. The constriction region is near the middle of the channel formed by the inward folded loop L3. In the constriction region, a strong electric field 20 is created as a consequence of the negatively charged residues on loop L3 (D105, D113 and E109) and the positively charged residues on the opposite side of the pore wall (K16, R37, R74 and R124) that are separated only by a

3



Figure 1: Structural features of the OmpC porin from E. coli. Loops L2, L3 and L4 are highlighted in blue, magenta and green, respectively. (A) Top view of the OmpC trimer in cartoon representation. (B) Side and (C) top views of an OmpC monomer in cartoon representation. The negatively charged residues (C atoms in green and O atoms in red) located on loop L3 (D105, D113 and E109) and the positively charged residues (C atoms in yellow and N atoms in blue) on the opposite barrel wall (K16, R37, R74 and R124) are shown as sticks.

short distance (see Figs. 1B and 1C). The present paper is organized as follows. Section 2 is devoted to providing the theoretical foundations underlying the algorithms within the BRODEA approach. Section 3 describes the novel features included into the BRODEA code while Section 4 provides the main results and is divided into three parts. In the first one, the equivalence between the simulation schemes is tested against the results from our previous scheme. 13 Thereafter, the BRODEA results for ion permeation are compared with those extracted from all-atom MD simulations and experimental results. In the third part, the BRODEA results for the ciprofloxacin and enrofloxacin permeation are compared with those obtained from a computational protocol that combines metadynamics and a zero-temperature string method. 22,23 Finally, conclusions and some remarks are given in Section 5.

2

Theory

As aforementioned, our main goal is to develop a computational approach that enables one to study ion permeation and substrate translocation across nanopores, i.e., artificial or

4


Page 4 of 37

Page 5 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


biological channels, embedded in a membrane and surrounded by aqueous salt solutions. To this end, the system under investigation is divided into implicit and explicit or so-called Brownian particles. The evolution for the Cartesian coordinates of N Brownian particles r N = {r1 , . . . , r3N } is treated on a coarse time scale as a Markov process, which is the basis of the theory of Brownian motion, 24 and determined by overdamped Langevin equations (see Section 2.1). Implicit particles are either described as a continuum dielectric, i.e., water molecules, or remain at fixed positions, e.g., the implicit pore and membrane atoms. As a consequence, implicit particles do not need to be propagated in time by the respective equations of motion but their effect should be properly incorporated into the many-body potential energy function as described in Section 2.1. Section 2.2 focuses on the fact that a proper non-equilibrium simulation algorithm needs to provide a satisfactory description of the thermal equilibrium state of the system in the absence of a net ion flux. For this reason, we put a particular emphasis on the available statistical ensembles in our BD framework as well as on those boundary conditions required to simulate these ensembles. Although the present hybrid MD-BD scheme significantly reduces the number of degrees of freedom of the system under investigation, the modeling of substrate diffusion across a nanopore is still a formidable task since it can take milliseconds or even more. Thus, a method for efficiently sampling the free energy landscape and determining the translocation pathways can be very helpful. A large variety of enhanced sampling techniques have been proposed in the literature. 25,26 In Section 2.3, we focus on a temperature accelerated method 15 as part of the BRODEA approach.

2.1

Equations of motion and many-body potential energy function

In this section we concentrate on the Smoluchowski equation, i.e., the Fokker-Planck equation in coordinate space, which describes the time evolution of the probability distribution of the Brownian particle positions, and the coupled Langevin equations which essentially generate the trajectories of these particles. Furthermore, we describe the many-body potential energy 5



Page 6 of 37

function that determines the properties of the system. Fokker-Planck equations for many interacting particles have been derived from microscopic considerations for systems of massive simple particles 27 and for polymeric systems. 28 From the Fokker-Planck equations one may derive the reduced differential equations in coordinate space known as Smoluchowski equations. 27–29 It is noteworthy that the existence of Fokker-Planck equations is sufficient, but apparently not necessary, to permit the derivation of Smoluchowski equations. 30,31 Let ρ(r N , t | r0N , 0) be the conditional space probability distribution that the particles adopt a configuration r N ≡ r N (t) at time t given the initial configuration r0N ≡ r N (0) at time t = 0. The time evolution of ρ(r N , t | r0N , 0) can be described by the Smoluchowski equation 27–29

∂ ˆ − Or ρ(r N , t | r0N , 0) = 0 ∂t

(1)

ρ(r N , 0 | r0N , 0) = δ(r N − r0N ) .

(2)

with initial condition

The Smoluchowski operator

ˆr = O

3N X ∂ ∂ Di − βfi ∂ri ∂ri i=1

(3)

acts on the variables r N . Note that Di is the self-diffusion constant of the Brownian particle associated with index i, fi is the force acting in direction i and β as usual the inverse of the thermal energy (β = 1/(kB T )). Eq. 1 can be deduced from a more general diffusion equation 27–29 after considering that: (i) hydrodynamic interactions among Brownian particles are neglected, and (ii) self-diffusion constants for Brownian particles depend on the position along the channel axis. The formal operator solution is given by

ˆ r t)δ(r N − r N ) . ρ(r N , t | r0N , 0) = exp(O 0

6


(4)

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


Thus, Eq. 1 is analogous to the Liouville equation in classical mechanics. 32 The Smoluchowski equation described in Eq. 1 is equivalent to the Langevin equation given by 33–35 p ∂ dri N N N = βDi (r )fi (r ) + CS Di (r ) + 2Di (r N )ηi (t) . dt ∂ri

(5)

Depending on the convention adopted for the stochastic differential equation (SDE), the constant CS either has the value of 1/2 in case of the Stratonovich interpretation or the value of 1 in case of the Îto interpretation. 36 The first term on the right hand side of Eq. 5 is denominated as coarse-grained velocity or drift term, while the second one is usually called Brownian velocity or diffusion term. It is important to realize that ηi is a Gaussian white noise process whose first and second moments are given by hηi (t)i = 0 and hηi (t0 )ηj (t00 )i = δ(t0 − t00 )δij , and that the integration of the SDE given by Eq. 5 leads to the same results independent of Stratonovich or Îto interpretation (see Section 1 in the Supporting Information (SI)). Several numerical methods for solving Eq. 5 are available in the BRODEA code (see Section 2 in the SI). In Eq. 5, the force fi is given by the negative derivative of the many-body potential energy function W with respect to ri . The potential is given by W (r N ) =

X

Kb (b − b0 )2 +

bonds

+

X

X

X

Kθ (θ − θ0 )2 +

angles

KS (S − S0 )2

U rey−Bradley

Kϕ (1 + cos (nϕ − δ)) +

X

Kω (ω − ω0 )2 +

impropers

dihedrals N X

X

UCM AP

residues

φrf (rα ) + Ucore (rα ) + qα φsf (rα ) + 2 α " 12 6 # N X σαβ σαβ qα q β − + wsr (rαβ ) + . + 4αβ r r 4πε αβ αβ 0 C(rαβ )rαβ β>α

(6)

As can be seen, the many-body potential energy function is composed of quite a number of terms. These include the internal terms, i.e., bond lengths (b), valence angles (θ), UreyBradley distances (S), dihedral angles (ϕ), improper angles (ω), and the CMAP backbone 7



torsional corrections 37,38 as used the CHARMM force field. 39 It should be noted that other force fields could be used whenever they share the same functional forms for their internal terms with those in the CHARMM force field. Further contributions due to external potentials need to be taken into account, including the electrostatic potential φsf from the implicit protein charges and the transmembrane potential, the reaction field potential φrf arising from the dielectric hetero-junction between solvent, implicit protein and lipid, and the core-repulsive steric potential Ucore from the implicit pore and membrane. Moreover, the many-body potential energy function includes non-bonded terms between Brownian particles, i.e., Lennard-Jones, water-mediated short-range, 40 and dielectric-screened Coulombic electrostatic terms. In addition, holonomic bond distance constraints (see Section 3 in the SI) and harmonic restraints (see Section 4 in the SI) can be imposed on the Brownian particles. Note that rα and qα denote the position and (fixed) charge of Brownian particle α, rαβ the distance between Brownian particles α and β, ε0 the vacuum permittivity, and C(r) a distance-dependent dielectric screening function. 41–44

2.2

Ensembles and boundary conditions

By definition, our prototypical system corresponds to the neighborhood of a channel embedded in a membrane and in contact with surrounding aqueous salt solutions. Thus, the system of interest constitutes of an open system with a fluctuating number of ions. For this reason, one requires to model the particle entry and exit at both channel ends during the simulations. In principle, a grand canonical ensemble is adequate to mimic an open system with a variable particle number. Accordingly, the GCMC algorithm was combined with the BD method to allow for fluctuating ion numbers. 14 In this algorithm, the particle creation and destruction is accepted based on the energy of the system and the chemical potential of the ion species. Furthermore, the (unphysical) creation and destruction is restricted to buffer regions or control cells defined at the boundaries of the simulation box and located away from the main simulation region. It is worth mentioning that reflecting boundary conditions or 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


elastic walls also need to be incorporated into the GCMC/BD scheme to preserve the finite size of the system (see Section 3.1). In spite of its efficiency in providing satisfactory concentrations even for low ion concentrations, the GCMC algorithm has some major drawbacks. As pointed out in previous studies, 14,45,46 the GCMC algorithm is CPU time consuming since the total system energy needs to be calculated in any attempt of inserting or annihilating ions in the control cells. Moreover, a certain percentage of these creation/annihilation attempts are rejected based on the Metropolis criterion which means CPU time is wasted in futile attempts. The use of different boundary conditions can be justified based on the screened nature of the electrostatics interactions in BD simulations, 45 which is quite different from those in MD. Namely, an ion’s electric field is efficiently screening by water and almost vanishing beyond a few Debye lengths due to the shielding of counterions. Provided the system boundaries are far away from the main simulation region, a simple implementation of the boundary conditions should then suffice. Indeed, it has been shown that the treatment of the system boundaries is mostly irrelevant to the conductance properties of a channel as long as these boundaries are at a reasonable distance, e.g., 3–4 Debye lengths from the channel. 45 A simple and practical approach to ensure continuity and to avoid the need to account for particle number fluctuations is to impose periodic boundary conditions (PBCs) in the BD simulations. This implies to replace the grand canonical ensemble by a canonical one. As compared to the GCMC algorithm, PBCs are conceptually simpler, involve less calculations, and are considerably faster at moderate and high ion concentrations. However, PBCs only allow for modeling symmetric concentrations of salt solutions at both sides of the channel. We can summarize by saying that the grand canonical ensemble requires GCMC boundaries that allow fluctuations in the number of ions, whereas the canonical ensemble requires PBCs that maintain a fixed number of Brownian particles in the system. A third type of boundaries, i.e., the particle counting (PACO) algorithm, 46 has been proposed recently to deal chiefly with non-equilibrium simulations. However, this kind of boundary conditions

9



Page 10 of 37

cannot easily be connected to a specific ensemble. Instead, they represent an intermediate case between grand canonical and canonical ensembles at equilibrium conditions.

2.3

Accelerated sampling for substrate translocation

Let θ(r N ) = (θ1 (r N ), . . . , θm (r N )) be a set of collective variables (CVs) which is adequate to describe the substrate translocation across a nanopore confinement. Let z = (z1 , . . . , zm ) and F (z) be a particular realization of these CVs and the free energy of the system, respectively. If the sampling of F (z) is hindered by free energy barriers between metastable basins higher than the thermal energy, then going over these barriers is a rare event and we cannot monitor the evolution of θ(r N ) from Eq. 5 to compute F (z). To overcome this difficulty, Maragliano and Vanden-Eijnden devised a temperature accelerated method 15 in which the CVs are treated as dynamical variables and an artificially high temperature for these CVs allows to accelerate the sampling of the free energy surface (FES). In particular, the many-body potential energy function is replaced by m

1X κj (zj − θj (r N ))2 Wκ (r , z) = W (r ) + 2 j=1 N

N

(7)

where {κj }j=1,...,m is a set of non-negative parameters, and the evolution of z(t) is described by the overdamped dynamics dzj ¯ j fj (z) + = β¯D dt

q

¯ j η z (t) 2D j

(8)

¯ j=1,...,m denotes a set of where β¯ = 1/(kB T¯) is the inverse of an artificial temperature T¯, {D} self-diffusion constants of the CVs, the forces fj (z) are given by fj (z) = −∂Wκ (r n , z)/∂zj = κj (θj (r N ) − zj ) and ηjz (t) is a white-noise independent of ηi (t) (see Eq. 5). When the self¯ j=1,...,m are much smaller than those from the standard BD scheme diffusion coefficients {D} {Di }i=1,...,N so that the CVs z(t) evolve much slowler than the r N and {κj }j=1,...,m are large

10


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


enough so that zj (t) ≈ θj (r N ) at all times, it has been demonstrated that the solution of Eq. 8 for the CVs z(t) satisfies approximately the effective equation dzj ¯ j ∂F (z) + = −β¯D dt ∂zj

q

¯ j η z (t) , 2D j

(9)

i.e., the CVs z(t) move by a dynamics consistent with the probability density function ¯ ¯ e−βF (z) . Therefore, the sampling of the FES can be greatly improved by adjusting β.

3

Implementation details

This section describes the novel features included into the BRODEA code. The BRODEA code has been developed based on the BROMOCEA code. 13 Thereby, BRODEA includes the functionalities and modules previously implemented in BROMOCEA.

3.1

Simulation schemes

In this subsection, we describe the simulation protocols implemented into the BRODEA code which are available for equilibrium conditions, i.e., in absence of a net ion flux, and for nonequilibrium conditions, i.e., under the effect of a transmembrane potential 3 (see Table 1). These simulation schemes share some common features with each other. First, the bonded interactions among Brownian particles are described by a force field, e.g., CHARMM 39 in our particular case. Second, non-bonded interactions between Brownian particles are composed of Lennard-Jones and dielectric-screened electrostatic terms. In addition, water-mediated short-range interactions can also be included to consider hydration effects to some extent. Third, the reaction field potential φrf is estimated either from a precomputed reaction field matrix 47 or from a set of reaction field parameters precomputed and stored on a mesh grid. Finally, Ucore is either precomputed and stored on a mesh grid 14 or calculated on the fly using Lennard-Jones potentials. 13 Although the GCMC/BD scheme was introduced in our previous work, 13 it is also described here for the sake of completeness and clarity. 11



Table 1: List of simulations schemes available in the BRODEA code. In case of the GCMC/BD scheme, the static field φsf is calculated using a modified Poisson-Boltzmann (PB-V) equation to account for the transmembrane potential V. In case of the PBC/BD and the PACO/BD schemes, φsf is determined using the Poisson equation and a homogeneous electric field (P + E-field). For equilibrium simulations the transmembrane potential is simply set to zero. Scheme GCMC/BD PBC/BD PACO/BD

System boundaries GCMC PBCs PACO

Static field PB-V P + E-field P + E-field

In the GCMC/BD scheme, the simulation domain is divided into a transport region and control cells located away from the two ends of the channel. The control cells are kept in equilibrium with the bulk solution with which they are in contact via the GCMC algorithm. The stochastic BD trajectories of the Brownian particles in the transport and control cell regions are generated by numerically solving the Langevin equations. To confine the positions of the Brownian particles inside the simulation box, reflecting boundary conditions are applied after every BD step. Moreover, the electrostatic potential φsf is calculated using a modified Poisson-Boltzmann equation to account for the membrane voltage, 3 which reduces to the Poisson-Boltzmann equation under equilibrium conditions. In addition, the reaction field φrf can also include the effect of the implicit salt outside the simulation box. 47 Note that one cycle of the GCMC/BD scheme can be composed of several GCMC and BD steps. Various geometries for the simulation domains, i.e., spherical, cylindrical or orthorhombic, are available in the BRODEA code for this scheme. In the PBC/BD scheme, the stochastic BD trajectories of Brownian particles are generated by solving the Langevin equations in an orthorhombic simulation box where PBCs are applied. Thus, the number of Brownian particles remains constant during the simulation. The static field potential φsf is composed of a potential arising from the implicit charge distribution of the channel, which is calculated using the Poisson equation, and of a linear potential along the channel axis whose slope determines the constant electric field. A homogeneous electric field perpendicular to the membrane plane is applied across the entire 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


system for describing non-equilibrium steady states. This method is equivalent to having surrounding infinite baths kept at a voltage difference via ion-exchange electrodes connected to an electromotive force and it has been tested extensively for MD simulations. 48,49 Although the contribution of φsf due to the implicit channel is computed on a finite mesh grid and hence is not fully consistent with the absence of walls as imposed by PBCs, this should not be a serious problem provided that the limits of the grid are located far enough from the main region. Recently, the PACO algorithm has been proposed for imposing boundary conditions in BD simulations. 46 As compared to the GCMC algorithm, the computational time saving of the PACO algorithm is due to two main reasons. First, ions are inserted into and deleted from the control cells only when needed. Second, the PACO algorithm only needs ion counting inside the control cells without the computationally expensive energy calculation. In the PACO/BD scheme, the simulation domain is divided into a transport region and control cells in an analogous fashion to the GCMC/BD scheme. In addition, the orthorhombic simulation box needs to be extended beyond the control cells so that outer regions are added at both sides of the channel and linked by PBCs. The outer regions form a region analogous to the transport region, whose purpose is to act as a buffer region between the control cells to keep them as thermodynamically independent as possible. 46 The concentration gradients in the outer regions are not a problem since these regions should be far from the transport region.

3.2

Temperature accelerated BD method

For improving the sampling of the FES and determine the substrate transition pathways across nanopores, the extended system of equations composed of Eqs. 5 and 8 has been

13



Page 14 of 37

discretized in the BRODEA code as follows ri (t + lδt) = ri (t + (l − 1)δt) + βDi (r N (t + (l − 1)δt))fi (r N (t + (l − 1)δt))δt ∂Di (r N (t + (l − 1)δt)) δt ∂ri p + 2Di (r N (t + (l − 1)δt))δtZ((t + (l − 1)δt)), l = 1, . . . , L , q ¯ ¯ ¯ j ∆tZ(t) , zj (t + ∆t) = zj (t) + β Dj hfj (z(t))i∆t + 2D +

(10)

where L ≥ 1 is the total number of BD steps for each single step in the CV space, δt is the BD time step, ∆t = Lδt, and Z is a standard normal random variable. Moreover, we have ∂Wκ (r N (t + (l − 1)δt), z(t)) ∂ri N ∂W (r (t + (l − 1)δt) =− ∂ri m X ∂θj (r N (t + (l − 1)δt)) + κj zj (t) − θj (r N (t + (l − 1)δt)) ∂ri j=1

fi (r N (t + (l − 1)δt)) = −

(11)

and L

1 X ∂Wκ (r N (t + (l − 1)δt), z(t)) hfj (z(t)i = − L l=1 ∂zj L κj X = θj (r N (t + (l − 1)δt)) − zj (t) . L l=1

(12)

It should be noted that the Milstein and the predictor-corrector methods (see Section 2 in the SI) are also available for solving the overdamped Langevin equation in the configurational space and that the Euler-Maruyama method is used in the CV space for sake of simplicity. Three different types of CVs have been implemented in the BRODEA code (see Section 5 in the SI). This set can easily be extended in the future to include other types if desired. In addition, full and half-harmonic restraints can be applied on selected CVs so that they remain located inside specific regions in the CV space or are restrained at a given CV value.

14


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


3.3

Milstein method for solving the Langevin equation and Verlet list for non-bonded interactions

In addition to the Euler-Maruyama approach and the predictor-corrector method, i.e., an improved version of the second-order Runge-Kutta algorithm, 13 the BRODEA code allows to integrate the Langevin equation using the Milstein method (see Section 2 in the SI). In principle, the Milstein scheme increases the accuracy of the integration as compared to the Euler-Maruyama scheme but it is less time demanding than the predictor-corrector scheme since only a single evaluation of the forces on Brownian particles is required. It is worth stressing that the Euler-Maruyama and Milstein methods are equivalent in case of non-position dependent self-diffusion constants. As does its predecessor, the BRODEA code employs a cutoff scheme, i.e., the atom-based force-switching method, for calculating the non-bonded interactions between Brownian particles. However, the calculation to determine which particle pairs fall within the cutoff distance can, itself, be time-consuming. Verlet 50 first introduced the idea of reducing the required frequency of this calculation by extending the spherical cutoff region about each particle with an additional volume shell. Thereby, particle pairs separated by a distance below an outer cutoff radius are stored in a non-bonded list but only those separated by a distance less than an inner cutoff radius are computed for non-bonded interactions. Between consecutive BD steps, particle positions do not change drastically unless unusual and unphysical long jumps 13 are not properly avoid. In principle, the same list may be used until a pair of particles moves from beyond the outer cutoff to within the inner cutoff. At the very least, one particle distance must have decreased by the width of the buffer shell before the list needs to be recalculated. It is worth noting that the use of a non-bonded list is efficient in the GCMC/BD and PACO/BD schemes only when the system is propagated several time steps by the BD algorithm before applying the ion insertion and/or deletion in the control cells.

15



3.4

Holonomic bond distance constraints and harmonic restraints

In the BRODEA code, two algorithms are available for holonomic bond distance constraints: a version of the SHAKE algorithm, which was previously implemented in the BROMOCEA code, 13 and a version of the LINCS algorithm (see Section 3 in the SI). Furthermore, BRODEA enables one to apply harmonic restraints on Cα atoms and/or a set of explicit atoms specified by the user (see Section 4 in the SI).

4

Results

This section is split into three parts. In the first part, we compare the different simulation schemes detailed above by modeling ion permeation using the same model system, i.e., OmpC, as in our previous study. 13 In the second part, the model of the OmpC system is substantially improved and the results for ion permeation are compared to all-atom MD simulations and experimental values. In the third part, the diffusion routes of ciprofloxacin and enrofloxacin molecules across OmpC are characterized using the temperature accelerated BD method and the results are compared to those recently provided by combining metadynamics and a zero-temperature string method. 22,23

4.1

Equivalence between simulation schemes

Since we are interested in comparing different BD simulation schemes and determining their efficiency, the specific model channel for these tests is of little concern. Thus, for simplicity, we use the same system setup and simulation details as described in our previous work for studying the ion permeation through an OmpC monomer that includes explicit atoms undergoing thermal fluctuations in the constriction region (see Section 5.6 in Ref. 13). Further details on the system setup and BD simulation protocols are given in Section 6.1 in the SI. Table 2 and S1 (see Section 6.2 in the SI) list the results for the ion currents and conductances obtained by applying different simulation procedures to describe a non-equilibrium 16


Page 16 of 37

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


steady state. As can be seen, the results are identical within the statistical errors irrespective of the simulation scheme, the algorithm for applying holonomic bond distance constraints and the integration method for solving the overdamped Langevin equations. Furthermore, Table 2 provides an estimation of the amount of CPU time required in each case. As compared to the performance in our previous work, 13 a significant reduction of CPU time, i.e., between 66% and 78%, is obtained in the GCMC/BD scheme when invoking the GCMC algorithm with a smaller frequency and using a Verlet list for nonbonded interactions. In addition, the CPU time requirements are also substantially decreased, i.e., between 17% and 29%, when the GCMC/BD scheme is replaced by the PBC/BD scheme. Despite the PACO algorithm being less time demanding than the GCMC algorithm, the PACO/BD simulations require more CPU time than improved GCMC/BD simulations. As noted earlier, PACO/BD simulations involve an extra pre-equilibration step and a greater simulation box that includes outer regions. Still, one can expect that PACO/BD simulations will require a smaller amount of CPU time when increasing the number of Brownian particles and/or the production times. It is worth mentioning that an optimization of control cells and outer regions sizes for improving performance is not the primary purpose of this work. The CPU times for the SHAKE and LINCS constraint algorithms are very similar for a given simulation scheme and integration method for solving the Langevin equations. However, the LINCS algorithm exhibits the advantage of being easily parallelizable. 51 Fig. 2 shows the average 1D multi-ion potential of mean forces (PMFs) along the channel axis 13 extracted from different simulation procedures using equilibrium conditions. Again we see how similar the results are within the statistical errors irrespective of the simulation scheme and the integration method for solving the overdamped Langevin equations. Moreover, the PMFs are identical to those extracted using the same conditions in our previous work (see Fig. 7 in Ref. 13). In addition, identical PMFs within the statistical errors were also obtained using the LINCS algorithm for the holonomic bond distance constraints (data not shown).

17



Page 18 of 37

Table 2: Average ion currents I, conductance values G and total CPU time T for simulations with a 300 mM KCl salt solution at a transmembrane potential of Vmp = 100 mV and at 300 K. The numerical results for the Euler-Maruyama (EM) solver of the Langevin equations are listed. The results for the other two solvers, i.e., Milstein (MLS) and predictor-corrector (PC), can be found in Section 6.1 in the SI. In addition, the errors are shown in terms of the standard deviation. SCHEME GCMC/BD (Ref. 13) GCMC/BD GCMC/BD PBC/BD PBC/BD PACO/BD PACO/BD

4.2

CONSTRAINT SHAKE SHAKE LINCS SHAKE LINCS SHAKE LINCS

IK+ [pA] 69±20 66±26 78±17 61±26 70±19 73±25 65±24

ICl− [pA] 3±2 4±3 3±2 2±2 2±1 5±2 3±3

G [nS] 0.73±0.21 0.70±0.27 0.81±0.17 0.63±0.26 0.72±0.20 0.78±0.26 0.68±0.25

T [hours] 155.59±3.86 34.21±0.47 33.62±0.17 24.18±0.25 24.00±0.20 47.68±0.54 52.83±3.55

Ion permeation

So far we have demonstrated the equivalence between different BD simulation schemes. Under thermal equilibrium conditions, this fact reflects the equivalence between different ensembles in statistical mechanics. In this section, we focus on comparing ion currents and 1D multi-ion PMFs with those extracted from all-atom MD simulations and to experimental conductance values. To this end, a more realistic OmpC system setup is employed (see Fig. 3) and atomic thermal fluctuations are properly adjusted to reproduce those observed in MD simulations. Further details on the system setup, simulation types regarding the atomic thermal fluctuations and BD simulation protocol are given in Section 7 in the SI. The results for average ion currents and conductances as obtained from flexible BD and all-atom MD simulations are given in Tables 3 and 4. As can be observed, the flexible BD simulations provide very similar conductance values to those extracted from all-atom MD simulations. In turn, the conductance values extracted from BD simulations are in excellent agreement with an experimental value (for one monomer) of 1.05 nS, 20 which was estimated at 300 K from a linear interpolation from values at different temperatures, when using the dielectric screening functions C2 and C3 (see Table S2 in the SI). This finding is an remarkable accomplishment for BD simulations which, in general, tend to overestimate

18


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


Figure 2: Average 1D multi-ion PMFs along the channel axis for 300 mM KCl salt solution at 300 K under equilibrium conditions (Vmp = 0) using different simulation schemes and numerical integration methods of Langevin equations, i.e., the EM, MLS and PC methods. The SHAKE algorithm was used for the holonomic bond distance constraints. Note that the overlap among lines for the different numerical integration methods is a consequence of the excellent convergence achieved in these cases.

the conductance by 30-50%. 3,52 The flexible BD simulations using the screening function C1 slightly underestimate the experimental conductance. This is not surprising because C1 assumes residues significantly buried inside the protein that is not the case for the considered ones. Still, one should keep in mind that a modification of the dielectric constant values for the implicit pore, membrane and solvent regions could also affect the conductance in a nontrivial way. Interestingly, all-atom MD simulations provide asymmetric conductance values with respect to the direction of the applied field. This trend seems to be confirmed by the BD simulations but should be taken with caution since the statistical errors usually overlap for both cases. The only minor discrepancies between BD and MD simulations are observed in the ratio between K+ and Cl− currents. In the BD simulations, the hydrodynamics model for ion diffusion is approximate as based on several assumptions. 13,53 In the MD simulations, the ion diffusion depends critically on the selected water model. 54 Average multi-ion PMFs from flexible BD and all-atom MD simulations are depicted in Fig. 4. These PMFs successfully reproduce the weak cation selectivity of the OmpC porin. Comparing the BD and MD simulations, ions should be able to overcome similar overall energy barriers during the permeation process. This finding is remarkable in view of the important differences between both simulation methodologies. To name a few: (i) the 19



Figure 3: The system setup for the BRODEA simulations is mainly shown in surface representation. Residues belonging to loops L2, L3 and L4 are highlighted in blue, magenta and green, respectively. The loops L2 and L4 from monomers 2 and 3 are shown in van der Waals representation. In addition, residues belonging to turn T1 from monomer 3 are also shown in van der Waals representation and highlighted in red. For turn T1, residues from the associated β-sheets, i.e., β160 and β1, are included and denoted as β160 − T 1 − β1.

Table 3: Average ion currents I and conductance values G for a simulations with a 1 M KCl salt solution at a transmembrane potential of Vmp = 200 mV and at 300 K. Flexible BD simulations using the dielectric screening functions C1 , C2 and C3 (see Table S2 in the SI) are labeled as BD− C1 , BD− C2 and BD− C3 , respectively. Also the results from all-atom MD simulations are given for a single monomer. The predictor-corrector method was used for the integration of Langevin equations. In addition, the errors are shown in terms of the the standard deviation. METHOD BD− C1 BD− C2 BD− C3 MD

IK+ [pA] 118±9 150±29 182±14 122±8

ICl− [pA] 29±6 33±4 30±5 45±6

IK+ /ICl− 4.2±1.0 4.6±1.2 6.2±1.2 2.7±0.4

G [nS] 0.74±0.06 0.92±0.14 1.06±0.07 0.83±0.05

solvent is replaced by a high-dielectric medium; (ii) the lipid bilayer is replaced by layers of dummy atoms; (iii) atomic fluctuations are only partially considered for a reduced set of Brownian particles; and (iv) the Ewald method for electrostatic interactions is replaced by mean-field approximations and cutoff schemes. It has to be emphasized that the energy differences between BD and MD PMFs remain in the same order as the thermal energy and that better statistics are expected from BD simulations considering that the PMFs are adequately estimated by using a very fine grid size of 0.5 Å.

20


Page 20 of 37

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


Table 4: Same as in Table 3 but with opposite voltage of Vmp = −200 mV. METHOD BD− C1 BD− C2 BD− C3 MD

IK+ [pA] -135±12 -161±13 -186±13 -132±14

ICl− [pA] -32±6 -31±3 -32±6 -64±3

IK+ /ICl− 4.4±0.8 5.2±0.9 6.1±1.5 2.1±0.2

G [nS] 0.83±0.06 0.96±0.06 1.09±0.07 0.98±0.07

Figure 4: Average 1D multi-ion PMFs along the channel axis for 1 M KCl salt solution at 300 K without applied field using flexible BD and MD simulations. The predictor-corrector method was employed for the integration of the Langevin equations.

4.3

Translocation of antibiotics molecules

In this section, we discuss the permeation pathways of ciprofloxacin and enrofloxacin molecules through the porin OmpC. The diffusion routes are inferred from 2D FESs obtained by performing temperature accelerated BD (TABD) simulations. Moreover, the findings are compared to all-atom level results extracted from a combination of metadynamics simulations and a zero-temperature string method. 22,23

4.3.1

System setup and simulation details

A monomer of the OmpC porin was considered as the starting structure, which includes additional loops and a turn from the adjacent monomers as described in Section 4.2. The ciprofloxacin or enrofloxacin molecule was placed at around 25 and 20 Å from the constriction region in the extracellular and the periplasmic vestibules, respectively. Three different kind of systems were created for the TABD simulations: (i) a so-called implicit system in which 21



the entire protein was treated implicitly and the antibiotic molecules were treated explicitly; (ii) a rigid explicit system where some protein residues located in the constriction region (the same ones as those described in Section 7 in the SI) were treated explicitly as well but their positions were fixed during the simulations; and (iii) a flexible explicit system where the position of the protein residues were allowed to fluctuate. Most of the simulation details were identical to those described in Section 7 in the SI. For this reason, we mention here only those features that are specific for the TABD simulations. An optimized CGenFF force field was used for describing the bonded and non-bonded parameters of the atoms belonging to the respective antibiotic molecules. 22,23 The electrostatic interactions between Brownian particles were damped using the dielectric screening function C3 (see Table S2 in the SI). Furthermore, the self-diffusion constants for the hydrogen atoms of the antibiotics molecule were manually adjusted to the value of the HB1 atom type of a protein. The CV z, which specifies the location of the antibiotic along the channel axis, was treated as dynamical variable in the TABD simulations. In addition, the CV zij , which describes the orientation of the antibiotic with respect to the channel axis, was monitored during the TABD simulations. A proper definition for these CVs is provided in Section 5 of the SI. Note that the CVs z and zij are equivalent to those employed in previous studies for ciprofloxacin 22 and enrofloxacin, 23 i.e., the same antibiotic atoms were involved. In addition, the antibiotic position was restricted inside the monomer by applying half-harmonic walls with a constant force of 0.25 kcal/mol/Å2 on the x, y and z variables, where x and y were defined in the same manner as the z coordinate based on the respective positions of protein and antibiotic molecule. Several parameters need to be adjusted prior to the TABD simulations. Four different values of the artificial temperature, i.e., 500, 1000, 1500 and 2000 K, were selected for improving the sampling in the z direction. The constant force associated to the CV z (see Eq. 7) was set to 15 kcal/mol/Å2 for simulations carried out at an artificial temperature of

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


500 K and 25 kcal/mol/Å2 for simulations carried out at the rest of the artificial temperatures. The self-diffusion constant associated with the CV z was set to 0.0118 Å2 /ps, which is five time smaller than the lowest diffusion constant found in the explicit atom list. For a given artificial temperature, a total of 20 independent TABD simulations were performed, where the initial position of the ciprofloxacin was located at the extracellular and periplasmic sides for 10 simulations each. For each simulation, 100 BD steps were used for each single step in CV space. The simulations were performed for 2 µs each which corresponds to a total simulation time of 40 µs for each artificial temperature. The values of the CVs z and zij were stored every 10 ps, which results in 4 million data points for each artificial temperature. Normalized histograms N (z, zij ) were estimated as a function of the CVs z and zij with a grid size of 0.5 and 0.4 Å, respectively. The FESs were computed by means of the expression 55,56 F (z, zij ) = −kB T¯ln(N (z, zij )) + C, where the constant C was adjusted such that F (z, zij ) = 0 at z = −25 Å and zij = 8 Å. 4.3.2

Free energy landscapes

The FESs for ciprofloxacin translocation obtained from the TABD simulations are shown in Fig. 5. The CV z describes the position of the ciprofloxacin molecule along the OmpC channel axis. The constriction region is located at z ∈ [−5, 5] Å between the extracellular (z < −5 Å) and the periplasmic (z > 5 Å) vestibules. Information concerning the orientation of the ciprofloxacin molecule with respect to the OmpC channel axis is provided by the CV zij . The limiting values of −8 and 8 Å correspond to those orientations that are parallel and anti-parallel to the channel axis. The carboxyl group of the ciprofloxacin is facing ahead in the z direction for zij = −8 Å, whereas the amino group is facing ahead in the z direction for zij = 8 Å. Moreover, the intermediate value zij = 0 Å corresponds to those orientations that are perpendicular to the channel axis. As shown in Fig. 5, the ciprofloxacin molecule enters into the extracellular mouth predominantly with its amino group ahead pointing towards the constriction region. In turn, the ciprofloxacin molecule exits from the

23



Figure 5: 2D FESs as function of the CVs z and zij from TABD simulations. The FESs for ciprofloxacin permeation through OmpC porin are shown for three different systems, i.e., an implicit, a rigid explicit and a flexible explicit setup, at two temperatures, i.e., 500 K (A), 1000 K (B). (C) Reweighted free energy landscapes as a function of the CVs z and zij from multiple walker WTmetaD simulation performed with different bias deposition rate β. The lowest-energy translocation pathway is depicted by a black line. (D) Three panels illustrating representative ciprofloxacin conformations along the lowest-energy translocation pathway as the molecule permeates from the EC to PP side. The ciprofloxacin molecule is shown in stick representation using different colors and the carboxyl and amino groups are highlighted in red and blue, respectively. These results shown in C and D have been extracted from our previous work. 22

periplasmic mouth with its carboxyl group facing ahead in the z direction. Therefore, a reorientation of the antibiotic is required inside the monomer. Two main diffusion pathways 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


can be seen: (i) a pathway that requires a reorientation of the ciprofloxacin molecule in the extracellular vestibule before crossing the constriction region with its carboxyl group ahead and (ii) a pathway that requires such reorientation in the periplasmic vestibule after crossing the constriction region. Pathway (i) is the energetically most favorable since it is the only pathway available at the lowest artificial temperature, i.e., at 500 K. This fact is in very good agreement with the results obtained from a computational protocol that combines all-atom metadynamics and a zero-temperature string method shown in Fig. 5C. 22 It has to be noted that already the metadynamics simulation do show clearly visible variation of the FESs when varying the steering parameter β. For a more detailed discussion on these variations we refer the interested reader to Ref. 22. The major FES discrepancies between the implicit and explicit systems are mainly located inside or near the constriction region, i.e., at positions where the electrostatic interactions are handled differently and become more pronounced when increasing the artificial temperature. In case of the implicit system, the pathway (i) remains the most favorable one at the full range of artificial temperatures since their related FESs do not undergo dramatic changes. In the case of the explicit systems, pathway (ii) becomes accessible for artificial temperatures above 500 K (see results for 1500 and 2000 K in the SI). Interestingly, the implicit system provides overall FESs in better agreement with those previously reported. 22 This fact seems to suggest an opposite trend as compared to the ion permeation 13 but further research has to be done in order to determine the effect of different partitions of pore atoms between (fixed) implicit and explicit ones (see Section 5). The above mentioned all-atom metadynamics combined with a zero-temperature string method has also been performed for the translocation of enrofloxacin through OmpC. 23 Thus, we tested also this system with the present approach. We restrained ourselves, however, to those schemes which seemed to be most promising in the case of ciprofloxacin translocation, i.e., a fully implicit description of the pore and artificial temperatures of 500 and 1000 K. The respective results are shown in Fig. 6. Again this antibiotic molecule has to change

25



Figure 6: 2D FESs as function of the CVs z and zij from TABD simulations for the enrofloxacin permeation through the porin OmpC obtained using the implicit model are shown at the two temperatures, i.e., 500 K (A) and 1000 K (B). (C) Reweighted free energy landscapes as a function of the CVs z and zij from multiple walker WTmetaD simulation with lowest-energy translocation pathway depicted as a black line. These results have been extracted for our previous work. 23

orientation during the translocation process as proposed by the all-atom simulations and reproduced by the TABD calculations. As clearly visible in the 1000 K TABD results, a higher barrier in the constriction zone close to z = 0 has to be overcome in case of enrofloxacin compare to ciprofloxacin. This finding is in very good qualitative agreement with the all-atom metadynamics simulations as well as experiments. 23 In summary, we demonstrate that the TABD method is able to provide qualitative predictions for the substrate diffusion routes across nanopores. Several issues need to be addressed in future studies. First, a certain degree of flexibility could be required for most of the water exposed residues along the channel axis with the aim of quantitatively improving the FES estimations. However, this may require an optimization of the code as the CPU time increases drastically with the number of Brownian particles (see Section 5). Second, the influence of various physiologically relevant ionic salts can modify the permeability significantly. 57 Third, the influence of a transmembrane potential might be relevant in the diffusion process. 57 Fortunately, these last two issues can be handled adequately using the BRODEA code. An interesting point is of course the speedup gained by using the BD instead of all-atom

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


Table 5: Core hours on an Intel Xeon CPU ES-2640 (2.50 GHz) for a simulation time of 4 µs to obtain the respective FESs. The numbers within the brackets give the speedup compare to the all-atom calculations. all-atom metadynamics ciprofloxacin 120810 enrofloxacin 124097

implicit BD 24 (5110) 33 (3760)

explicit rigid BD explicit flexible BD 188 (640) 521 (230)

MD simulations to obtain FESs. Table 5 list the number of core hours for a simulation time of 4 µs. Note that the MD simulations are highly parallelized and the mentioned calculations were performed on 192 cores in parallel. Especially for the implicit BD version in which the antibiotics molecules are still treated on an all-atom level, the speedup is several thousand-fold. Details in the description are lost but the computational gain is impressive while yielding the same qualitative outcome.

5

Summary and outlook

The results in this study are based on a new approach for including explicit atoms into BD simulations called BRODEA. This approach can be employed for studying ion permeation and substrate translocation across nanopores. In order to achieve this aim, several schemes were employed by combining different boundary conditions and electrostatic models. In addition, a temperature accelerated BD method was used for sampling free energy landscape and characterizing substrate diffusion routes across a nanopore confinement. The equivalence between the different schemes, the integration methods for solving the Langevin equations and the algorithms for applying bond distance constraints was demonstrated. BRODEA simulations were able to reproduce the experimental conductance at high ion concentration yielding values similar to those of previous all-atom MD simulations. In addition, BD simulations were able to yield reasonable conductance estimates at low ion concentrations for which MD simulations are computationally very demanding. In spite of some relevant differences between the BD and MD methodologies, the BRODEA simulations were able to

27



provide free ion energy profiles that were directly comparable to those from MD simulations. Moreover and notably, BRODEA simulations were able to characterize the ciprofloxacin and enrofloxacin permeation pathways in qualitative close agreement with a recent studies that combine metadynamics and a zero-temperature string method. To our knowledge, this is the first time such antibiotic permeation pathways has been characterized by a technique based on Brownian dynamics. Despite all these positive results, one should not forget that BD simulations contain considerable approximation when compared to MD simulations. Thus, one cannot expect a quantitative agreement between both methods, e.g., for FESs. Nevertheless, a good qualitative agreement can be very helpful. In the case of antibiotics translocation, already the knowledge of orientation during the translocation process can be extremely valuable when considering chemical modifications which should not decrease or even increase this translocation process. Moreover, relative barriers like the once for ciprofloxacin and enrofloxacin are of interest to get an insight into the transport process. In this aspect, the results for a fully implicit description of the nanopore combined with an explicit description of the antibiotic has shown to be powerful and computationally very efficient. This technique would actually allow for a first fast screening of many compounds potentially followed by a more detailed atomistic analysis for the best hits, i.e., the fastest translocating molecules. A couple of points for improving the BRODEA performance can be mentioned. On the one hand, a more systematic study of how to determine which nanopore atoms that should be considered explicitly would be helpful. As a primary rule, Brownian particles have to be surrounded by water and display significant dynamics in the relevant time scales. However, the modification of the set of explicit atoms could affect the conductance and other properties in a non-trivial way since two different conceptual frameworks are used for the treatment of the electrostatic effects, i.e., different macroscopic media for implicit regions separated by discontinuous boundaries and characterized by specific dielectric constants in combination with electrostatic interactions between Brownian particles that are mainly damped by nonlinear distance-dependent screening functions of sigmoidal form. On the other hand, the

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


computation of the non-bonded forces increases critically as a function of the total number of Brownian particles. Several algorithms has been proposed to speed-up the force and minimum image distance computations, e.g., the linked-list cell method, 58–61 ghost atoms approach 61 and techniques for splitting up the neighbor lists , 62,63 which can also be incorporated into the MD-BD framework. Finally, a parallelization of the BRODEA code ought to be an important goal in the near future.

Acknowledgement We gratefully acknowledge Sergei Yu. Noskov for helpful discussions. The research leading to these results was conducted as part of the Translocation consortium (www.translocation.eu) and has received support from the Innovative Medicines Joint Undertaking under Grant Agreement No. 115525, resources which are composed of financial contribution from the European Union’s seventh framework programme (FP7/2007- 2013) and EFPIA companies in kind contribution.

Supporting Information Available Supporting material is available on the details of the stochastic differential equations, numerical integration schemes, the LINCS algorithm as well as on details of the BD and MD simulations. In addition, a movie of a flexible BD simulation with a 1 M KCl salt solution at a transmembrane potential of -200 mV is accessible ( K+ ions in green, Cl− ions in magenta). This information is available free of charge via the Internet at http://pubs.acs.org.

References (1) Majd, S.; Yusko, E. C.; Billeh, Y. N.; Macrae, M. X.; Yang, J.; Mayer, M. Applications of Biological Pores in Nanomedicine, Sensing, and Nanoelectronics. Curr. Opin. 29



Biotechnol. 2010, 21, 439–476. (2) Kuyucak, S.; Andersen, O. S.; Chung, S.-H. Models of Permeation in Ion Channels. Rep. Prog. Phys. 2001, 64, 1427–1472. (3) Roux, B.; Allen, T.; Bernéche, S.; Im, W. Theoretical and Computational Models of Biological Ion Channels. Q. Rev. Biophys. 2004, 37, 15–103. (4) 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. (5) Maffeo, C.; Bhattacharya, S.; Yoo, J.; Wells, D.; Aksimentiev, A. Modeling and Simulation of Ion Channels. Chem. Rev. 2012, 112, 6250–6284. (6) Corry, B.; Thomas, M. Mechanism of Ion Permeation and Selectivity in a Voltage Gated Sodium Channel. J. Am. Chem. Soc. 2012, 134, 1840–1846. (7) Gordon, D.; Chung, S.-H. Extension of Brownian Dynamics for Studying Blockers of Ion Channels. J. Phys. Chem. B 2012, 116, 14285–14294. (8) Berti, C.; Furini, S.; Gillespie, D.; Boda, D.; Eisenberg, R. S.; Sangiorgi, E.; Fiegna, C. Three-Dimensional Brownian Dynamics Simulator for the Study of Ion Permeation through Membrane Pores. J. Chem. Theory Comput. 2014, 10, 2911–2926. (9) Parkin, J.; Chavent, M.; Khalid, S. Molecular Simulations of Gram-Negative Bacterial Membranes: A Vignette of Some Recent Successes. Biophys. J. 2015, 109, 461–468. (10) Pohorille, A.; Wilson, M. A.; Wei, C. Validity of the Electrodiffusion Model for Calculating Conductance of Simple Ion Channels. J. Phys. Chem. B 2017, 121, 3607–3819. (11) Kopec, W.; Köpfer, D. A.; Vickery, O. N.; Bondarenko, A. S.; Jansen, T. L. C.; de Groot, B. L.; Zachariae, U. Direct Knock-On of Desolvated Ions Governs Strict Ion Selectivity in K+ Channels. Nat. Chem. 2018, 10, 813. 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


(12) Pothula, K. R.; Solano, C. J.; Kleinekathöfer, U. Simulations of Outer Membrane Channels and Their Permeability. Biochimica et Biophysica Acta (BBA) - Biomembranes 2016, 1858, 1760–1771. (13) Solano, C. J. F.; Pothula, K. R.; Prajapati, J. D.; De Biase, P. M.; Noskov, S. Y.; Kleinekathöfer, U. BROMOCEA Code:

An Improved Grand Canonical Monte

Carlo/Brownian Dynamics Algorithm Including Explicit Atoms. J. Chem. Theory Comput. 2016, 12, 2401–2417. (14) Im, W.; Seefeld, S.; Roux, B. A Grand Canonical Monte Carlo-Brownian Dynamics Algorithm for Simulating Ion Channels. Biophys. J. 2000, 79, 788–801. (15) Maragliano, L.; Vanden-Eijnden, E. A Temperature Accelerated Method for Sampling Free Energy and Determining Reaction Pathways in Rare Events Simulations. Chem. Phys. Lett. 2006, 426, 168–175. (16) Benz, R. Structure and Function of Porins from Gram-Negative Bacteria. Annu. Rev. Microbiol. 1988, 42, 359–393. (17) Baslé, A.; Rummel, G.; Storici, P.; Rosenbusch, J. P.; Schirmer, T. Crystal Structure of Osmoporin OmpC from E. coli at 2.0 Å. J. Mol. Biol. 2006, 362, 933–942. (18) Masi, M.; Pagès, J.-M. Structure, Function and Regulation of Outer Membrane Proteins Involved in Drug Transport in Enterobactericeae: The OmpF/C–TolC Case. Open Microbiol. J. 2013, 7, 22. (19) Yamashita, E.; Zhalnina, M. V.; Zakharov, S. D.; Sharma, O.; Cramer, W. A. Crystal Structures of the OmpF Porin: Function in a Colicin Translocon. EMBO J. 2008, 27, 2171–2180. (20) Biro, I.; Pezeshki, S.; Weingart, H.; Winterhalter, M.; Kleinekathöfer, U. Comparing the

31



Temperature-Dependent Conductance of the Two Structurally Similar E. coli Porins OmpC and OmpF. Biophys. J. 2010, 98, 1830–1839. (21) Mahendran, K. R.; Kreir, M.; Weingart, H.; Fertig, N.; Winterhalter, M. Permeation of Antibiotics through Escherichia coli OmpF and OmpC Porins Screening for Influx on a Single-Molecule Level. J. Biomol. Screen. 2010, 15, 302–307. (22) Prajapati, J. D.; Solano, C. J. F.; Winterhalter, M.; Kleinekathöfer, U. Characterization of Ciprofloxacin Permeation Pathways across the Porin OmpC Using Metadynamics and a String Method. J. Chem. Theory Comput. 2017, 13, 4553–4566. (23) Prajapati, J. D.; Solano, C. J. F.; Winterhalter, M.; Kleinekathöfer, U. Enrofloxacin Permeation Pathways across the Porin OmpC. J. Phys. Chem. B 2018, 122, 1417–1426. (24) van Kampen, N. Stochastic Processes in Physics and Chemistry, 2nd ed.; NorthHolland, 1992. (25) Pohorille, A.; Jarzynski, C.; Chipot, C. Good Practices in Free-Energy Calculations. J. Phys. Chem. B 2010, 114, 10235–10253. (26) Bussi, G.; Branduardi, D. Reviews in Computational Chemistry; Wiley-Blackwell, 2015; Vol. 28; pp 1–49. (27) Murphy, T. J.; Aguirre, J. L. Brownian Motion of N Interacting Particles. I. Extension of the Einstein Diffusion Relation to the N -Particle Case. J. Chem. Phys. 1972, 57, 2098–2104. (28) Yamakawa, H.; Tanaka, G.; Stockmayer, W. H. Correlation Function Formalism for the Intrinsic Viscosity of Polymers. J. Chem. Phys. 1974, 61, 4535–4539. (29) Wilemski, G. On the Derivation of Smoluchowski Equations with Corrections in the Classical Theory of Brownian Motion. J. Stat. Phys. 1976, 14, 153–169.

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


(30) Gray, P. Projection Operators in Classical Kinetic Theory. Mol. Phys. 1971, 21, 675– 687. (31) Falkenhagen, H.; Ebeling, W. Statistical derivation of diffusion equations according to the Zwanzig method. Phys. Lett. 1965, 15, 131–132. (32) Zwanzig, R. Nonequilibrium Statistical Mechanics; Oxford University Press, 2001. (33) Lax, M. Classical Noise Iv: Langevin Methods. Rev. Mod. Phys. 1966, 38, 541. (34) Hess, W.; Klein, R. Dynamical Properties of Colloidal Systems: 1. Derivation of Stochastic Transport Equations. Physica A 1978, 94, 71–90. (35) Tough, R. J. A.; Pusey, P. N.; Lekkerkerker, H. N. W.; Van Den Broeck, C. Stochastic Descriptions of the Dynamics of Interacting Brownian Particles. Mol. Phys. 1986, 59, 595–619. (36) Gardiner, C. W. Handbook of Stochastic Methods, 4th ed.; Springer, 2009. (37) A. D. MacKerell, J.; Feig, M.; Brooks, C. L. Improved Treatment of the Protein Backbone in Empirical Force Fields. J. Am. Chem. Soc. 2004, 126, 698. (38) A. D. MacKerell, J.; Feig, M.; Brooks, C. L. Extending the Treatment of Backbone Energetics in Protein Force Fields: Limitations of Gas-phase Quantum Mechanics in Reproducing Protein Conformational Distributions in Molecular Dynamics Simulations. J. Comp. Chem. 2004, 25, 1400. (39) Brooks, B. R.; Brooks, C. L.; Mackerell, A. D.; Nilsson, L.; Petrella, R. J.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch, S.; Caflisch, A.; Caves, L.; Cui, Q.; Dinner, A. R.; Feig, M.; Fischer, S.; Gao, J.; Hodoscek, M.; Im, W.; Kuczera, K.; Lazaridis, T.; Ma, J.; Ovchinnikov, V.; Paci, E.; Pastor, R. W.; Post, C. B.; Pu, J. Z.; Schaefer, M.; Tidor, B.; Venable, R. M.; Woodcock, H. L.; Wu, X.; Yang, W.;

33



York, D. M.; Karplus, M. CHARMM: The Biomolecular Simulation Program. J. Comp. Chem. 2009, 30, 1545–1614. (40) Im, W.; Roux, B. Ion Permeation and Selectivity of OmpF Porin: A Theoretical Study Based on Molecular Dynamics, Brownian Dynamics, and Continuum Electrodiffusion Theory. J. Mol. Biol. 2002, 322, 851–869. (41) Mehler, E. L.; Guarnieri, F. Biophys. J. 1999, 75, 3. (42) Hassan, S. A.; Guarnieri, F.; Mehler, E. L. A General Treatment of Solvent Effects Based on Screened Coulomb Potentials. J. Phys. Chem. B 2000, 104, 6478–6489. (43) Hassan, S. A.; Mehler, E. L. A Critical Analysis of Continuum Electrostatics: The Screened Coulomb Potential-Implicit Solvent Model and the Study of the Alanine Dipeptide and Discrimination of Misfolded Structures of Proteins. Proteins 2002, 47, 45–61. (44) Li, X.; Hassan, S. A.; Mehler, E. L. Long dynamics simulations of proteins using atomistic force fields and a continuum representation of solvent effects: calculation of structural and dynamic properties. Proteins 2005, 60, 464–484. (45) Corry, B.; Hoyles, M.; Allen, T. W.; Walker, M.; Kuyucak, S.; Chung, S.-H. Reservoir Boundaries in Brownian Dynamics Simulations of Ion Channels. Biophys. J. 2002, 82, 1975–1984. (46) Berti, C.; Furini, S.; Gillespie, D. PACO: PArticle COunting Method to Enforce Concentrations in Dynamic Simulations. J. Chem. Theory Comput. 2016, 12, 925–929. (47) Im, W.; Roux, B. Brownian Dynamics Solutions of Ions Channels: A General Treatment of Electrostatic Reaction Fields for Molecular Pores of Arbitrary Geometry. J. Chem. Phys. 2001, 115, 4850.

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


(48) Gumbart, J.; Khalili Araghi, F.; Sotomayor, M.; Roux, B. Constant Electric Field Simulations of the Membrane Potential Illustrated with Simple Systems. Biochim. Biophys. Acta, Biomembr. 2012, 1818, 294–302. (49) Melcr, J.; Bonhenry, D.; Timr, S.; Jungwirth, P. Transmembrane Potential Modeling: Comparison between Methods of Constant Electric Field and Ion Imbalance. J. Chem. Theory Comput. 2016, 12, 2418–2425. (50) Verlet, L. Computer "Experiments" on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules. Phys. Rev. 1967, 159, 98. (51) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS: A Linear Constraint Solver for Molecular Simulations. J. Comput. Chem. 1997, 18, 1463–1472. (52) Zhekova, H. R.; Ngo, V.; da Silva, M. C.; Salahub, D.; Noskov, S. Selective Ion Binding and Transport by Membrane Proteins–A Computational Perspective. Coord. Chem. Rev. 2017, (53) Noskov, S.; Im, W.; Roux, B. Ion Permeation through the α-Hemolysin Channel: Theoretical Studies Based on Brownian Dynamics and Poisson-Nernst-Plank Electrodiffusion Theory. Biophys. J. 2004, 87, 2299–2309. (54) Modi, N.;

Singh, P. R.;

Mahendran, K. R.;

Schulz, R.;

Winterhalter, M.;

Kleinekathöfer, U. Probing the Transport of Ionic Liquids in Aqueous Solution Through Nanopores. J. Phys. Chem. Lett. 2011, 2, 2331–2336. (55) Cuendet, M. A.; Tuckerman, M. E. Free Energy Reconstruction from Metadynamics or Adiabatic Free Energy Dynamics Simulations. J. Chem. Theory Comput. 2014, 10, 2975–2986. (56) Awasthi, S.; Nair, N. N. Exploring High Dimensional Free Energy Landscapes: Temperature Accelerated Sliced Sampling. J. Chem. Phys. 2017, 146, 094108. 35



(57) Bhamidimarri, S. P.; Prajapati, J. D.; van den Berg, B.; Winterhalter, M.; Kleinekathöfer, U. Role of Electroosmosis in the Permeation of Neutral Molecules: CymA and Cyclodextrin As an Example. Biophys. J. 2016, 110, 600–611. (58) Mattson, W.; Rice, B. M. Near-Neighbor Calculations Using a Modified Cell-linked List Method. Comput. Phys. Commun. 1999, 119, 135–148. (59) Yao, Z.; Wang, J.-S.; Liu, G.-R.; Cheng, M. Improved Neighbor List Algorithm in Molecular Simulations Using Cell Decomposition and Data Sorting Method. Comput. Phys. Commun. 2004, 161, 27–35. (60) Heinz, T. N.; Hünenberger, P. H. A Fast Pairlist-Construction Algorithm for Molecular Simulations under Periodic Boundary Conditions. J. Comput. Chem. 2004, 25, 1474– 1486. (61) Gonnet, P. A Simple Algorithm to Accelerate the Computation of Non-Bonded Interactions in Cell-Based Molecular Dynamics Simulations. J. Comput. Chem. 2007, 28, 570–573. (62) van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. GROMACS: Fast, Flexible, and Free. J. Comput. Chem. 2005, 26, 1701–1718. (63) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435–447.

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


Graphical TOC Entry

37


Brownian Dynamics Approach Including Explicit Atoms for Studying

Recommend Documents