BROMOCEA Code: An Improved Grand Canonical Monte Carlo

Subscriber access provided by MONASH UNIVERSITY

Article

The BROMOCEA code: an improved grand canonical Monte Carlo/Brownian dynamics algorithm including explicit atoms Carlos Jose Fernandez Solano, Karunakar Reddy Pothula , Jigneshkumar Dahyabhai Prajapati, Pablo Martin De Biase, Sergei Yu. Noskov, and Ulrich Kleinekathöfer J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.5b01196 • Publication Date (Web): 18 Apr 2016 Downloaded from http://pubs.acs.org on April 18, 2016

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

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

Page 1 of 51

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

The BROMOCEA code: an improved grand canonical Monte Carlo/Brownian dynamics algorithm including explicit atoms Carlos J. F. Solano,∗,† Karunakar R. Pothula,† Jigneshkumar D. Prajapati,† Pablo M. De Biase,‡ Sergei Yu. Noskov,‡ and Ulrich Kleinekathöfer∗,† †Department of Physics and Earth Sciences, Jacobs University Bremen, Campus Ring 1, 28759 Bremen, Germany ‡Centre for Molecular Simulation, Department of Biological Sciences, University of Calgary, Calgary, AB, Canada, T2N 1N4 E-mail: [email protected]; [email protected]

Abstract All-atom molecular dynamics simulations have a long history of applications studying ion and substrate permeation across biological and artificial pores. While offering unprecedented insights into the underpinning transport processes, MD simulations are limited in time-scales, ability to simulate physiological membrane potentials or asymmetric salt solutions and require substantial computational power. While several approaches to circumvent all of these limitations were developed, Brownian dynamics simulations remain an attractive option to the field. The main limitation, however, is an apparent lack of protein flexibility important for the accurate description of permeation events. In the present contribution, we report an extension of the Brownian dynamics scheme which includes conformational dynamics. To achieve this goal, the dynamics

1

ACS Paragon Plus Environment


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

of amino-acid residues was incorporated into the many-body potential of mean force and into the Langevin equations of motion. The developed software solution, called BROMOCEA, was applied to ion transport through OmpC as a test case. Compared to fully atomistic simulations, the results show a clear improvement in the ratio of permeating anions and cations. The present tests strongly indicate that pore flexibility can enhance permeation properties which will become even more important in future applications to substrate translocation.

1

Introduction

Protein pores with diameter below 10 nm have attracted great deal of attention in recent years. 1 Biological channels of this size are ubiquitous among all kingdoms of life including baterial porins, mitochondrial voltage-dependent channels or translocation systems. Channels with diameter less than 10 nm are involved in a broad variety of transport processes including antibiotic translocation, toxin activity and selective metabolite transport, to name a few. 2,3 The mid- and wide-sized nanopores can also be manufactured artificially with potential applications in targeted delivery, DNA sequencing and for uses as nanosensors. 4 Therefore, the research directed at the unraveling key principles responsible for selective ion and substrate dynamics in nanoscale confinements is crucial both for human health and development of novel biotechnological systems. Molecular dynamics (MD) simulations have been proven to be an important tool for studies of ion transport. 5–7 Among the investigated membrane pores are porins like OmpF 8,9 and OmpC, 10 more narrow and specific pores like OmpG, 11 OmpA, 12 and OprP 13 and also channels selective for large substrates. 14 These atomic simulations greatly improved the understanding of the ion transport processes and resulted in numerous leads for experimental studies. 14–17 In spite of the great utility of MD simulations in understanding transport processes, the time scales relevant for ion and substrate translocation represent a serious problem. All-atom MD simulations require integration over the fast motions such as vibra2


Page 2 of 51

Page 3 of 51

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


tions of hydrogen atoms. At the same time, direct modeling of a metabolite or antibiotic transport may require micro- to milliseconds simulations. The other conceptual challenge comes from the apparent need in modeling of translocation events in the presence of membrane voltage and/or asymmetric compositions of buffers. 18 Both challenges are conceptually solvable, for example Roux and co-workers reported recently on the extension of periodic boundary conditions to asymmetric buffer solutions, 19 but require a number of assumptions to be in place. A viable alternative can be found in implicit models effectively integrating out degrees of freedom not directly involved into transport process. Coarse-grained bead models are one way to approach larger time-scales inaccessible to all-atom models. 20 For problems such as the transport of ions across membrane pores, Brownian dynamics (BD) simulations have been shown to be an attractive computational alternative for simulating processes over long timescales without having to treat all the degrees of freedom explicitly. 21,22 The partial inclusion of the relevant degrees of freedom due to a protein by, for example, atoms forming narrow permeation pathways, is possible within dynamic Monte Carlo schemes, thus expanding horizons of the method. 23 A computational scheme based on grand canonical Monte Carlo (GCMC) and BD has been developed to simulate the movement of ions in porins. 5,21,24–26 The applicability of GCMC/BD has been validated by previous simulation studies for the outer membrane porin F (OmpF), 25 α-Hemolysin (α-HL) 26,27 and the voltage dependent anion channel (VDAC). 28 Later, an extension of this theory was developed to model ion and DNA translocation across a nanopore confinement driven by an applied electric potential. 29–31 Moreover, using three-dimensional potential of mean force maps determined from all-atom molecular dynamics simulation, Comer and Aksimentiev obtained accurate results from BD simulations close to those from MD simulations. 32 Despite the efficiency of the GCMC/BD approach in studies of wide β-barrel pores, the approach has several major limitations. The solvent is represented by a featureless dielectric medium and therefore the application to narrow pores is challenging. It can be argued that

3



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

ion permeation across a pore confinment cannot be described sufficiently accurate within a mean-field approximation. 33 In narrow channels, significant effects arise from the dynamics of the permeation pathway and the ability to form hydrogen bonds between solvent molecules and a protein. In addition, the protein dielectric constant has a small effect for large and medium-sized pores, where one can approximate the pore water to resemble bulk water. The water in confinements shows molecular structuring, which is dramatically different from that observed in the bulk or wide pores. That is, there are a number of situations where the average structure assumption is no longer valid and one has to model explicitly protein and/or water dynamics in addition to mobile ions. 33–35 Moreover, charged or polar residues lining the pore interact electrostatically with permeating ions and substrates, adapting its conformation to the position of the substrate inside the pore. It is very likely that this dynamic ability to adapt to a permanent conformation will have a major impact onto the dynamics of mobile ions and substrates across the channel. A first approach to include some flexibility back into BD simulations has been performed by Chung and Corry. 36 In this approach, the movements of the carbon and oxygen atoms belonging to carbonyl groups were calculated explicitly in the selectivity filter of the KcsA channel with very promising results. Accordingly, an improved GCMC/BD algorithm was developed, tested and implemented in the new program code “BROMOCEA” (Brownian Dynamics/Monte-Carlo program including explicit atoms). The GCMC/BD software solution is based on the earlier open-source BROMOC code. 29–31 The main goal of the enhanced and significantly altered code is to provide a useful tool for ion permeation and substrate translocation studies, where the substrate and some atoms forming the permeation pathway or confinment, i.e., residues in the constriction zone as well as other residues along the channel surface, can be treated explicitly. This new approach is tested in terms of ion conductance and selectivity across the OmpC channel. The paper is organized as follows. In Section 2, after a brief review of the original

4


Page 4 of 51

Page 5 of 51

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


GCMC/BD algorithm, we introduce the key steps in the improved GCMC/BD algorithm, i.e., the incorporation of explicit atoms into the many-body potential of mean force. In Section 3, we describe some relevant implementation features included into the GCMC/BD algorithm which make an efficient framework with the ability to describe the ion permeation and substrate translocation across membrane channels. In Section 4, we expose a brief review of BROMOCEA. In Section 5, we present the results of computational simulations performed with BROMOCEA on the OmpC porin which we used to test and validate the new features of the code. We conclude with remarks in Section 6.

2

Theory

In this section, an improved GCMC/BD algorithm including internal terms from the CHARMM force field 37 is described. Other force field sets can be included in a similar fashion though details might differ due to slightly different functional forms. In the original GCMC/BD algorithm, 21,24–26 only the mobile ions are explicitly simulated, while the influence from the protein, lipid, and water is implicitly taken into account via effective time-invariant potential and electrostatic maps obtained by solving a modified Poisson-Boltzmann equation 38 (see below). This approximation represents the solvent (water) region as an ion-accessible high dielectric medium and as a source of random diffusive force. The lipid membrane and the channel protein are represented as ion-inaccessible low dielectric regions. In turn, the protein contains a distribution of fixed partial atomic charges and its boundary region is constructed from the molecular structure of the protein. The GCMC algorithm controls the electrochemical potential by stochastically inserting and annihilating ions in the buffer regions that are usually located at the boundaries of the system. The GCMC/BD method rigorously separates the electrostatic forces into static and reaction field contributions for enhanced computational efficiency which allows off-lattice dynamics of the explicit mobile ions and imposes proper thermodynamic boundary conditions on the ion concentrations.

5



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

Page 6 of 51

Assuming that the dynamics of the mobile ions is overdamped due to the high friction in the solvent, a position Langevin equation can be deduced to the Brownian equation of motion, which is in agreement with the process described by the Smoluchowski equation under the assumptions that (i) hydrodynamic interactions (i.e., the effect of one ion on another through the flow of solvent around them) are negligible and that (ii) the ion self-diffusion constants are position-dependent (along the channel axis). 39,40 Thereby, the overdamped Langevin equation is given by dri Di (ri ) =− ∇i W (r1 , . . . , rN ) + ∇i Di (ri ) + ζ i (t), dt kB T

(1)

where ri denotes the position of ion i, kB the Boltzmann constant, T the temperature, W the many-body potential of mean force (PMF), Di (ri ) is the diffusion constant, and ζ i the random Gaussian noise representing water collisions with the mobile ions. Two important properties of this random noise are hζ i (t)i = 0 and hζ i (t) · ζ i (t′ )i = 6Di δ(t − t′ ). A version of the Ermark and McCammon (EM) algorithm 41 obeying the above assumptions is recovered from Eq. 1 when the Euler scheme is applied, i.e., the EM algorithm is a first-order Euler method. The forces acting on the mobile ions are calculated from the spatial gradients of the PMF X N N X φrf (ri ) W (r1 , . . . , rN ) = qi φsf (ri ) + + Ucore (ri ) + uij (rij ), 2 i j>i

(2)

where φsf denotes the electrostatic potential from the protein charges and the transmembrane potential, φrf the reaction field potential arising from ions present in the outer region and the dielectric hetero-junction between solvent, protein, and lipid, while Ucore describes the core-repulsive steric potential from the ion-inaccessible region (protein and lipid). Note that N in Eq. 2 denotes the total number of ions, rij the distance between ions i and j, and qi the (fixed) charge of ion i. The ion-ion potential including the Lennard-Jones (LJ), water-mediated short-range, and dielectric-screened Coulombic electrostatic terms is given 6


Page 7 of 51

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


by uij (rij ) = 4ǫij

"

σij rij

12

−

σij rij

6 #

+ wsr (rij ) +

qi qj , 4πε0 εbulk rij

(3)

where ε0 denotes the vacuum permittivity and εbulk the dielectric constant of bulk water. Moreover, the ion-ion interactions are updated every BD time step while the other potential terms, Ucore , φsf , and φrf , are calculated from look-up tables computed once before the GCMC/BD simulation. In addition, the water-mediated short-range potential is given by

wsr (rij ) = c0 e

c1 −rij c2

cos (c3 (c1 − rij )π) + c4

c1 rij

6

,

(4)

where the coefficients ci are empirically adjusted until a reasonable agreement is obtained between the ion-ion radial distributions from MC simulations using the potential described in Eq. 3 and those from MD simulations. 25 It should be stressed that the number of ions N may be different among successive time steps because of the application of the GCMC algorithm in the buffer regions. Since our main goal is to study ion permeation and substrate translocation where dynamics of the confinement has to be modelled as part of the process, a number of explicit atoms must be properly incorporated into the algorithm. If one assumes that the atomic velocities relax to equilibrium much faster than their positions, i.e., we assume the diffusive regime, and that the hydrodynamics effects between atoms can be neglected, Eq. 1 can be extended to include explicit atoms without additional modifications. On the other hand, the inclusion of atoms in the PMF can be addressed through the incorporation of a force field approach in the description. Henceforth, ions and explicit atoms are interchangeably

7



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

Page 8 of 51

denominated as Brownian particles. The modified PMF can be written as W (r1 , . . . , rN ) =

X

bonds

+

Kb (b − b0 )2 +

X

dihedrals

X

angles

Kθ (θ − θ0 )2 +

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

X

X

U rey−Bradley

impropers

KS (S − S0 )2

Kω (ω − ω0 )2 +

X N N X φrf (ri ) + Ucore (ri ) + uij (rij ). + qi φsf (ri ) + 2 j>i i

X

residues

UCM AP (5)

As shown in this equation, the PMF is now a sum over internal terms for atoms treated explicitly and nonbonded terms for all particles. The internal terms include terms for bonds (b), valence angles (θ), Urey-Bradley (S), dihedral angles (ϕ), improper angles (ω), and CMAP backbone torsional corrections. 42,43 These terms are the same as the internal terms, e.g., in the CHARMM force field. 37 The nonbonded terms include the terms shown in Eq. 2 but in Eq. 5 the summation also includes explicit atoms since N is the total number of Brownian particles.

3

Implementation details

In this section, some relevant improvements included into the GCMC/BD algorithm are described. For the sake of completeness, additional implementation details are also described in the Supplementary Information.

3.1

CHARMM force field for explicit atoms

The internal terms from the CHARMM force field 37 have been implemented in BROMOCEA code. All internal terms are taken to be harmonic, except the dihedral angle term, which is given by a sinusoidal expression. The all-atom implementations of the CHARMM force field include all possible valence and dihedral angles for bonded atoms, and the dihedral angle terms concerning a given bond may be expanded in a Fourier series of up to six terms.

8


Page 9 of 51

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


In addition, the CMAP backbone torsional correction 42,43 has been implemented. Both the Urey-Bradley and improper dihedral terms are used to optimize the fit to vibrational spectra and out-of-plane motions. While the improper dihedral term is used very generally in the CHARMM force field, the Urey-Bradley term tends to be used only in special cases. In principle, the code can also be used to include united atom force fields though this possibility is not detailed here.

3.2

Distance-dependent dielectric constant

As shown in Eq. 3, the particle-particle potential includes an electrostatic term which is dielectrically screened using the dielectric constant of bulk walter εbulk . In the improved GCMC/BD algorithm, a distance-dependent dielectric constant ε(rij ) is available instead of εbulk . Thereby, one can simultaneously model two effects: a substantially-damped electrostatic field when the i and j particles are separated by moderate distances, and a diminished dielectric screening, approaching toward vacuum values, when the particles are in close proximity. A sigmodial form for ε(rij ) has been incorporated into a procedure to calculate pH-dependent properties in proteins with excellent accuracy. 44 This sigmoidal form is now also implemented in the BROMOCEA code.

3.3

The atom-based force-switching method

The nonbonded force inserted on particle i by particle j can be written as

nb fij,α = 24ǫij

σij rij

6 "

1−2

σij rij

6 #

rij,α rij2

6 ) cos (c3 (c1 − rij )π) c4 c1 rij,α c3 π sin (c3 (c1 − rij )π) − −6 + c0 e c2 rij rij rij 1 rij dε(rij ) qi qj rij,α , + − 2 4πε0 ε(rij )rij rij ε(rij ) drij (

c1 −rij c2

9


(6)


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

Page 10 of 51

where the first term corresponds to the force due to the LJ interactions, the second term matches the force due to the water-mediated short-range potential and the third term corresponds to the force due to the Coulombic interactions. Note that rij,α = rj,α − ri,α with α = {x, y, z} are the Cartesian coordinates for the interparticle vector. It should be stressed that the reaction field potential φrf includes the long range effect of the ions present in the outer region (i.e., the region outside of the simulation box where the ions are treated implicitly). To gain computational efficiency, the nonbonded interactions between Brownian particles can be turned off at a cutoff distance to reduce the number of interacting pairs. The atombased force-switching method 45 is known to fully account for short-range, damping longrange component of forces in the system monotonically to zero in the interval between ron and rof f . In the BROMOCEA code, the pairwise nonbonded force is given by nb,sw nb fij,α = fij,α S(rij ),

where S(rij ) =

   1,  

2 −r 2 )2 (r 2 +2r 2 −3r 2 ) (rof on ij ij f of f 2 −r 2 )3 (rof on f

(7)

rij ≤ ron

(8)

, ron < rij ≤ rof f

is a switching function which satisfies S(ron ) = 1 and S(rof f ) = 0.

3.4

Alternative core-repulsive potential scheme

In the original GCMC/BD algorithm, 21,24–26 the core-repulsive potential preventing corecore overlap with the channel and membrane is calculated on a discrete grid and stored for computational efficiency. The entire simulation region is scanned and a large positive energy is assigned to the potential if the distance to a channel atom is less than the combined LJ radius. However, this potential exhibits abrupt variations which may introduce unrealistic particle displacements.

10


Page 11 of 51

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


An alternative repulsive potential scheme has been implemented in the BROMOCEA code. In this scheme, the core-repulsive potential from a pre-computed look-up table is replaced by on the fly LJ potentials where fixed channel atoms are explicitly considered. Although this alternative scheme is more computational demanding than the original one, LJ interactions between particles and fixed channel atoms may be turned off at a cutoff distance using the previous atom-based force-switching method so that computational efficiency is gained with respect to a non-cutoff distance scheme. It is worth mentioning that the alternative core-repulsive scheme leads to a significant improvement compared to the original scheme. Indeed, the core repulsive potentials from the original scheme are replaced by LJ potentials which include attractive dispersion terms in addition to the repulsive ones. The modeling the implicit membrane in the alternative core-repulsive potential scheme needs to be done carefully. The membrane is mimicked by two planes at both sides of the membrane not allowing any atoms to enter the region between them. Each plane is composed of identical closely packed dummy atoms. Care needs to be taken so that no dummy atoms overlap with the channel protein.

3.5

Second-order stochastic Runge-Kutta algorithm for the Brownian dynamics propagation

The EM algorithm 41 for the position Langevin equation defined in Eq. 1 is given by

ri (t + ∆t) = ri (t) + ∆ri (∆t),

(9)

where ∆ri (∆t) = ∆rsi (∆t) + and ∆rsi (∆t) = ∇i Di (ri (t))∆t +

p

2Di (ri (t))∆tη i

Di (ri (t)) Fi (r1 (t), . . . , rN (t))∆t . kB T

11


(10)

(11)


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

Page 12 of 51

In this equation Fi (r1 (t), . . . , rN (t)) = −∇i W (r1 (t), . . . , rN (t)) denotes the systematic forces on the particle i and η i a vector of the random displacement on the particle i sampled from a Gaussian distribution with zero mean (i.e, hη i i = 0) and variance obeying hηi,α ηj,β i = δij δαβ with α, β = {x, y, z}. A second-order stochastic Runge-Kutta (SRK) algorithm can be developed for positiondependent self-diffusion constants. The SRK method updates the particle positions according to Eq. 9 with auxiliary positions rs,a and rs,b where s,b ∆rs,a i + ∆ri ∆ri (∆t) = + 2

p p 2Di (rai )∆t + 2Di (rbi )∆t ηi. 2

(12)

The calculation of the displacement takes place at two points in time. First, one calculates a ∆rs,a i = ∇i Di (ri )∆t +

Di (rai ) Fi (ra1 , . . . , raN )∆t kB T

(13)

where rai = ri (t). Subsequently, a corrector step given by rbi = ri (t) + ∆rs,a i +

p 2Di (rai )∆tη i

(14)

is employed where the random displacement η i is the same as in Eq. 12. Finally, one calculates b ∆rs,b i = ∇i Di (ri )∆t +

Di (rbi ) Fi (rb1 , . . . , rbN )∆t . kB T

(15)

In case of position-independent self-diffusion constants, this algorithm is equivalent to the one described in Ref. 46.

3.6

Constraints

There are two main reasons for using bond-length constraints in simulations. First, removing the rapid vibrations by constraining the bonds allows to increase the simulation time

12


Page 13 of 51

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


step, and thus a reduction of the computational effort for obtaining a trajectory of a given length. Second, treating bonds as constraints is probably a better approximation of their quantum-mechanical behavior than treating them as classical harmonic oscillators 47 since high-vibrational bond-streching frequencies need to be treated quantum mechanically at room temperature. The SHAKE method 48 has been adapted for using in the BD framework (see the Supplementary Information) and has been implemented in the BROMOCEA code. Bond-angle constraints are not implemented since there are several reasons why applying bond-length and bond-angle constraints at the same time are not recommended. 49 On the one hand, bond-angle constraints may alter the dihedral transition rates since bond angle frequencies may be close to dihedral frequencies. On the other hand, the application of the SHAKE algorithm does not improve the computational efficiency when bond angles are also constrained because the strong interdependence of the constraints slows down the convergence of the SHAKE algorithm.

4

The BROMOCEA code

The BROMOCEA code has been implemented based on the BROMOC-D code 29 which was developed to model ion and DNA translocation across a nanopore confinement under an applied electric potential. Thereby, BROMOCEA includes the functionalities and modules previously implemented in BROMOC-D but many of them have been improved. Henceforth, we focus the discussion on the new functionalities and modules implemented in BROMOCEA since DNA translocations using a DNA coarse-grain polymer model 50 has been extensively discussed. 29–31 A typical research project with BROMOCEA can be described in very general terms based on the information flow in the program, which is schematically illustrated in Fig. 1. Before running a BROMOCEA simulation, the parameters for electrostatic potential φsf , reaction field potential φrf , and core-repulsive steric potentital Ucore need to be calculated on a rectan-

13



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

Page 14 of 51

gular grid and stored in (binary) files. Note that Ucore should be calculated on a rectangular grid and stored only when the original core-repulsive potential scheme is used. The BROMOCEA code supports the binary files generated by the Poisson-Boltzmann (PB)/ PoissonNernst-Planck (PNP) program (available at http://thallium.bsd.uchicago.edu/RouxLab/pbeq.html) for these potentials. 25 The user begins a BROMOCEA project by providing the input files required for the simulation. For the explicit atom description which employs the CHARMM force field, one should provide the force field parameters file (PRM), the protein structure file (PSF) for topology/connectivity and the coordinates from a DCD file. The BROMOCEA code allows for describing explicit atoms using several chains (or frames). Thereby, the user can provide PSF and DCD files for each chain. The explicit atoms belonging to the pore are joined to the rigid region through atoms at fixed position. BROMOCEA supports any set of force field parameters properly defined so that extended atoms and/or other coarse-grained particles could be incorporated into this set. The description of the ions is included in the general BROMOCEA input file which, in turn, describes the rest of the parameters necessaries for a proper simulation, selects the functionalities and modules required for the simulation, and specifies the output files and the output parameters to be printed in the general output file. Once the input files are provided, the project enters the production stage, during which internal CHARMM parameters, initial coordinates and the bonded list (and exclude list) are generated, and the BD and GCMC algorithms can be invoked several times in the trajectory builder. BROMOCEA allows the simulation of ion channels with a realistic implementation of boundary conditions of concentration and transmembrane potential under equilibrium and non-equilibrium conditions of ion diffusion and permeation.

5

Modelling of ion permeation across biological pores

In this section, some computational illustrations are described where the main outputs from BROMOCEA simulations, i.e., conductance and one-dimensional (1D) multi-ion PMFs from

14


Page 15 of 51


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 ACS Paragon Plus Environment


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

5.1

Page 16 of 51

Diffusion models

In the GCMC/BD algorithm, the self-diffusion constants for ions and explicit atoms are part of the input to a simulation. It represents diffusive or frictional effects in a dilute system. While ion diffusion constants in bulk solutions may be estimated, there exists no generally accepted model for space-dependent ion diffusion constants in the vicinity of or inside the pore. Two different ion space-dependent diffusion models were employed in this work. The single-valued (SV) diffusion model assumes a uniform ion diffusion constant inside the pore, whose value is the bulk value scaled by a factor 0 < f ≤ 1, that is connected to the bulk value using a switching function. 25 As in previous works, 25,51 this scaling factor is set to f = 0.5 (see Fig. 2). The hydrodynamics (HD) based model of ion diffusion is derived from the assumption of a rigid sphere moving through a rigid cylindrical pore 26 (see Fig. 2). If one assumes that the channel is aligned along the z-axis, the z-dependent ion diffusion constant is related to the ion bulk diffusion constant by the scaling factor

fz (β) =

D(z) 1 = . (β/C) Dbulk A + Be + De(β/E)

(16)

In this equation, β(z) denotes the ratio between the ion radius and the pore radius as a function of z, and A, B, C, D and E are a set of parameters. 26 Note that the scaling factor decreases as the ion radii increase. Here the Lennard-Jones radii were used as the ion radii and the pore radius was estimated using the HOLE algorithm. 52 Although it would be possible to fit ad hoc the scaling factor in the SV diffusion model for improving the agreement with some experimental results such as conductance measurements, we have rejected this possibility since it would involve arbitrary choices. It should be noted, however, that this factor might be fitted to generate ion diffusion profiles similar to those from all-atom MD simulations. Next the diffusion coefficients of the explicit atoms need to be defined. To this end, we

16


Page 17 of 51

!"#&$ !"#& (!)*&+,-.

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: Self-diffusion constant vs channel axis z. Shown are the single-valued (SV, dashed line) and hydrodynamics (HD, solid line) based diffusion values for K+ (red) and Cl− (blue) ions. The HD values for flexible GCMC/BD simulations (see Section 5.6) are shown as dotted lines.

17



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

Page 18 of 51

remind that the Stokes-Einstein model for the diffusion constant in solution is given by

D=

kB T , 6πηrα

(17)

where rα is the radius of a spherical particle moving in a continuous fluid of viscosity η. The Stokes-Einstein equation is valid for a sphere moving through a homogeneous, continuous, incompressible fluid. On a molecular level the continuum property implies that rα ≫ rβ , where rβ is the radius of a solvent molecule at a low Reynolds number. It has been found that for molecules below a certain size the experimental diffusion coefficient is systematically higher than determined from Eq. 17 and that the discrepancy increases as the molecules get smaller. 53 More generally, the conventional Stokes-Einstein formula can be written as 54

D=C

kB T , ηrα⋆

(18)

where C is a numerical constant and rα⋆ is an effective hydrodynamic radius. The interactions between diffusing particles and the environment serves to couple the solute motion to the solvent flow and eventually leads to the viscous dissipation of momentum and energy. The effective hydrodynamic radius should be regarded as a measure of the strength of this interaction. In the improved GCMC/BD algorithm, Eq. 18 has been implemented for estimating the diffusion coefficients of explicit atoms. The LJ radii were used as effective hydrodynamic radii. Moreover, C depends on the pore radius and is an input parameter which should be specified by the user. Since the pre-defined values for the effective hydrodynamic radii might not be adequate in some cases, BROMOCEA also allows the user to optionally specify values for the diffusion constants of selected atom types.

5.2

GCMC/BD system setup

As an example pore we use the OmpC porin of Escherichia coli in our simulations. This pore is one of the most abundant pores in the outer membrane of Gram-negative entero18


Page 19 of 51




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 input files for the fields φsf , φrf and Ucore (original scheme) were generated by the web interface of the GCMC/BD Simulator within CHARMM-GUI. 51,64 It should be noted that the system setup of a trimeric porin would be quite complicated compared to monomer treatment. This complication is related to the generation of several ion-accessible high-dielectric and ion-inaccessible low dielectric regions for a suitable computation of the potentials φsf , φrf and Ucore . Because the φsf and φrf potentials are relatively smooth, second-order B-splines were used to calculate the electrostatic energy and forces in the simulations. Third-order B-splines were used for the original core potential Ucore due to its abrupt variations. The OmpC channel was placed in a 30 Å-thick implicit lipid membrane and at the top and bottom of the simulation box, 5 Å-thick buffer regions of constant concentration were located at roughly 20 Å away from the protein in the z direction. All simulations were performed in 300 mM KCl solution at 300 K. According to above values, the GCMC buffer regions are at a reasonable distance from the channel, i.e., between 3-4 Debye lengths, such that any particle creation or destruction is screened from the channel. 65 The implicit membrane was located between −15 Å and 15 Å. In the alternative core-repulsive potential scheme, the LJ parameters for the membrane dummy atoms were set to εd = 0.01 Kcal/mol and σd = 6.236291 Å yielding a radius of 3.5 Å. To gain computational efficiency, the forces associated with LJ interactions between ions and dummy atoms were damped monotonically to zero in the interval between ron = 5.0 and rof f = 7.0 Å using the atom-based force-switching method. Following previous studies, 25,28 the dielectric constants for water, protein and lipid bilayer were set to 80, 2 and 2, respectively. The grid spacing for look-up potential maps was set to 0.5 Å. The ion parameters, i.e., LJ parameters, bulk diffusion constants, excess chemical potentials and coefficients for the water-mediated short-range potentials, were taken from previous studies. 21,25 For describing the bonded and nonbonded explicit atom parameters, the CHARMM36 force field 42,66 was used. This procedure represent a significantly different approach as compared with a previous study 36 in which the force constants were modified as a consequence of simulating the motion of only a few explicit protein atoms attached to

20


Page 20 of 51

Page 21 of 51

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

5.3

GCMC/BD simulation schemes

The new framework allows for different simulation schemes as described in Table 1. In the original GCMC/BD simulations, the OmpC channel is described as a low dielectric medium and all its atoms are implicitly considered. The improved GCMC/BD framework allows for combining two complementary a priori schemes in the description of ion permeation through a fixed channel configuration. In the restrained GCMC/BD simulations, some residues lining the pore with fixed positions are explicitly considered while the rest of the residues are implicitly treated and their effects are introduced via electrostatic, reaction field and corerepulsive steric potentials. Since the explicit atoms are maintained at the fixed positions, only the dynamics of the mobile ions is governed by the Langevin equation. However, the PMF includes steric and electrostatic contributions from the implicit OmpC atoms as well as shortrange and electrostatic interactions between the explicit OmpC atoms and the ions. The framework allows for incorporating thermal fluctuations in key residues lining the pore. In the fully-flexible GCMC/BD simulations, the dynamics of mobile particles is also described by the Langevin equation, while the multi-particle PMF includes steric and electrostatic contributions from implicit OmpC atoms as well as short-range and electrostatic interactions between all Brownian particle pairs. In addition to the GCMC/BD simulations, all-atom applied-field MD simulations were run for the OmpC channel under the same conditions as those for the GCMC/BD simulations (see details in the Supplementary Information). These simulations are used as reference and allow us to compare with the results obtained from GCMC/BD simulations. For each of the different GCMC/BD simulation schemes, the improved GCMC/BD algorithm allows for using the SV and HD ion diffusion models as well as the EM and SRK BD propagators. For the sake of simplicity, some abbreviations for the possible combinations of diffusion models and BD propagators are being defined in Table 2 which will be used in the 21



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 1: Different simulation schemes. Simulation scheme original GCMC/BD restrained GCMC/BD flexible GCMC/BD applied-field MD

Explicit atoms none restrained flexible flexible

Remainder of the pore implicit implicit implicit none

remainder of the paper. In the GCMC/BD simulations for the different schemes described in the next subsections, the SV and HD ion diffusion models were used together with the EM and SRK BD propagation algorithms. In each GCMC/BD cycle, one step of GCMC was performed for each step of BD. The atom-based force-switching method was employed to calculate the nonbonded forces in each BD step. Ten indepedent runs were carried to obtain statistically well-converged results. For the transport properties, i.e., conductivity estimates and the study described in Section 5.4.1, a transmembrane potential of Vmp = 0.1 V was applied. We note that throughout this study, the −z side of the membrane has higher electric potential than the +z side, such that the resulting electric field in the membrane region is along the +z direction, i.e., the direction from the periplasmic to the extracellular region. Furthermore, equilibrium simulations for estimating energetic properties were run under the same conditions as the analogous non-equilibrium simulations for estimating transport properties but with vanishing transmembrane potential, i.e., Vmp = 0. On the other hand, it should be stressed that the distance-dependent dielectric constant was only activated for ion-explicit atom and explicit atom-explicit atom interactions since the water-mediated short-range potential was adjusted to properly incorporate the short-range ion-ion behavior. 25 Water-mediated short-range potentials for ion-explicit atom and explicit atom-explicit atom interactions were not considered in the numerical simulations though available in the BROMOCEA code. In the case of ions, these ad hoc interactions were introduced to take into account hydration effects which have a short-range nature. For the case of explicit atoms, the procedure to obtain the coefficients for water-mediated short-range potentials seems hardly generalizable since explicit atoms are mutually interconnected generating complex structures

22


Page 22 of 51

Page 23 of 51

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


(see further discussion in Section 6). Table 2: Abbreviations defined for the possible combinations of diffusion models and BD propagators in a given GCMC/BD simulation scheme. Abbreviation Ion diffusion model BD propagator SV_EM SV EM SV_SRK SV SRK HD_EM HD EM HD_SRK HD SRK

5.4 5.4.1

Original GCMC/BD simulations Long jump exceptions

In several papers on BD simulations, the so-called long jump exceptions (LJE) for ions have been reported. 67–70 Briefly, during a BD time step ions can get very close to each other, due to the random forces. In such a case, the short-range repulsive force evaluated at the next time step becomes so large that its application during the subsequent time step leads to unphysically long jumps of individual ions. Several ad-hoc fixes can be introduced to deal with these unphysical ion configurations in the BROMOCEA code. A maximum displacement can be imposed to the displacements due to the systematic forces, i.e., forces obtained from the spatial gradient of the PMF. A collection of increasing repulsive layers inside ion-inaccessible low dielectric regions can be applied in combination with the original core-repulsive potential scheme so that these layers generate a repulsion gradient to aid ions in exiting the protein and membrane regions. 31 As discussed in Section 3, an alternative corerepulsive potential scheme is also available. It should be mentioned that a formal solution of the LJE is accomplished only when the hydrodynamic interactions among particles are properly incorporated into BD. 71 To study in detail the LJE, a set of simulations were performed using the original corerepulsive potential scheme in combination with a collection of repulsive layers as described above. Since spurious LJE may be caused by the application of the GCMC algorithm in the 23



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

Page 24 of 51

buffer regions, these regions are not considered in the LJE computations. The ron and rof f parameters for the switching function were set to 15 and 20 Å, respectively. Additionally, some tests were run without the atom-based force-switching method. Although large values were used for ron and rof f , the CPU time was increased by a factor of between 3 and 4 when atom-based switching was not applied (data not shown). The original GCMC/BD simulations were generated with a time step of ∆t = 10, 50, and 100 fs. Furthermore, the systems were equilibrated for 100 ns. The simulation times of the actual production runs were 0.1, 0.5 and 1 µs for time steps of 10, 50, and 100 fs, respectively. Although different simulation times were used for each time step, the same number of data points was used for obtaining the LJE since the total number of GCMC/BD cycles was the same in all cases. By taking ∆t small but non-zero, the overdamped Langevin equation models the process described by the Smoluchowski equation to an accuracy of order ∆t. To first order in ∆t, the solution of the Smoluchowski equation in the absence of hydrodynamics interactions is uniquely defined by the moments

∆ri (∆t)i =

Di Fi + ∇i Di ∆t kB T

(19)

and (20)

h∆ri2 (∆t)i = 6Di ∆t. Thereby, an estimation of the root-mean square displacement (RMSD) is given by

√

6Di ∆t.

The threshold displacement between consecutive steps defining a LJE, rLJE , depends on the time step and the ion diffusion coefficients. To consider this dependence, the RMSD from bulk ion diffusion was used as a guide to establish correlated rLJE at different time steps (see Table 3). Note that a rLJE = 1.0 Å was chosen arbitrarily to ∆t = 100 fs and the rLJE for the rest of time steps were assigned in correlation with this value. A maximum displacement due to the systematic forces of 10 Å was used to avoid unusual but problematic very long jumps while its effect was properly accounted for in the LJE estimations. The results for 24


Page 25 of 51

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 3: RMSD and rLJE for the different time steps ∆t. ∆t [fs] RMSD [Å] rLJE [Å] 100 0.3 1.0 50 0.2 0.7 10 0.1 0.3 Table 4: LJE for the original GCMC/BD simulations. The time step is denoted by ∆t. √ Uncertainties were calculated using σ/ n where σ denotes the standard deviation and n = 10 the number of simulations. ∆t [fs]

LJE [%] +

SV_EM SV_EM SV_EM SV_SRK SV_SRK SV_SRK HD_EM HD_EM HD_EM HD_SRK HD_SRK HD_SRK

10 50 100 10 50 100 10 50 100 10 50 100

K 0.003±0.000 0.033±0.005 0.386±0.065 0.003±0.000 0.007±0.000 0.028±0.000 0.003±0.000 0.027±0.000 0.085±0.010 0.003±0.000 0.007±0.000 0.027±0.000

Cl− 0.005±0.000 0.035±0.000 0.069±0.003 0.005±0.000 0.010±0.000 0.035±0.000 0.005±0.000 0.033±0.000 0.054±0.000 0.005±0.000 0.010±0.000 0.034±0.000

LJE are shown in Table 4. At a time step of 10 fs, the LJE percentage exhibits the same low-limit values, which are negligible, for the different diffusion models and BD propagators. At time steps of 50 and 100 fs, the low-limit values, i.e., LJE ∼ 0.01% and 0.03%, were achieved for simulations using the SRK propagation algorithm. The LJE percentage for ions is in all cases small and below 1% gradually increasing with increasing time step. To the best of our knowledge, this is the first study that investigates the LJE in the GCMC/BD framework. The crowding effects can become important for ion motion in confined spaces, i.e., in the selectivity filter of an ion channel. It can be argued that short time steps may be necessary to deal properly with the LJE. 70 However, we have obtained small LJE percentages in the original GCMC/BD simulations even for long time steps, i.e., 100 fs. This seems to be a consequence of including space-dependent ion diffusion terms in the BD algorithms implemented in the BROMOCEA code since self-diffusion constants 25



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

are significantly reduced in confined environments. Further, LJE could be more critical for studies of permeation through narrow pores, where the application of continuum models is questionable anyway. 5.4.2

Energetic and transport properties

To obtain converged conductance estimate, a set of simulations were run using the original and the alternative scheme to treat core-repulsive potentials. We used the same system setup with the production time of 1 µs. The maximum displacements due to the systematic forces were set to 1, 0.7 and 0.3 Å for simulations with time steps of 100, 50 and 10 fs, respectively. Thus, the total production time for several independent runs was 10 µs for each simulation type with a given time step, ion diffusion model, BD propagator and core-repulsive scheme. The results for the conductance and the average number of ions in the channel are shown in Table 5 (original core-repulsive potential scheme) and Table 6 (alternative core-repulsive potential scheme). The channel conductance was determined by counting the number of ions that cross the x-y plane at z = 0. The EM and SRK BD propagators combined with the alternative core-repulsive scheme require around 2 and 3 times more CPU time than the EM and SRK BD propagators combined with the original scheme. As expected, the CPU time increases linearly as the time step increases. The conductance estimates from simulations using different time steps for a given diffusion model, BD propagator and repulsive scheme, do agree remarkably well. Since only nonbonded ion-ion interactions are explicitly considered, the BD algorithms are stable even at long time steps. The conductance estimation is very sensitive to the applied ion diffusion model but the conductance difference between SV and HD diffusion models is roughly constant (i.e., around 0.3 nS) irrespective of the BD propagator and of the core-repulsive scheme. Further, the conductance estimation is also sensitive to the applied BD propagator so that the conductance difference between EM and SRK BD propagators is roughly constant (i.e., around 0.1 nS) irrespective of the ion diffusion model and of the core-repulsive scheme. The alternative core-repulsive scheme

26


Page 26 of 51

Page 27 of 51

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


slightly increases the number of ions in the channel but the conductance estimation is quite insensitive to this modification. An experimental conductance estimation from a linear interpolation from values at different temperatures was found to be 1.24 nS 10 which means roughly 0.41 nS for each monomer. In view of this experimental value, the better estimations for the conductance are achieved when the HD diffusion model and the SRK algorithm are used. In principle, this seems reasonable since the HD model represents the diffusion profile inside the pore more accurately by including the variation of pore size effect and the SRK is a second-order algorithm. In addition to physical measurable properties, reference simulations provide a complementary insights in the validation of the method. When one compares the average number of ions inside the channel with those from applied-field MD simulations (see Fig. S6), it is evident that the GCMC/BD simulations using rigid and implicit OmpC pore significantly underestimate the number of ions inside the channel although the cation preference of the pore is correctly reproduced. Indeed, most of the values for the average ratio NK + /NCl− are significantly higher than those reported for applied-field MD simulations (see Table 7). It is noteworthy that a few values seem to be correctly predicted. Nonetheless, this is an accidental consequence of the wide variance in the NK + /NCl− values, mainly due to the very low number of chloride anions inside the channel. Ions exhibit a uniform three-dimensional density or concentration in the bulk reservoirs. Inside a channel, however, their motions perpendicular to the pore axis are highly restricted by the channel walls. Therefore, the thermodynamics of the ions in the channel can approximately be described by 1D PMF along the pore axis denoted as W1D . 25 The PMF profile allows to estimate the work needed to translocate an ion across the porin. It shows where the main transport barriers exit but also where the wells, which attract ions, are located. For multion systems, the one-dimensional profiles can be calculated as

W1D,α = −kB T ln

27

ρα (z) ρbulk α

,


(21)


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

Page 28 of 51

Table 5: Average number of ions in the channel N , average ratio NK + /NCl− and average conductance G for the original GCMC/BD simulations using the original core-repulsive potential scheme. ∆t [fs] SV_EM SV_EM SV_EM SV_SRK SV_SRK SV_SRK HD_EM HD_EM HD_EM HD_SRK HD_SRK HD_SRK

10 50 100 10 50 100 10 50 100 10 50 100

NK + /NCl−

N +

K 3.01±0.01 3.01±0.01 3.01±0.01 2.77±0.01 2.81±0.01 2.82±0.01 3.06±0.01 3.06±0.01 3.09±0.01 2.68±0.01 2.70±0.01 2.74±0.01

G [nS]

−

Cl 0.14±0.01 23.08±1.35 0.92±0.01 0.11±0.00 28.01±1.12 0.92±0.01 0.10±0.01 31.70±1.33 0.91±0.01 0.07±0.00 38.33±2.14 0.84±0.01 0.07±0.00 39.59±2.12 0.85±0.01 0.07±0.01 44.61±4.75 0.86±0.01 0.14±0.01 22.82±1.36 0.57±0.01 0.12±0.01 25.47±1.45 0.57±0.01 0.12±0.01 27.07±2.75 0.59±0.01 0.04±0.00 64.94±5.16 0.51±0.01 0.04±0.00 60.00±6.74 0.52±0.00 0.04±0.01 74.66±11.47 0.53±0.01

where ρα (z) denotes the density of ion type α at position z and ρbulk its bulk density. The α obtained energy profiles reveal the change of the free energy of the ions in dependence of the z position along the channel axis. The free energy in the bulk is used as reference energy. The 1D PMFs along the pore axis were extracted from simulations using equilibrium conditions. The results are displayed in Figs. 4 and 5. The volume effect causes the barriers in all PMF profiles near the pore entrances, i.e., z ∼ ±15 Å. A cation selectivity can be inferred from these profiles. Since the protein-ion interactions are the main determinant for the ion distribution inside the channel, slightly lower average values of the ions in the channel than those reported from the conductance study, but in excellent correlation, are obtained in the equilibrium simulations. As observed, an excellent convergence is achieved for the K+ PMFs irrespective of the ion diffusion model, BD propagator and time step. Small discrepancies are seen for the Cl− PMFs which are more pronounced in the vicinity of the constriction zone. This finding is due to the larger statistical error in the estimation of the Cl− PMF profiles as a consequence of the extremely low number of Cl− ions inside the channel. On the other hand, similar K+ PMFs are obtained from different core-repulsive schemes where only

28


Page 29 of 51

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 6: Same as in Table 5 but using the alternative core-repulsive potential scheme. ∆t [fs] SV_EM SV_EM SV_EM SV_SRK SV_SRK SV_SRK HD_EM HD_EM HD_EM HD_SRK HD_SRK HD_SRK

10 50 100 10 50 100 10 50 100 10 50 100

N K+ 3.33±0.01 3.31±0.01 3.28±0.01 3.06±0.01 3.08±0.01 3.04±0.01 3.47±0.02 3.42±0.02 3.44±0.02 2.93±0.01 2.94±0.01 2.99±0.01

Cl− 0.34±0.02 0.30±0.01 0.28±0.01 0.21±0.01 0.22±0.01 0.16±0.01 0.48±0.03 0.41±0.02 0.41±0.02 0.16±0.01 0.14±0.01 0.04±0.01

NK + /NCl−

G [nS]

9.86±0.40 11.29±0.46 12.05±0.51 14.86±0.65 13.95±0.61 19.96±1.64 7.56±0.43 8.59±0.49 8.59±0.40 19.68±1.83 22.29±2.01 20.19±1.76

0.94±0.01 0.94±0.01 0.94±0.01 0.85±0.01 0.87±0.01 0.86±0.00 0.62±0.01 0.61±0.01 0.62±0.01 0.53±0.00 0.52±0.01 0.54±0.01

small differences are observed in the vicinity and inside of the constriction zone. However, important differences are observed between the Cl− PMFs from the difference core-repulsive schemes for z between 0 Å and 15 Å. This finding proves that (i) ion profiles are fairly sensitive to the choice of the core-repulsive scheme within the framework of a rigid and implicit pore model and that (ii) similar conductance estimates can be obtained from different rigid and implicit pore models which exhibit significant variations in the thermodynamic properties. Our simulations show that the conductance is predominantly determined by the ion diffusion inside the narrow transmembrane region while the observed PMF differences between different repulsive schemes are defined by forces acting outside of this region. In a next step, we explored the effect of different cutoff values in the atom-based force switching method on the energetic, i.e., PMF, and permeation properties. The same simulation scheme was employed with a time step of ∆t = 10 fs. The switching distances ron and rof f were set to 8 Å and 12 Å, respectively. It has been reported that values of rof f ≥ 12 Å were able to the reproduce no-cutoff results for MD simulations of carboxymyoglobin in water where a 4 Å switching region was used. 45 Taking into account that a (water) dielectric medium is employed in the GCMC/BD simulations and that the electrostatic polarization of the implicit salt in the outer region is included in the reaction field contribution, these cutoff 29



6"#+7#&'!()*# +,-#./'.0102'/

"(

"(

"'

"'

"&

"&

)!*"+,-./012/34

)!*"+,-./012/34

0"#+7#&'!()*# 45#./'.0102'/

"% "$ "! "#

"% "$ "! "#

ï! ï'# ï&# ï%# ï$# ï!#

ï! ï'# ï&# ï%# ï$# ï!#

"# "!# "$# "%# "&# "'# 5"+64

3"#$%#&'!()*# 45#./'.0102'/ "(

"'

"'

"&

"&

"% "$ "! "# ï! ï'# ï&# ï%# ï$# ï!#

"# "!# "$# "%# "&# "'# 5"+64

!"#$%#&'!()*# +,-#./'.0102'/

"(

)!*"+,-./012/34

)!*"+,-./012/34

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

Page 30 of 51

"% "$ "! "#

"#

ï! ï'# ï&# ï%# ï$# ï!#

"!# "$# "%# "&# "'#

5"+64

"#

"!# "$# "%# "&# "'#

5"+64

Figure 4: 1D multi-ion average PMF extracted from 10 independent original GCMC/BD simulations under equilibrium conditions (Vmp = 0) using the original core-repulsive scheme. The results for the K + ions are represented in red while for the Cl− ions in blue. The time steps ∆t = 100, 50, 10 fs are represented by solid, dashed and dotted lines, respectively. Notice that the overlapping among lines for different time steps of K+ PMFs is a consequence of the excellent convergence achieved in these cases. values are also very reasonable for GCMC/BD simulations. The results for the conductance and the average number of ions in the channel are shown in Table S1 while the PMF profiles are displayed in Figs. S1 and S2. As observed, a very good convergence is achieved in modelling conductance properties, i.e., the difference between conductances obtained with different cutoff schemes is lower than 0.1 nS in all studied cases (see Tables 5 and 6). It should be noted that the minor discrepancies in the Cl− PMFs around the constriction zone can be attributed to the same factors as ones discussed above. The statistical uncertainity in the PMF calculation for Cl− is a consequence of the low amount of anions inside the 30


Page 31 of 51

channel leading to a relatively poor statistics. Interestingly, it has been observed that the computational effort is quite similar to that needed using larger ron and rof f parameters, i.e., ron = 15 Å and rof f = 20 Å. However, a significant saving in CPU time is expected in those cases where higher ion concentrations and/or explicit atoms are chosen. For this reason, the values ron = 8 Å and rof f = 12 Å were used in the simulations described in Sections 5.5 and 5.6. Tests (data not shown) using restrained and flexible GCMC/BD simulations and the above sets of values for ron and rof f showed no significant differences in the obtained results. 2"#$%#&'()*+# $34#./'.!0!1'/

"(

"(

"'

"'

"&

"&

)!*"+,-./012/34

)!*"+,-./012/34

!"#$%#&'()*+# ,-#./'.!0!1'/

"% "$ "! "#

"(

"% "$ "! "#

ï! ï'# ï&# ï%# ï$# ï!#

ï! ï'# ï&# ï%# ï$# ï!#

"# "!# "$# "%# "&# "'# 5"+64

5"#67#&'()*+# ,-#./'.!0!1'/

"# "!# "$# "%# "&# "'# 5"+64

("#67#&'()*+# $34#./'.!0!1'/ "(

"' "'

"& )!*"+,-./012/34

)!*"+,-./012/34

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


"% "$ "! "# ï! ï'# ï&# ï%# ï$# ï!#

"& "% "$ "! "#

"#

ï! ï'# ï&# ï%# ï$# ï!#

"!# "$# "%# "&# "'#

5"+64

"# "!# "$# "%# "&# "'# 5"+64

Figure 5: Same as in Fig. 4 but using the alternative core-repulsive potential scheme.

31



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

5.5

Restrained GCMC/BD simulations

Mutations of residues on the L3 loop have been shown to significantly affect the pore selectivity. 72 The L3 loop is an integral part of the constriction zone and takes part in the conformational changes involved in channel opening and closing. In the restrained GCMC/BD simulations, several residues on the L3 loop and other residues located on the opposite barrel wall were explicitly considered. Specifically, the set of explicit atoms was composed of those atoms belonging to the L3 loop residues 103 to 115 which include the important negatively charged residues D105, E109 and D113 as well as those atoms belonging to the positively charged residues K16, R37, R74 and R124 on the opposite side of the wall. Thereby, the total number of explicit atoms was set to 289, including all hydrogen atoms. The criteria for selecting the set of explicit atoms was a compromise between computational efficiency and a proper description of the interactions in the constriction zone. To study the conductance, a set of simulations was run using two different core-repulsive schemes. The restrained GCMC/BD simulations were generated with a time step of 10 fs. The equilibration time was 20 ns, the production time 0.1 µs and a maximum displacement due to the systematic forces of 0.3 Å was used. Thus, the total production time for several independent runs was 1 µs for each simulation type with a given ion diffusion model, BD propagator and core-repulsive scheme. The results for the conductance, the average number of ions in the channel and the average ratio NK + /NCl− are shown in Table 7 (see also Table S2). To facilitate comparison, the average number of ions in the channel and the average ratio NK + /NCl− for applied-field MD simulations are also included in this table. Moreover, an estimate for the conductance value from MD simulations is included 10 which agrees nicely with the experimental value. The average GCMD/BD conductance values and their uncertainties are greater than the corresponding values for the original GCMC/BD simulations (see Tables 7 and S1). It can be seen that the average conductance increment is around 0.4 nS and 0.2 nS for the SV and HD ion diffusion models, respectively. Since the resolution of the conductance estimation is defined by the timescale of the simulations, 32


Page 32 of 51

Page 33 of 51

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 7: Average number of ions in the channel N , average ratio NK + /NCl− and average conductance G for original, restrained and flexible GCMC/BD simulations using the HD ion diffusion model and the SRK BD propagator. The MD results together with their standard deviations (in parentheses) are listed for comparison and detailed in the Supplementary Information concerning the ratio NK + /NCl− and in Ref. 10 concerning the conductance. SIMULATION original GCMC/BD original GCMC/BD restrained GCMC/BD restrained GCMC/BD flexible GCMC/BD flexible GCMC/BD applied-field MD

REPULSIVE ∆t [fs] SCHEME original 10 alternative 10 original 10 alternative 10 alternative 10 alternative 5

NK + /NCl−

N K+ 2.80±0.01 3.07±0.01 6.88±0.10 7.21±0.15 6.59±0.06 6.79±0.09 5.95 (1.64)

Cl− 0.05±0.00 0.17±0.02 0.46±0.12 0.76±0.17 0.79±0.07 0.95±0.10 1.87 (1.06)

63.71±6.50 0.55±0.01 20.01±2.03 0.56±0.01 14.01±2.73 0.77±0.03 11.03±3.42 0.76±0.02 9.08±0.83 0.63±0.04 8.60±1.51 0.71±0.07 5.42 (3.95) 0.41

it seems that a less accurate conductance estimate compared to the experimental value is a consequence of the shorter simulation time. As in previous simulations, the more accurate agreement with the experimental conductance is achieved when the HD diffusion model and the SRK algorithm are employed. On the other hand, the average number of ions in the channel is significantly greater than those for the original GCMC/BD simulations and is similar to that one obtained from applied-field MD simulations (see Fig. S6). The original GCMC/BD scheme generally yields average ratio NK + /NCl− which are higher than those obtained from applied-field MD simulations. The restrained GCMC/BD scheme, however, provides NK + /NCl− values closer to those estimated from the MD simulations as expected from an improved estimation of the average number of ions inside the channel. This finding is very remarkable in view of the important differences between the MD and restrained GCMC/BD simulations which include the facts that (i) the solvent is replaced by a highdielectric medium, that (ii) the flexible protein structure is replaced by a frozen configuration where only some atoms are explicitly considered and that (iii) the lipid bilayer is replaced by a slab formed either by a homogeneous dielectric material (original repulsive scheme) or by a layer of dummy atoms (alternative repulsive scheme). The 1D ion PMFs along the pore axis were extracted from simulations using equilibrium

33


G [nS]


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

conditions. The PMF profiles are displayed in Fig. 6 (see also Figs. S3 and S4). Important differences are observed between the profiles obtained from the original and restrained GCMC/BD simulations. For the potassium PMFs, one can observe an energy reduction in the barriers and wells as a first effect of including explicit atoms. Moreover, trapping of K+ due to less screened interactions between K+ and L3 loop (mainly negatively charge residues D105, E109 and D113) in the vicinity of and inside the constriction zone is also observed. The chloride PMFs from the restrained GCMC/BD simulations exhibit a deep well in the constriction zone. This may be explained in terms of two complementary effects. First, the interactions between Cl− and positively charge residues K16, R37, R74 and R124 are less screened as a consequence of considering the atoms belonging to these residues explicitly. Second, the trapped K+ cations also attract nearby Cl− anions. In an earlier study, 10 Pezeshki et al. extracted PMF profiles from equilibrium MD simulations. As compared to those from the present GCMC/BD simulations, these free energy profiles have the drawback of a poor statistic due to shorter simulations times. However, it is interesting to note that these profiles also predict a long-lived trapping of chloride anions in the constriction zone. In fact, the PMFs from MD and restrained GCMC/BD simulations reproduce successfully the weak cation selectivity of the OmpC porin. These results are consistent with those from average ratio NK + /NCl− shown in Table 7. In contrast, the original GCMC/BD simulations generate chloride PMFs in which this trapping effect is not observed and, thus, a high cation selectivity is erroneously predicted. The previous observations validate the treatment of interactions between explicit atoms and ions as implemented in the improved GCMC/BD scheme. In the original scheme, the channel protein is represented as an ion-inaccessible low dielectric region which contains a distribution of fixed partial atomic charges. As discussed earlier, 33 these approximations become problematic for narrow pores, e.g., due to uncertainties concerning the values of the dielectric constants. Defining dielectric interfaces in narrow pores is actually problematic by itself. As our results strongly suggest, the assumption of an implicit pore description

34


Page 34 of 51

Page 35 of 51




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

Because the inclusion of explicit atoms significantly increases the number of nonbonded interactions, even though these atoms are fixed, the CPU time increases significantly by a factor of about 16 compared to those times for the respective simulations where the OmpC pore is only considered implicitly.

5.6 5.6.1

Flexible GCMC/BD simulations Set of flexible atoms and adjustment of interdependent input parameters

The criteria for selecting a set of flexible channel atoms should consider two important aspects. First, the flexible atoms should obviously include those atoms which undergo significant fluctuations. Second, additional atoms should also be included which eliminate artefacts in the dynamics of nearby groups of flexible atoms and allow for a proper description of the interactions with those groups. In Fig. S8, the root-mean square fluctuations (RMSF) of the Cα atoms belonging to the residues forming the L3 loop are depicted for all-atom appliedfield MD simulations. In view of these RMSF values and the previous criteria, the list of flexible explicit atoms was composed of those atoms belonging to the residues 104 to 114 in the L3 loop. Furthermore, the residues K16, R37, R74 and R124 on the opposite side of the wall were also included in the list of flexible atoms since their inclusion allows for a proper description of the interactions with the flexible atoms belonging to the L3 loop. Finally, some fixed atoms were added to the list of explicit atoms such that the flexible ones can be fixed to the implicit rigid channel. Specifically, the list of fixed atoms was composed of the atoms belonging to the residues 103 and 115 in the L3 loop as well as a few additional atoms adjacent to each of the residues K16, R37, R74 and R124. It should be noted that the total list of explicit atoms was set to 313 atoms (244 mobile and 69 fixed atoms), including all hydrogen atoms and, thus, it was the same as in the previous subsection except for those additional atoms adjacent to the residues K16, R37, R74 and R124. In the flexible GCMC/BD simulations, the explicit atom self-diffusion constants, time step, maximum displacement for systematic forces and relative geometric tolerance of the 36


Page 36 of 51

Page 37 of 51

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


SHAKE algorithm for bond constraints are strongly interdependent. The self-diffusion constants for explicit atoms were estimated from Eq. 18 and two sets of diffusion coefficients were tested by fixing the C parameter to values of 1/(6π) and 1/(12π), respectively. The set generated by C = 1/(6π) corresponds to the conventional Stokes-Einstein model (see Eq. 17) for a spherical particle in an unbounded fluid. The decrease of self-diffusion constants inside channel confinements is mimicked by using a reduced value of C = 1/(12π). Because polar hydrogens involved in hydrogen bonds, i.e., the H and HC atom types in the CHARMM36 force field, show anomalously high values for their diffusion constants, these diffusion constants were manually adjusted to the values of backbone hydrogens, i.e., the HB1 atom type in the CHARMM36 force field. It should be noted that these anomalies are a consequence of the association between small LJ radii for H and HC atom types and the effective hydrodynamic ones. This fact reveals the limitations of the present model for estimating explicit atom diffusion constants and points out the necessity of developing a more realistic procedure for estimating these coefficients (see further discussion in the next section). For each set of explicit atom self-diffusion constants, the time step was selected according to a proper convergence of the SHAKE algorithm using a geometric tolerance of 10−4 as usual in MD simulations. Thereby, the time steps were set to ∆t = 10 fs and 5 fs for C = 1/(12π) and 1/(6π), respectively. A maximum displacement for systematic forces of 0.3 Å was selected in all cases. This value allows to discard unrealistic long displacements but still preserves stable structures. 5.6.2

Simulation procedure

Because the flexible GCMC/BD simulations are CPU time demanding (i.e., the CPU time increased by a factor of about 25 compared to the original GCMC/BD scheme for the example of the OmpC pore), one should optimize the simulation procedure to achieve a reasonable agreement between accuracy and computational effort. In equilibrium and non-equilibrium simulations, a procedure composed by three different stages was devised. In the first stage,

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

an initial configuration of the simulation system was obtained using a restrained GCMC/BD simulation. Thus, the same conditions as in Section 5.5 were employed with a time step of 100 fs, a maximum displacement due to the systematic forces of 1.0 Å and a simulation time of 50 ns. The next two stages were the equilibration and production runs. In these stages, the selected set of explicit atoms was allowed to fluctuate. All bond distances were constrained using the SHAKE algorithm. The equilibration time was 1 ns irrespective of the time step. The production time was 0.1 and 0.05 µs for the time steps of 10 and 5 fs, respectively. Thus, the total production time for several independent runs was 1 and 0.5 µs for a time step of 10 and 5 fs, respectively. It should be stressed that the HD ion diffusion profile was slightly modified since the pore radius was also modified (see Fig. 2). The pore radius was modified because the explicit atoms were excluded from the implicit pore region. Further, SV and HD diffusion models were not used for explicit atoms since their diffusion profiles remain unaltered in the constriction region, i.e., the region where they are located. The position of the L3 loop in OmpC may be determined by three types of interactions between the loop and the barrel wall: (i) a hydrogen bond network between the tip of the L3 loop and the adjacent barrel wall, (ii) salt bridges at the root of L3 and (iii) a transverse electrostatic field between negatively charged residues on the L3 loop and a cluster of positive charges on the opposite barrel wall. The OmpC porins are relatively insensitive to membrane potential differences and in experiment close only when this potential exceeds 0.25 V. 73 When forces due to the reaction field potential were applied to the mobile OmpC explicit atoms, the flexible GCMC/BD simulations exhibited a gating process where the L3 loop was bent into the opposite wall residues (data not shown). Since the transmembrane potential was 0 (equilibrium simulations) or 0.1 V (non-equilibrium simulations), this gating process seems to be an artifact of our simulations. Interestingly, the SHAKE algorithm, which is very sensitive to the self-diffusion and time step values, was able to maintain the bond distance constraints during the simulations. An important aspect of the reaction fields arising from the dielectric boundaries is its manifestation as a repulsive force acting on the particles as they diffuse in the

38


Page 38 of 51

Page 39 of 51

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


high dielectric bulk solution and approach a region with low dielectric constant (membrane or protein). Near a dielectric interface, the reaction field is equivalent to the interaction between the real particle charges and a fictitious charges located symmetrically on the opposite side of the interface. Due to the proximity of some mobile atoms to the interface created by the implicit pore atoms, detailed observations revealed that these atoms undergo unrealistic and abrupt reaction field forces which push the L3 loop inwards the constriction zone. This movement could be a consequence of the approximate estimation of the reaction field which would be subjected to important inaccuracies for mobile atoms. However, it is more likely that this problem is closely related to the unrealistic assumption that mobile explicit atoms and implicit ones belong to very different dielectric regions. In fact, they belong to the same protein structure. For these reasons, forces due to the reaction field were deactivated for mobile explicit atoms in the simulations described below. It should be noted, however, that the atom charges were still considered in the estimation of the reaction field for the ions. Similar as occurred with the reaction field, the flexible GCMC/BD simulations showed an unrealistic gating process when the original core-repulsive steric potential was applied. The original core-repulsive model only includes repulsive terms. In other words, it only takes into account in a rough way the repulsive part of the van der Waals interactions. However, the attractive dispersion terms are ignored although they can be important, together with the Coulomb interactions, in the stabilization of different conformations. In the previous subsections, we could observe important differences in PMF profiles between simulations using the two different repulsive schemes. Now, it can be observed that the original repulsive scheme introduces unrealistic displacements in the mobile atoms located close to the implicit OmpC surface. For these reasons, only the alternative core-repulsive steric potential was employed in the simulations described below.

39



"( "' )!*"+,-./012/34

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

"& "% "$ "! "# ï! ï'# ï&# ï%# ï$# ï!#

"#

"!# "$# "%# "&# "'#

5"+64

Figure 7: Multi-ion average PMF extracted from 10 independent flexible GCMC/BD simulations under equilibrium conditions (Vmp = 0) using the improved core-repulsive scheme, the HD ion diffusion model and the SRK BD propagator. The K+ ions are represented in red and Cl− ions in blue. GCMC/BD simulations using a time step of 10 fs are plotted in solid lines and as dotted lines using a time step of 5 fs. 5.6.3

Energetic and transport properties

The results for the conductance and the average number of ions in the channel are shown in Table 7 (see also Table S3). The results from simulations using a time step of 10 fs can be directly compared with those extracted from restrained GCMC/BD simulations since the same time step and production time were used in both cases. As can be seen, the conductance values are reduced to ∼ 0.35 nS and 0.13 nS for the SV and HD ion diffusion models, respectively, when flexibility is incorporated into the protein description. The improved conductance estimates are accompanied by enhanced predictions of the average numbers of ions inside the channel and of the average ratio NK + /NCl− as compared to the MD 40


Page 40 of 51

Page 41 of 51

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


results (see Table 7 and Fig. S6). Thus, the inclusion of flexibility in the constriction region of OmpC improves the calculation of transport properties compared to the rigid OmpC picture. On the other hand, the average number of ions and the conductance estimates from simulation using a time step of 5 fs and higher explicit atom diffusion constants are, within the statistical error, identical to those obtained from a time step of 10 fs and lower diffusion constants. However, slightly higher uncertainties are generally reported as a consequence of the lower resolution due to the shorter production time. Thus, an increase of the thermal atomic fluctuations seems not to significantly affect the transport properties. Moreover, the ion PMF profile is displayed in Fig. 7 (see also Fig. S5). A remarkable convergence is observed between the K+ PMFs from different explicit atom diffusion coefficients (and time steps). Small differences between the Cl− PMFs from different atom diffusion coefficients are related to the different resolution associated with each of simulations. Thus, an increase of the thermal atomic fluctuations seems not to significantly affect the mechanism of ion conduction. However, flexible GCMC/BD simulations yield PMF profiles which exhibit an improved cation selectivity compared to those from restrained GCMC/BD simulations (see Fig. 6). As pointed out previously, 33 it is expected that atomic fluctuations of the pore-lining atoms have a large effect on the ion permeation. Clearly, the flexibility of proteins allows for local displacements in response to the perturbing influence of an ion which in turn influences the ion permeation as in the case of OmpC.

6

Summary and future work

To bridge the gap between biological structure and function, it is highly desirable to develop computational techniques allowing direct accurate and quantitative modelling of experimental observables. The computer simulations of ion transport has a long track of successful studies of structure-function relationships and evolved into a powerful tool for studies of ion permeation mechanisms. At an all-atom level, one can introduce a voltage bias or an

41



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

electrochemical gradient in MD simulations and determine the corresponding electric current directly from the observed ion flux. 18,74–79 However, to obtain statistically significant results, one have to observe the passage of multiple ions through the channel during the simulation, which is computational demanding for channels. To reduce dimensionality of the problem at hand, BD simulations can be used as an attractive computational alternative for simulating process over long time scales since fast degrees of freedom are replaced by random Gaussian noise and electrostatic maps. The problem, however, is that some of the degrees of freedom represented by either static maps or random noise are essential for permeation process and thus have to be treated explicitly. The inclusion of explicit atoms into the GCMC/BD framework allows to relax the rigid channel assumption. This paper presents an improved GCMC/BD algorithm that properly incorporates position fluctuations of the protein atoms into the multi-particle PMFs and general Langevin equations. A new program package called BROMOCEA was introduced which includes force field energy potentials for explicit atoms and allows the simulation of permeation process across nanochannels with a realistic implementation of boundary conditions of concentration and transmembrane potential under equilibrium and non-equilibrium conditions. The improved GCMC/BD algorithm was applied under conditions, i.e., low ion concentration and transmembrane voltage, where all-atom applied-field MD fails to provide a reasonable conductance estimation unless extremly long simulations are performed (see examples in the Supplementary Information). The implemented features have been tested by studying the ion permeation across an OmpC pore and displayed excellent agreement to the experimental observations. It was found that the SRK propagator algorithm provides better conductance estimates as compared to the EM algorithm. The alternative core-repulsive potential developed in our work appears to be a valuable tool for stable simulations with mobile ions and dynamic treatment of protein residues. Moreover, the bond distance constraints and the atom-based force-switching method are able to reduce the computational

42


Page 42 of 51

Page 43 of 51

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


cost while preserving a proper representation of the system’s dynamics. The importance of including explicit atoms in the description of pores such as the OmpC porin has been demonstrated above. While the original GCMC/BD scheme fails to reproduce the weak cation selectivity of OmpC, the improved GCMC/BD scheme is able to overcome this deficiency and to provide a more realistic picture for ion permeation. Furthermore, we investigated the LJE in the GCMC/BD framework. A marginal percentage of LJE was obtained in all investigated cases as a consequence of using BD propagator algorithms which incorporate the variation of the self-diffusion constants inside the pore in a natural way. Several points for improving the BROMOCEA code in the future can be mentioned. First, it seems desirable to gain extra computational efficiency since the number of nonbonded interactions N (N − 1)/2 increases drastically with the number of particles N . Although cutoff distance techniques have been implemented to reduce the number of interacting pairs, the computational effort can still be prohibitive for systems including thousands of atoms. Two lines of action are currently being pursued in this direction. On the one hand, a multiple time step method 80 has been proposed within the GCMC/BD framework. On the other hand, a parallelization using OpenMP (Open Multi-Processing) seem a reasonable option. Second, a more realistic procedure for the estimation of atom diffusion constants is desirable. A procedure for extracting self-diffusion constants from all-atom MD simulations is currently in development. Third, it seems advisable to consider the molecular nature of the solvent in the interactions involving explicit atoms. Effective, solvent-mediated potentials for ion-explicit atom and explicit atom-explicit atom interactions might be generated from inverse Monte Carlo (MC) techniques. 81 Using radial distribution functions from fully atomistic MD simulations, the effective potentials are determined in an iterative process which requires MC simulations and solving a system of linear equations. These potentials have not a specific functional form requiring a grid approximation. Nonetheless, this method is not void of limitations in our particular case. Effective potentials should be generated for a large number of different atom types making an iterative process difficult. Furthermore,

43



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 inverse problem might be not well-defined in some cases. Last but not least, effective potentials should be generated for each different pore and, even more, for different regions in the same pore. As a fourth point, the inclusion of hydrodynamics effects among particles is desirable since such effects make BD more realistic and self-consistent. However, a drawback would be that the random variables for different particles are no longer statistically independent and generating these random variables could be very time-consuming. Fifth, one of the authors has developed a theoretical framework where the potential energy due to dipole polarization is included in the PMF using an atom dipole polarizability model. The inclusion of these effects improves the description of the electrostatic interactions in a medium with a high dielectric constant such as water, but it also makes the BD computations more time-consuming as well as it could make the parallelization difficult in some cases. Finally, the inclusion of explicit atoms in the GCMC/BD simulations allows direct and, more importantly, accurate simulations of the substrate translocation events. Preliminary simulations about the dynamics of the ciprofloxacin antibiotic inside of a OmpC porin provide results very encouraging. However, some enhanced sampling techniques should be incorporated into the GCMC/BD framework to deal with transitions between long-lived metastable states in the transport of molecules through nanopores.

Acknowledgement 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. The work in Calgary was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) (Discovery Grant RGPIN-315019 to SYN) and the Alberta Innovates Technology Futures (AITF) Strategic

44


Page 44 of 51

Page 45 of 51

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


Chair in BioMolecular Simulations (Centre for Molecular Simulation).

Supporting Information Available Implementation details concerning the electrostatics, the reaction field and the original corerepulsive steric potentials are given in the supporting information. The adoption of the SHAKE algorithm to the BD framework and the employed random number generator are described as well. Moreover, additional results from GCMC/BD and MD simulations are reported together with a movie from a flexible GCMC/BD simulation with 0.3 M KCl at a transmembrane voltage of 0.1 V. This information is available free of charge via the Internet at http://pubs.acs.org/.

References (1) Nature Nanotechnology Focus on Nanopores. http://www.nature.com/nnano/focus/nanopores, 2012. (2) Stavenger, R. A.; Winterhalter, M. Sci. Transl. Med 2014, 6, 228ed7. (3) Winterhalter, M.; Ceccarelli, M. Eur. J. Pharm. Biopharm. 2015, 95, 63–67. (4) Majd, S.; Yusko, E. C.; Billeh, Y. N.; Macrae, M. X.; Yang, J.; Mayer, M. Curr. Opin. Biotechnol. 2010, 21, 439–476. (5) Roux, B.; Allen, T.; Bernéche, S.; Im, W. Q. Rev. Biophys. 2004, 37, 15–103. (6) Maffeo, C.; Bhattacharya, S.; Yoo, J.; Wells, D.; Aksimentiev, A. Chem. Rev. 2012, 112, 6250–6284. (7) Modi, N.; Winterhalter, M.; Kleinekathöfer, U. Nanoscale 2012, 4, 6166–6180. (8) Im, W.; Roux, B. J. Mol. Biol. 2002, 319, 1177–1197. 45



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

(9) Pezeshki, S.; Chimerel, C.; Bessenov, A.; Winterhalter, M.; Kleinekathöfer, U. Biophys. J. 2009, 97, 1898–1906. (10) Biro, I.; Pezeshki, S.; Weingart, H.; Winterhalter, M.; Kleinekathöfer, U. Biophys. J. 2010, 98, 1830–1839. (11) Chen, M.; Khalid, S.; Sansom, M. S. P.; Bayley, H. Proc. Natl. Acad. Sci. USA 2008, 105, 6272–6277. (12) Khalid, S.; Bond, P. J.; Carpenter, T.; Sansom, M. S. Biochimica et Biophysica Acta (BBA)-Biomembranes 2008, 1778, 1871–1880. (13) Modi, N.; Bárcena-Uribarri, I.; Bains, M.; Benz, R.; Hancock, R. E.; Kleinekathöfer, U. ACS Chem. Biol. 2015, 10, 441–451. (14) van den Berg, B.; Bhamidimarri, P. S.; Prajapati, D. J.; Kleinekathöfer, U.; Winterhalter, M. Proc. Natl. Acad. Sci. USA 2015, 112, E2991–E2999. (15) Hajjar, E.; Bessonov, A.; Molitor, A.; Kumar, A.; Mahendran, K. R.; Winterhalter, M.; Pages, J.-M.; Ruggerone, P.; Ceccarelli, M. Biochemistry 2010, 49, 6928–6935. (16) Parkin, J.; Khalid, S. Biophys. J. 2014, 107, 1853–1861. (17) Acosta Gutierrez, S.; Scorciapino, M. A.; Bodrenko, I.; Ceccarelli, M. J. Phys. Chem. Lett. 2015, 6, 1807–1812. (18) Roux, B. Biophys. J. 2008, 95, 4205–4216. (19) Khalili Araghi, F.; Ziervogel, B.; Gumbart, J. C.; Roux, B. J. Gen. Physiol. 2013, 142, 465–475. (20) Ingolfsson, H. I.; Lopez, C. A.; Uusitalo, J. J.; de Jong, D. H.; Gopal, S. M.; Periole, X.; Marrink, S. J. WIREs Comput Mol Sci 2014, 4, 225–248. (21) Im, W.; Seefeld, S.; Roux, B. Biophys. J. 2000, 79, 788–801. 46


Page 46 of 51

Page 47 of 51

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


(22) Kuyucak, S.; Andersen, O. S.; Chung, S.-H. Rep. Prog. Phys. 2001, 64, 1427–1472. (23) Boda, E., D.and Csanyi; Gillespie, D.; Kristof, T. J. Phys. Chem. C 2014, 118, 700– 707. (24) Im, W.; Roux, B. J. Chem. Phys. 2001, 115, 4850. (25) Im, W.; Roux, B. J. Mol. Biol. 2002, 322, 851–869. (26) Noskov, S.; Im, W.; Roux, B. Biophys. J. 2004, 87, 2299–2309. (27) Egwolf, B.; Luo, Y.; Walters, D. E.; Roux, B. J. Phys. Chem. B 2010, 114, 2901. (28) Lee, K. I.; Rui, H.; Pastor, R. W.; Im, W. Biophys. J. 2011, 100, 611–619. (29) De Biase, P. M.; Solano, C. J. F.; Markosyan, S.; Czapla, L.; Noskov, S. Y. J. Chem. Theory Comput. 2012, 8, 2540–2551, PMID: 22798730. (30) De Biase, P. M.; Markosyan, S.; Noskov, S. J. Comput. Chem. 2014, 35, 711–721. (31) De Biase, P. M.; Markosyan, S.; Noskov, S. J. Comput. Chem. 2015, 36, 264–271. (32) Comer, J.; Aksimentiev, A. J. Phys. Chem. C 2012, 116, 3376–3393. (33) Allen, T. W.; Andersen, O. S.; Roux, B. J. Gen. Physiol. 2004, 124, 679. (34) Feenstra, K. A.; Hess, B.; Berendsen, J. C. J. Comp. Chem. 1999, 20, 786–798. (35) Noskov, S.; Bernéche, S.; Roux, B. Nature 2004, 431, 830–834. (36) Chung, S. H.; Corry, B. Biophys. J. 2007, 93, 44–53. (37) Brooks et al. , B. R. J. Comp. Chem. 2009, 30, 1545. (38) Roux, B. Biophys. J. 1997, 73, 2980–2989. (39) Murphy, T. J.; Aguirre, J. L. J. Chem. Phys. 1972, 57, 2098–2104. 47



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

(40) Tough, R. J.; Pusey, P. N.; Lekkerkerker, H. N. W.; van den Broeck, C. Mol. Phys. 1986, 59, 595–619. (41) Ermak, D. L.; McCammon, J. A. J. Chem. Phys. 1978, 69, 1352. (42) A. D. MacKerell, J.; Feig, M.; Brooks, C. L. J. Am. Chem. Soc. 2004, 126, 698. (43) A. D. MacKerell, J.; Feig, M.; Brooks, C. L. J. Comp. Chem. 2004, 25, 1400. (44) Mehler, E. L.; Guarnieri, F. Biophys. J. 1999, 75, 3. (45) Steinbach, P. J.; Brooks, B. R. J. Comp. Chem. 1994, 15, 667. (46) Heyes, D. M.; Brańka, A. C. Mol. Phys. 2000, 98, 1949. (47) Tironi, I. G.; Brunne, R. M.; van Gusteren, W. F. Chem. Phys. Lett. 1996, 250, 19–24. (48) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327–341. (49) van Gunsteren, W. F. Mol. Phys. 1980, 40, 1015–1019. (50) Knotts IV, T. A.; Rathore, N.; Schwartz, D. C.; de Pablo, J. J. J. Chem. Phys. 2007, 126, 084901. (51) Lee, K. I.; Jo, S.; Rui, H.; Egwolf, B.; Roux, B.; Pastor, R. W.; Im, W. J. Comput. Chem. 2012, 33, 331–339. (52) Smart, O. S.; Neduvelil, J. G.; Wang, X.; Wallace, B.; Sansom, M. S. J. Mol. Graph. 1996, 14, 354–360. (53) John T. Edward, J. Chem. Educ. 1970, 47, 261. (54) Zwanzig, R.; Harrison, A. K. J. Chem. Phys. 1985, 83, 5861. (55) Benz, R.; Schmid, A.; Hancock, R. E. W. J. Bacteriol. 1985, 162, 722–727. (56) Seltmann, G.; Holst, O. The Bacterial Cell Wall ; Springer: Berlin, 2002. 48


Page 48 of 51

Page 49 of 51

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


(57) Baslé, A.; Rummel, G.; Storici, P.; Rosenbusch, J. P.; Schirmer, T. J. Mol. Biol. 2006, 362, 933–942. (58) Suenaga, A.; Komeiji, Y.; Uebayasi, M.; Meguro, T.; Saito, M.; Yamato, I. Biosci. Rep. 1998, 18, 39–48. (59) Zachariae, U.; Koumanov, A.; Engelhardt, H.; Karshikoff, A. Protein Sci. 2002, 11, 1309–1319. (60) Zachariae, U.; Helms, V.; Engelhardt, H. Biophys. J. 2003, 85, 954–962. (61) Danelon, C.; Suenaga, A.; Winterhalter, M.; Yamato, I. Biophys. Chem. 2003, 104, 591–603. (62) Nestorovich, E. M.; Danelon, C.; Winterhalter, M.; Bezrukov, S. M. Proc. Natl. Acad. Sci. USA 2002, 99, 9789–9794. (63) Rostovtseva, T. K.; Nestorovich, E. M.; Bezrukov, S. M. Biophys. J. 2002, 82, 160169. (64) Jo, S.; Kim, T.; Iyer, V. G.; Im, W. J. Comput. Chem. 2008, 29, 1859–1865. (65) Corry, B.; Hoyles, M.; Allen, T. W.; Walker, M.; Kuyucak, S.; Chung, S.-H. Biophys. J. 2002, 82, 1975–1984. (66) Klauda, J. B.; Venable, R. M.; Freites, J. A.; O’Connor, J. W.; Tobias, D. J.; Mondragon Ramirez, C.; Vorobyov, I.; MacKerell Jr, A. D.; Pastor, R. W. J. Phys. Chem. B 2010, 114, 7830–7843. (67) Heyes, D. M.; Brańka, A. C. Mol. Phys. 1998, 94, 447–454. (68) Jardat, M.; Bernard, O.; Turq, P.; Kneller, G. R. J. Chem. Phys. 1999, 110, 7993. (69) Batôt, G.; Mériguet, V.; Louis, A. A.; Jardat, M. Phys. Rev. E 2013, 88, 043304.

49



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

(70) Berti, C.; Furini, S.; Gillespie, D.; Boda, D.; Eisenberg, R. S.; Sangiorgi, E.; Fiegna, C. J. Chem. Theory Comput. 2014, 10, 2911–2926. (71) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: New York, 1987. (72) Liu, N.; Delcour, A. H. Protein. Eng. 1998, 11, 797–802. (73) Mobasheri, H.; Lea, E. J. Eur. Biophys. J. 2002, 31, 389–399. (74) Aksimentiev, A.; Schulten, K. Biophys. J. 2005, 88, 3745–3761. (75) Khalili Araghi, F.; Tajkhorshid, E.; Schulten, K. Biophys. J. 2006, 91, L72–74. (76) Sotomayor, M.; Vasquez, V.; Perozo, E.; Schulten, K. Biophys. J. 2007, 92, 886–902. (77) Jensen, M.; Borhani, D.; Lindorff Larsen, K.; Maragakis, P.; Jogini, V.; Eastwood, M.; Dror, R.; Shaw, D. Proc. Natl. Acad. Sci. USA 2010, 107, 5833. (78) Kutzner, C.; Grubmüller, H.; de Groot, B. L.; Zachariae, U. Biophys. J. 2011, 101, 809–817. (79) Rui, H.; Lee, K. I.; Pastor, R. W.; Im, W. Biophys. J. 2011, 100, 602–610. (80) Streett, W. B.; Tildesley, D. J. Mol. Phys. 1978, 35, 639–648. (81) Lyubartsev, A. P.; Laaksonen, A. Phys. Rev. E 1995, 52, 3730–3737.

50


Page 50 of 51

Page 51 of 51

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

51


BROMOCEA Code: An Improved Grand Canonical Monte Carlo

Recommend Documents