Assessment of Common Simulation Protocols for Simulations of
Nanopores, Membrane Proteins, and Channels
Jirasak Wong-ekkabut*,†
and Mikko Karttunen*,‡
†
Department of Physics, Faculty of Science, Kasetsart University, 50 Phahon Yothin Road, Chatuchak, Bangkok 10900, Thailand
‡
Department of Chemistry, University of Waterloo, 200 University Avenue West, Waterloo, Ontario, Canada N2L 3G1
ABSTRACT: Molecular dynamics (MD) simulation has become a common technique to study biological systems. Transport of
small molecules through carbon nanotubes and membrane proteins has been an intensely studied topic, and MD simulations
have been able to provide valuable predictions, many of which have later been experimentally proven. Simulations of such
systems pose challenges, and unexpected problems in commonly used protocols and methods have been found in the past few
years. The two main reasons why some were not found before are that most of these newly discovered errors do not lead to
unstable simulations. Furthermore, some of them manifest themselves only after relatively long simulation times. We assessed the
reliability of the most common simulations protocols by MD and stochastic dynamics (SD) or Langevin dynamics, simulations of
an alpha hemolysin nanochannel embedded in a palmitoyloleoylphosphatidylcholine (POPC) lipid bilayer. Our ﬁndings are that
(a) reaction ﬁeld electrostatics should not be used in simulations of such systems, (b) local thermostats should be preferred over
global ones since the latter may lead to an unphysical temperature distribution, (c) neighbor lists should be updated at all time
steps, and (d) charge groups should be used with care and never in conjunction with reaction ﬁeld electrostatics.
■ INTRODUCTION
Molecular dynamics (MD) simulation is a powerful technique
to study the structures, dynamics, and interactions of biological
molecules.1−3
It has been successful in predicting physical
phenomena such as the ability of water molecules to penetrate
carbon nanotubes4,5
and the collective nature of lipid diﬀusion.6
Simulations can be quantitative and predictive: Both of the
above predictions and many others have been later veriﬁed by
experiments.7
Simulations can also tell how probes used in
experiments may disturb8,9
the system and hence help to
provide a more accurate interpretation of experiments.
Like all techniques, MD simulation has both advantages and
disadvantages. Algorithmic properties are one concern,10
but
over the past years force ﬁelds11−21
and methodological issues
such as how to include long-range electrostatic interac-
tions22−37
accurately and reliably and the long-time stability
of simulations have been the main concerns. Eﬀorts for more
quantitative and accurate modeling have led to improved
scrutiny of the methods, and the past couple of years have
provided surprises: Bonthuis et al.38,39
found that truncation of
the Lennard−Jones interactions may in some cases lead to
unexpected serious artifacts, and Wong-ekkabut et al.40
found
that some of the commonly used protocols may lead to
deceivingly beautiful but unphysical behavior, such as the
spontaneous ﬂow of water in systems containing narrow
channels such as nanotubes or nanopores. Similar conclusions
using diﬀerent software were reached by Zuo et al.41
In
addition, Ni et al.42
found that inappropriate treatment of
electrostatic interactions produces artiﬁcial repulsions between
charged residues in simulations containing DNA. This was
identiﬁed as the reason behind artiﬁcial melting of doublestranded
DNA. In this article, we use a membrane protein as
our test system and study how the method for computing
electrostatic interactions, reaction ﬁeld vs particle-mesh Ewald,
choice of thermostats, and using neighbor lists and charge
groups to speed up computations may change the physical
behavior of the system. At the end, we will provide some
recommendations for choosing a safe simulation protocol.
The most recent artifacts38,40,43
have been found in systems
that contain narrow channels through which molecules, such as
water, ions, or DNA, can move. These systems are
characterized by the presence of highly anisotropic areas of
higher and low density of molecules. Since such systems are of
great interest in both biology and nanotechnology, it is
important to assess how current methods can be used reliably
in modeling them. In this study, we ﬁrst tested how sensitive
the simulations of a membrane protein are to parameter choices
and then studied permeation of water molecules through the
channel. Our results are applicable to all systems that contain
narrow channels, for example, aquaporin,44−46
carbon nano-
tubes,4,47
and nanocontainers.48
As our test system, we use the bacterial toxin alpha
hemolysin (AHL) secreted by the human pathogen Staphylococus
aureus. AHL is heptameric, and its 33.2 kDa watersoluble
monomers bind to susceptible cell and self-assemble to
form a 232.4 kDa transmembrane pore. High-resolution X-ray
crystallography49
shows that the protein has a mushroom shape
and is about 10 nm high. The main structure consists of a cap
and beta-barrel stem, see Figure 1. The largest diameter inside
the pore is about 4.6 nm and is located in the cap. The
narrowest part is in the stem with a diameter of about 1.6 nm,50
Figure 1. The pore is nonselective and thus allows for transport
of water, ions, and small molecules.
The nonselective nature of transport exhibited by AHL is
appealing for applications in biotechnology, nanotechnology,
Received: February 15, 2012
Published: June 26, 2012
Article
pubs.acs.org/JCTC
© 2012 American Chemical Society 2905 dx.doi.org/10.1021/ct3001359 | J. Chem. Theory Comput. 2012, 8, 2905−2911
and medicine. Suggested applications include bombardment of
cancer cells with proteins that would damage cells’ outer
membranes, making them susceptible to chemotherapy,51
biosensors, and DNA sequencing.52,53
The wild-type AHL
channel shows a weak anion conductivity,54,55
and site-directed
mutations inside the channel,55,56
such as K147N or M113P,
can be made to tune the channel to be moderately cation
selective. AHL can hence be used as a good and physically
relevant test system.
■ MOLECULAR DYNAMICS SIMULATIONS
We performed classical MD simulations of AHL49
embedded
into a palmitoyloleoylphosphatidylcholine (POPC) lipid bilayer
to study the permeation of water into the protein pore (Figure
1). The structure was taken from the protein data bank (PDB
ID 7AHL).49,57
The Gromos 53A5 force ﬁeld19
was used. We
also performed simulations with Gromos53a619
and observed
no changes in the systems’ behaviors.58
AHL is protonated at neutral pH,59
resulting in a positively
(+7) charged structure. To maintain overall charge neutrality, 7
Cl−
counterions were added.60
A pre-equilibrated lipid bilayer
structure was used,61
and the double-bond parameters of
POPC lipids62
were corrected following Bachar et al.11,63
The
system consisted of one AHL, 401 POPC lipids, 7 Cl−
counterions, and 36 747 water molecules. For water, we used
the simple point charge model (SPC).64
The box size was 12.9
nm × 12.8 nm × 12.0 nm. After energy minimization using the
steepest decent algorithm to remove cavities and close contacts,
an equilibration MD run was performed over 100 ns. The
production runs for each of the independent systems were run
for 40−80 ns, and the initial structure was taken from the last
state of equilibration run. This same structure (and identical
setup) was used in all production simulations. Simulations were
carried out with the GROMACS package version 4.0.565
in the
NVT (constant particle number, volume, and temperature)
ensemble. The integration time step was set to 1 fs, and the
trajectories were saved every 2 ps for analysis. Periodic
boundary conditions were applied in all directions. Protein,
lipids, and water molecules were thermostatted separately at
323 K (ions were grouped with water molecules). Simulation
times and details are listed in Tables 1 and 2.
It has been shown that simulations of nanotube systems may
lead to an unphysical ﬂow.40,41
To eliminate possible errors and
see where they originate from, we tested diﬀerent algorithms
and procedures including the frequency of neighbor list
updates, long-range electrostatic interactions, thermostats, and
charge groups to calculate interaction cutoﬀ. Neighbor list
update and charge groups are speciﬁc to Gromacs, but the rest
apply to all simulation software; the concept of charge groups
in Gromacs are explained in detail in the section “Domain
Decomposition” in Hess et al.,65
and a detailed study of their
eﬀect on physical properties has been done by Wohlert and
Edholm.66
The lists of simulation parameters used in MD and
stochastic dynamics/Langevin dynamics (SD)67
are shown in
Tables 1 and 2, respectively. Visualizations were done using
Visual Molecular Dynamics (VMD) software.68
Figure 1. Snapshot from an MD simulation (a) from the side and (b) the same system as seen from above. Protein (radius gradient color ribbon) is
embedded in a lipid bilayer, and water molecules (red and white) hydrate the system with a density of about 1000 kg/m3
. Protein pore was occupied
by 1426 ± 38 water molecules. (Inset) Water molecules that are inside the pore.
Table 1. List of Simulation Parameters in the Molecular
Dynamics Simulationa
interaction
basis
temperature
coupling
electrostatics
calculation
frequency
of neighbor
list update
(nlist)
time
(ns)
system 1 charge
group
Berendsen reaction ﬁeld 10 50
system 2 charge
group
Berendsen reaction ﬁeld 1 40
system 3 charge
group
Berendsen PME 10 80
system 4 charge
group
Berendsen PME 1 80
system 5 charge
group
V-rescale PME 10 80
system 6 charge
group
V-rescale PME 1 80
system 7 atom Berendsen reaction ﬁeld 10 50
system 8 atom Berendsen reaction ﬁeld 1 15
a
The frequency of neighbor list updates (nlist; speciﬁc to Gromacs).
Long-range electrostatic interactions using reaction ﬁeld69
and particlemesh
Ewald (PME).70,71
Time refers to the length of the production
simulation after 100 ns equilibration. Thermostats: Berendsen weak
coupling76
and the Parrinello−Donadio−Bussi velocity rescale.81,82
Charge group and atomic group42,65
were used as a basis to compute
the interaction cut oﬀ. For charge group, default values were used.
Journal of Chemical Theory and Computation Article
dx.doi.org/10.1021/ct3001359 | J. Chem. Theory Comput. 2012, 8, 2905−29112906
Lennard−Jones interactions with a cutoﬀ radius of 1.0 nm
were applied in all simulations. For electrostatics interactions,
we compared the reaction ﬁeld69
and particle-mesh Ewald
(PME) methods.70,71
In simulations using the reaction ﬁeld
method, all electrostatic interactions were explicitly computed
within a cutoﬀ radius of 2.0 nm. Beyond the cutoﬀ a dielectric
constant of 62 was used.72
The cutoﬀ used here is very
conservative as very commonly either 1.2 and 1.8 nm is used.
Simulations with shorter cutoﬀs than the one used here
accentuate the observed artifacts. These issues are general to all
simulation software using reaction ﬁeld electrostatics.
In the simulations using the PME method, a direct space
cutoﬀ of 1.0 nm was applied. In the calculation of long-range
contribution in reciprocal space, we used cubic interpolation of
order four using a maximum spacing of 1.2 Å for the FFT grid
and the relative strength of the electrostatic interaction at the
cutoﬀ was 10−5
.
■ SIMULATION RESULTS AND DISCUSSIONS
1. Unphysical Unidirectional Flow. We ﬁrst examine the
possible appearance of unphysical ﬂow and its dependence on
simulation parameters as seen in an analogous carbon nanotube
system.40
Issues concerning the diﬀerent thermostats and
methods for computing the electrostatic interactions are
common to all simulation software, such as Gromacs,65
CHARMM,73
Amber,74
Espresso,75
and so on. As mentioned
above, charge groups and tunable neighbor list update
frequency are speciﬁc to Gromacs.65
In brief, a charge group
is a small set of nearby atoms that are grouped together.
Computation of electrostatic interactions is based on the
distance from the geometric center of this group: If a neighbor
that does not belong to a group is farther away from the
geometric center of the charge group than the cutoﬀ distance,
its direct electrostatic interactions are not computed with any of
the atoms belonging to the group even if the distance to some
of the members of the group is shorter than the cutoﬀ length.
This can provide a signiﬁcant speed up and is conceptually
simple, but charge groups have to be used very carefully in
order to avoid physical artifacts.40
As a default, Gromacs
neighbor lists are updated every 10 time steps (parameter nlist
= 10). Although this is often a safe choice, it may lead to
unphysical behavior in highly anisotropic systems since
interactions could be missing from the calculation when lists
were not suﬃciently updated. The problem is that wrong
choices do not lead to an unstable simulation: The algorithm
remains stable, but the results are physically wrong. The only
way to ensure the correctness of parameters is to perform tests
of signiﬁcant length. In our experience, the minimum length of
such tests is at least tens of nanoseconds.
Case 1: Eﬀect of Neighbor Lists and Reaction Field
Electrostatics.69
We will examine each system in detail. The
labels “system 1” and so on refer to the systems in Table 1.
• System 1: Neighbor list update was set to 10 (default
value in Gromacs). The Berendsen weak coupling
thermostat76
and reaction ﬁeld electrostatics69
were
applied. Default charge groups were used. We tracked
the number of water molecules ﬂowing from cis to trans
and vice versa. Figure 2 shows that this combination
leads to an artiﬁcial ﬂow. In our production simulation,
942 water molecules were transported from the cis to the
trans side and only 148 moved in the opposite direction.
• System 2: As in system 1, but instead of setting nlist =
10, we updated the neighbor lists at every time step (nlist
= 1). This is the typical setup in software other than
Gromacs. Somewhat surprisingly, the net ﬂux was found
to be even larger than for nlist = 10. We could not clearly
identify the reason for the increase of the net ﬂux, but the
main reason for it is in reaction ﬁeld electrostatics.
Case 2: Reaction Field vs Particle-Mesh Ewald. The next
possible error source is calculation of the long-range electrostatic
interactions. There are three main approaches to treating
electrostatic interactions: truncated electrostatics, Ewald
summation-based methods, and reaction ﬁeld which is
essentially truncation plus mean ﬁeld correction. A review of
diﬀerent methods is provided by Karttunen et al.24
First,
possible artiﬁcial periodicity is largely attenuated in the reaction
ﬁeld approach, but this has been debated in the past several
years.77
Second, the reaction ﬁeld method for electrostatics
oﬀers a reduction in the computer time required by assuming a
uniform dielectric constant of water beyond the short-range
cutoﬀ value. As parallelization of the FFT has become more
eﬃcient, PME has become very competitive against the
Table 2. List of Simulation Parameters Used in SD
Simulations67
(Langevin dynamics)a
temperature
coupling
inverse
friction
constant
(ps)
electrostatics
calculation
frequency of
neighbor list
update (nlist)
time
(ns)
system 9 Langevin 2 reaction ﬁeld 10 50
system
10
Langevin 0.1 reaction ﬁeld 10 50
system
11
Langevin 2 PME 10 50
system
12
Langevin 0.1 PME 10 50
system
13
Langevin 2 PME 1 50
system
14
Langevin 0.1 PME 1 50
a
In all simulations, the Langevin thermostat was used to control the
temperature with the inverse friction constants of 2 and 0.1 ps and
default values of charge group were used. Production simulations were
run for 50 ns, and the frequency of neighbor list updates (nlist) was
varied between 1 and 10 steps. Long-range electrostatic interactions
using reaction ﬁeld69
and particle-mesh Ewald (PME)70,71
were tested.
Figure 2. Net number of water molecules (n) successfully permeated
from cis to trans as a function of time in MD simulations. Simulation
parameters for each system are provided in Table 1.
Journal of Chemical Theory and Computation Article
dx.doi.org/10.1021/ct3001359 | J. Chem. Theory Comput. 2012, 8, 2905−29112907
reaction ﬁeld78
if a long cutoﬀ is used in reaction ﬁeld, PME
become computationally faster.78
The reaction ﬁeld approach
requires the dielectric constant as an input. In real systems the
dielectric constant is, however, not homogeneous and may lead
to errors. At the cutoﬀ there is also a discontinuity in the
electrostatic potential, and this may create an artiﬁcial force.79,80
Therefore, the particle-mesh Ewald (PME) technique70,71
should be used instead of the reaction ﬁeld method.24,27,78
• System 3: The simulation was run with the same
parameters as in system 1 (nlist = 10), except that PME
was used instead of the reaction ﬁeld method. The results
show that the net ﬂux is signiﬁcantly reduced compared
to the reaction ﬁeld simulations. The net ﬂux in the
production simulation was 423 water molecules transported
from the cis to the trans side.
• System 4: nlist = 1. As with system 2, this setup is typical
in software other than Gromacs. The net ﬂux nearly
disappeared (Figure 2): PME electrostatics is clearly
preferable to reaction ﬁeld.
Case 3: Eﬀect of Thermostat. In this case, we used the same
setups as in systems 3 and 4 but instead of the Berendsen weak
coupling thermostat76
we used the Parrinello−Donadio−Bussi
velocity rescale thermostat.81,82
• System 5: nlist = 10, PME electrostatics. The behavior
was very similar to system 3.
• System 6: nlist = 1, PME electrostatics. The net ﬂux
disappeared. Comparison with system 4 shows that the
thermostat clearly has an eﬀect and that the velocity
rescale thermostat is the preferred choice. Comparison of
systems 5 and 3 shows that errors due to the neighbor
list updates and thermostat are not easily separable into
distinct contributions.
Case 4: Charge Groups vs No Charge Groups. We
compared the eﬀect of removing charge groups using the
worst cases from the above in order to gain more insight into
the diﬀerent error sources. We used the systems from case 1,
i.e., reaction ﬁeld electrostatics and the Berendsen weak
coupling thermostat. In terms of computational eﬃciency,
simulations using atom-based truncation with nlist = 1 are
remarkably expensive, approximately an order of magnitude
slower than using nlist = 10 (with atom based, i.e., no charge
groups).
• System 7: nlist = 10, no charge groups. The observed
behavior was very similar to system 1; the net ﬂux
converged approximately to the same value, although the
ﬂux started later.
• System 8: nlist = 1, no charge groups. We did not run
this simulation further than about 15 ns, but there is a
clear improvement compared to system 2 (the same
system with default charge groups).
Case 5: Simulations Using Stochastic Dynamics. Finally,
we performed SD simulations (Figure 3) in order to ensure that
the velocities of particles in the system were generated in the
correct ensemble.
• Systems 9 and 10: Parameters similar to system 1 were
used (reaction ﬁeld electrostatics, nlist = 10). Inverse
friction coeﬃcients of 2 and 0.1 ps were applied to the
systems 9 and 10, respectively. The observed behavior
was very similar to the MD simulations of systems 1 and
2, i.e., a net ﬂux appeared. The ﬂux was smaller with
smaller inverse friction coeﬃcient; the dynamics of the
water molecules is damped by a high friction coeﬃcient
(smaller inverse friction).
• Systems 11 and 12: PME electrostatics was used instead
of reaction ﬁeld. The neighbor lists were updated every
10 time steps, and again, inverse friction coeﬃcients of 2
and 0.1 ps were used in systems 11 and 12, respectively.
Similar to the previous MD simulations, the net ﬂux was
signiﬁcantly reduced in comparison to systems 9 and 10,
and it almost disappeared when an inverse friction
coeﬃcient of 0.1 was used.
• Systems 13 and 14: The neighbor lists were updated
every time step. Inverse friction coeﬃcients of 2 and 0.1
ps were used in systems 13 and 14, respectively. The net
ﬂux disappeared in system 14. This is consistent with the
results from system 6, MD with the velocity rescale
thermostat. This is also similar to the results reported by
Zuo et al.41
using NAMD software.83
2. Pore Shape and Water in Conﬁned Geometry. As a
brief application we studied the density of water inside the
pore. On the basis of the previous section, we chose charge
group-based electrostatic interactions, V-rescale thermostat,
PME, and nlist = 1 (system 6) to ensure that the simulation
setup is physically correct. We used charge groups instead of
atoms based for speed since system 6 showed no artifacts. We
started by analyzing the water density within the pore. The
results are based on an 80 ns production simulation and shown
in Figure 4. The error bars are given by the standard deviation.
The proﬁle in Figure 4A is consistent with previous work.54,84
It
is also comparable with the open state of the mechanosensitive
channel of a large conductance protein (MSL) channel.85
The
AHL protein pore has a funnel-like shape with a wide cap (6.3
< z < 10.5 nm) and a narrow stem (1.5 < z < 6.3 nm). The
narrowest and widest radii inside the pore were measured to be
0.79 ± 0.08 and 2.07 ± 0.09 nm at z = 6.3 and 7.5 nm from the
simulation box bottom (or a distance of 4.8 and 6.0 nm from
the end of the protein stem), respectively. The protein pore was
occupied by a total of 1426 ± 38 water molecules. Water
density in the cap region is close to the bulk density of water,
but in the stem the density is signiﬁcantly lower. Moreover, in
the stem region there are two density minima: 548 ± 93 kg/m3
at z = 2.7 nm and 526 ± 86 kg/m3
at z = 5.7 nm. These are the
locations of the hydrophobic residues, leucine (Leu135,
Leu428, Leu721, Leu1014, Leu1307, Leu1600, and Leu1893)
Figure 3. Net number of water molecules (n) successfully permeated
from cis to trans as a function of time in SD simulations. Simulation
parameters for each system can be seen in Table 2.
Journal of Chemical Theory and Computation Article
dx.doi.org/10.1021/ct3001359 | J. Chem. Theory Comput. 2012, 8, 2905−29112908
and methionine (Met113, Met406, Met699, Met992, Met1285,
Met1578, and Met1871), respectively. These two minima were
also observed by Cozmuta et al.84
A decrease in water density
has also been shown in neutral and irregularly charged carbon
nanotubes.86−88
According to MD simulations,87,88
the carbon
nanotube channel was ﬁlled with a tubular-like water structure
at a central water core having near-randomly distributed
molecules due to thermal movement. It has also been found
that water density decreases as the pore diameter decreases.
The water density in the stem is equal to the density of water in
a 2 nm diameter carbon nanotube.87
We did not, however,
observe the reported tubular-like87,88
water structure inside this
protein channel.
■ CONCLUSIONS
We performed atomistic MD and SD simulations to investigate
error sources and deﬁne a safe protocol for simulations of
systems containing pores and channels. We used an AHL
protein embedded into a lipid bilayer as model and investigated
the permeation of water molecules through the protein channel.
We found that an artiﬁcial unidirectional ﬂow was mainly
related to neighbor list update, calculation of electrostatic
interactions, and choice of thermostat. The unidirectional ﬂow
was signiﬁcantly reduced when PME was used, and neighbor
lists were updated at every time step. When the Parrinello−
Donadio−Bussi V-rescale81,82
thermostat was used (with the
above protocol) instead of the Berendsen weak coupling
thermostat76
the unidirectional ﬂow disappeared. On the basis
of our ﬁndings, we recommend the following.
• Tests should be run over a suﬃciently long time. Some
of the errors manifest themselves only after tens of
nanoseconds: The results in Figure 1 were preceded by
100 ns of equilibration, and even then there was an oﬀset
time.
• PME electrostatics should always be preferred over
reaction ﬁeld, and cutoﬀ should be avoided altogether.
Even with a long cutoﬀ, such as that tested here, reaction
ﬁeld electrostatics leads to unpredictable behavior; mean
ﬁeld is not a good description for systems containing
pore-like structures.
• Stochastic or local thermostats should be used. Global
thermostats, these apply the same change everywhere at
the same moment, are ﬁne in bilayer systems, but in
systems containing areas with small numbers of water,
such as in pores, they seem to introduce deviations and
lead to creation of artiﬁcial temperature gradients. The
Parrinello−Donadio−Bussi V-rescale thermostat81,82
combines a stochastic component with a global
(Berendsen weak coupling76
) thermostat. Both here
and in other studies40
it has shown excellent behavior. It
is already implemented in Gromacs.65
NAMD83
includes
the Lowe−Andersen thermostat89
that is local and
momentum conserving and has been shown to perform
well in dissipative particle dynamics simulations.89,90
• If an option for neighbor list updates exists, it should be
set to 1 (nlist = 1). In Gromacs the default is nlist = 10,
that is, in our experience, safe in lipid bilayer simulations
but not for systems containing pores.
• Charge groups: This is speciﬁc to the GROMOS force
ﬁeld and Gromacs; to our knowledge, GROMOS
simulation software91
also uses them (we did not have
access to GROMOS software and could not test it).
NAMD83
and LAMPPS92
are able to use the GROMOS
force ﬁeld, but they do not utilize charge groups. We
observed that if neighbor lists are updated at each time
step, the default charge groups in GROMOS 53A519
do
not introduce artifacts. With nlist = 10, however,
signiﬁcant unphysical artifacts were observed. When
creating topologies for diﬀerent systems, great care
should be applied in keeping charge groups small as too
large charge groups can lead to very signiﬁcant errors as
reported by Wong-ekkabut et al.40
Finally, we would like to point out that the software and
force ﬁelds used here, as well as the other ones mentioned here,
have been tested in other studies, and they have been shown to
be capable of producing reliable results.19,58,93,94
Force-ﬁeld
development and reliability are among the core issues in
achieving more reliable simulations.21,95
It is, however, always
the user’s responsibility to check the correctness of the
protocols and methods even if it is time consuming.
Figure 4. Pore radius (A) and water density inside the channel (B). Red dotted line: bulk water density. Blue dot−dashed lines represent the average
phosphate locations of the lipids in the lipid bilayer. Constriction (the narrowest pore radius) is indicated by the green dotted line.
Journal of Chemical Theory and Computation Article
dx.doi.org/10.1021/ct3001359 | J. Chem. Theory Comput. 2012, 8, 2905−29112909
■ AUTHOR INFORMATION
Corresponding Author
*E-mail: jirasak.w@ku.ac.th; mikko.karttunen@uwaterloo.ca.
Notes
The authors declare no competing ﬁnancial interest.
■ ACKNOWLEDGMENTS
We would like to thank Dr. Markus Miettinen for helpful
discussions. Financial support was provided by the Thailand
Research Fund (TRF; J.W.), Kasetsart University Research and
Development Institute (KURDI; J.W.), the Faculty of Science
at Kasetsart University (J.W.), the Natural Sciences and
Engineering Research Council of Canada (M.K.), and the
Ontario Early Researcher Award program (M.K.). Computational
resources were provided by SHARCNET (www.sharcnet.
ca) and the Department of Physics, Faculty of Science,
Kasetsart University.
■ REFERENCES
(1) Klepeis, J. L.; Lindorff-Larsen, K.; Dror, R. O.; Shaw, D. E. Curr.
Opin. Struct. Biol. 2009, 19, 120−127.
(2) Schlick, T.; Collepardo-Guevara, R.; Halvorsen, L. A.; Jung, S.;
Xiao, X. Q. Rev. Biophys. 2011, 44, 191−228.
(3) Lyubartsev, A. P.; Rabinovich, A. L. Soft Matter 2011, 7, 25−39.
(4) Hummer, G.; Rasaiah, J. C.; Noworyta, J. P. Nature 2001, 414,
188−190.
(5) Busch, S.; Smuda, C.; Pardo, L. C.; Unruh, T. J. Am. Chem. Soc.
2010, 132, 3232−3233.
(6) Falck, E.; Rog, T.; Karttunen, M.; Vattulainen, I. J. Am. Chem. Soc.
2008, 130, 44−45.
(7) Holt, J. K.; Park, H. G.; Wang, Y. M.; Stadermann, M.; Artyukhin,
A. B.; Grigoropoulos, C. P.; Noy, A.; Bakajin, O. Science 2006, 312,
1034−1037.
(8) Stimson, L.; Dong, L.; Karttunen, M.; Wisniewska, A.; Dutka, M.;
Rog, T. J. Phys. Chem. B 2007, 111, 12447−12453.
(9) Repakova, J.; Holopainen, J. M.; Karttunen, M.; Vattulainen, I. J.
Phys. Chem. B 2006, 110, 15403−15410.
(10) Lippert, R. A.; Bowers, K. J.; Dror, R. O.; Eastwood, M. P.;
Gregersen, B. A.; Klepeis, J. L.; Kolossvary, I.; Shaw, D. E. J. Chem.
Phys. 2007, 126, 046101.
(11) Martinez-Seara, H.; Rog, T.; Karttunen, M.; Reigada, R.;
Vattulainen, I. J. Chem. Phys. 2008, 129, 105103.
(12) Ono, S.; Kuroda, M.; Higo, J.; Kamiya, N.; Nakajima, N.;
Nakamura, H. J. Biol. Phys. 2002, 28, 427−437.
(13) Hegefeld, W. A.; Chen, S. E.; DeLeon, K. Y.; Kuczera, K.; Jas, G.
S. J. Phys. Chem. A 2010, 114, 12391−12402.
(14) Soares, T. A.; Daura, X.; Oostenbrink, C.; Smith, L. J.; van
Gunsteren, W. F. J. Biomol. NMR 2004, 30, 407−422.
(15) Schmid, N.; Eichenberger, A. P.; Choutko, A.; Riniker, S.;
Winger, M.; Mark, A. E.; van Gunsteren, W. F. Eur. Biophys. J. 2011,
40, 843−856.
(16) Auffinger, P.; Cheatham, T. E.; Vaiana, A. C. J. Chem. Theory
Comput. 2007, 3, 1851−1859.
(17) Best, R. B.; Buchete, N. V.; Hummer, G. Biophys. J. 2008, 95, 7−
9.
(18) Cao, Z.; Liu, L.; Wang, J. J. Biomol. Struct. Dyn. 2011, 29, 527−
539.
(19) Oostenbrink, C.; Villa, A.; Mark, A. E.; Van Gunsteren, W. F. J.
Comput. Chem. 2004, 25, 1656−1676.
(20) Ulmschneider, J. P.; Ulmschneider, M. B. J. Chem. Theory
Comput. 2009, 5, 1803−1813.
(21) Piana, S.; Lindorff-Larsen, K.; Shaw, D. E. Biophys. J. 2011, 100,
47−49.
(22) Anezo, C.; de Vries, A. H.; Holtje, H. D.; Tieleman, D. P.;
Marrink, S. J. J. Phys. Chem. B 2003, 107, 9424−9433.
(23) Norberg, J.; Nilsson, L. Biophys. J. 2000, 79, 1537−1553.
(24) Karttunen, M.; Rottler, J.; Vattulainen, I.; Sagui, C. Curr. Top.
Membr. 2008, 60, 49−89.
(25) Fennell, C. J.; Gezelter, J. D. J. Chem. Phys. 2006, 124, 234104.
(26) Feller, S. E.; Pastor, R. W.; Rojnuckarin, A.; Bogusz, S.; Brooks,
B. R. J. Phys. Chem. 1996, 100, 17011−17020.
(27) Patra, M.; Karttunen, M.; Hyvonen, M. T.; Falck, E.;
Vattulainen, I. J. Phys. Chem. B 2004, 108, 4485−4494.
(28) Hardy, D. J.; Stone, J. E.; Schulten, K. Parallel Comput. 2009, 35,
164−177.
(29) Yonetani, Y. Chem. Phys. Lett. 2005, 406, 49−53.
(30) Patra, M.; Karttunen, M.; Hyvonen, M. T.; Falck, E.; Lindqvist,
P.; Vattulainen, I. Biophys. J. 2003, 84, 3636−3645.
(31) Tobias, D. J. Curr. Opin. Struct. Biol. 2001, 11, 253−261.
(32) Sagui, C.; Pedersen, L. G.; Darden, T. A. J. Chem. Phys. 2004,
120, 73−87.
(33) Ballenegger, V.; Arnold, A.; Cerda, J. J. J. Chem. Phys. 2009, 131,
094107.
(34) Mark, P.; Nilsson, L. J. Comput. Chem. 2002, 23, 1211−1219.
(35) Gurtovenko, A. A.; Vattulainen, I. J. Chem. Phys. 2009, 130,
215107.
(36) Gargallo, R.; Oliva, B.; Querol, E.; Aviles, F. X. Protein Eng.
2000, 13, 21−26.
(37) Gargallo, R.; Hunenberger, P. H.; Aviles, F. X.; Oliva, B. Protein
Sci. 2003, 12, 2161−2172.
(38) Bonthuis, D. J.; Rinne, K. F.; Falk, K.; Nadir Kaplan, C.;
Horinek, D.; Nihat Berker, A.; Bocquet, L.; Netz, R. R. J. Phys.:
Condens. Matter 2011, 23, 184110.
(39) Bonthuis, D. J.; Horinek, D.; Bocquet, L.; Netz, R. R. Phys. Rev.
Lett. 2009, 103, 144503.
(40) Wong-Ekkabut, J.; Miettinen, M. S.; Dias, C.; Karttunen, M.
Nat. Nanotechnol. 2010, 5, 555−557.
(41) Zuo, G.; Shen, R.; Ma, S.; Guo, W. ACS Nano 2010, 4, 205−
210.
(42) Ni, B.; Baumketner, A. J. Mol. Model. 2011, 1−11.
(43) Bonthuis, D. J.; Falk, K.; Kaplan, C. N.; Horinek, D.; Berker, A.
N.; Bocquet, L.; Netz, R. R. Phys. Rev. Lett. 2010, 105, 209401.
(44) de Groot, B. L.; Grubmuller, H. Science 2001, 294, 2353−2357.
(45) Tajkhorshid, E.; Nollert, P.; Jensen, M. O.; Miercke, L. J.;
O’Connell, J.; Stroud, R. M.; Schulten, K. Science 2002, 296, 525−530.
(46) Zhu, F. Q.; Tajkhorshid, E.; Schulten, K. Biophys. J. 2004, 86,
50−57.
(47) Beckstein, O.; Sansom, M. S. P. Proc. Natl. Acad. Sci. U.S.A.
2003, 100, 7063−7068.
(48) Cisse, I.; Okumus, B.; Joo, C.; Ha, T. J. Proc. Natl. Acad. Sci.
U.S.A. 2007, 104, 12646−12650.
(49) Song, L. Z.; Hobaugh, M. R.; Shustak, C.; Cheley, S.; Bayley, H.;
Gouaux, J. E. Science 1996, 274, 1859−1866.
(50) Gouaux, E. J. Struct. Biol. 1998, 121, 110−122.
(51) Panchal, R. G.; Cusack, E.; Cheley, S.; Bayley, H. Nat. Biotechnol.
1996, 14, 852−856.
(52) Dekker, C. Nat. Nanotechnol. 2007, 2, 209−215.
(53) Kasianowicz, J. J.; Brandin, E.; Branton, D.; Deamer, D. W. Proc.
Natl. Acad. Sci. U.S.A. 1996, 93, 13770−13773.
(54) Aksimentiev, A.; Schulten, K. Biophys. J. 2005, 88, 3745−3761.
(55) Noskov, S. Y.; Im, W.; Roux, B. Biophys. J. 2004, 87, 2299−
2309.
(56) Gu, L. Q.; Cheley, S.; Bayley, H. J. Gen. Physiol. 2001, 118, 481−
493.
(57) Bernstein, F. C.; Koetzle, T. F.; Williams, G. J. B.; Meyer, E. F.;
Brice, M. D.; Rodgers, J. R.; Kennard, O.; Shimanouchi, T.; Tasumi,
M. J. Mol. Biol. 1977, 112, 535−542.
(58) Olubiyi, O. O.; Strodel, B. J. Phys. Chem. B 2012, 116, 3280−
3291.
(59) Kasianowicz, J. J.; Misakian, M. J. Membr. Biol. 2003, 195, 137−
146.
(60) Patra, M.; Karttunen, M. J. Comput. Chem. 2004, 25, 678−689.
(61) Tieleman, D. P.; Berendsen, H. J. C. J. Chem. Phys. 1996, 105,
4871−4880.
Journal of Chemical Theory and Computation Article
dx.doi.org/10.1021/ct3001359 | J. Chem. Theory Comput. 2012, 8, 2905−29112910
(62) Berger, O.; Edholm, O.; Jahnig, F. Biophys. J. 1997, 72, 2002−
2013.
(63) Bachar, M.; Brunelle, P.; Tieleman, D. P.; Rauk, A. J. Phys. Chem.
B 2004, 108, 7170−7179.
(64) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.;
Hermans, J. In Intermolecular Forces; Pullman, B., Ed.; D. Reidel:
Dordrecht, The Netherlands, 1981; pp 331−342.
(65) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. J. Chem.
Theory Comput. 2008, 4, 435−447.
(66) Wohlert, J.; Edholm, O. Biophys. J. 2004, 87, 2433−2445.
(67) van Gunsteren, W. F.; Berendsen, H. J. C. Mol. Simul. 1988, 1,
173−185.
(68) Humphrey, W.; Dalke, A.; Schulten, K. J. Mol. Graphics 1996,
14, 33−38.
(69) Tironi, I. G.; Sperb, R.; Smith, P. E.; Vangunsteren, W. F. J.
Chem. Phys. 1995, 102, 5451−5459.
(70) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.;
Pedersen, L. G. J. Chem. Phys. 1995, 103, 8577−8593.
(71) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98,
10089−10092.
(72) Ng, J. A.; Vora, T.; Krishnamurthy, V.; Chung, S. H. Eur.
Biophys. J. 2008, 37, 213−222.
(73) Brooks, B. R.; Brooks, C. L.; Mackerell, A. D.; Nilsson, L.;
Petrella, R. J.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch,
S.; Caflisch, A.; Caves, L.; Cui, Q.; Dinner, A. R.; Feig, M.; Fischer, S.;
Gao, J.; Hodoscek, M.; Im, W.; Kuczera, K.; Lazaridis, T.; Ma, J.;
Ovchinnikov, V.; Paci, E.; Pastor, R. W.; Post, C. B.; Pu, J. Z.; Schaefer,
M.; Tidor, B.; Venable, R. M.; Woodcock, H. L.; Wu, X.; Yang, W.;
York, D. M.; Karplus, M. J. Comput. Chem. 2009, 30, 1545−1614.
(74) Case, D. A.; Cheatham, T. E.; Darden, T.; Gohlke, H.; Luo, R.;
Merz, K. M.; Onufriev, A.; Simmerling, C.; Wang, B.; Woods, R. J. J.
Comput. Chem. 2005, 26, 1668−1688.
(75) Limbach, H. J.; Arnold, A.; Mann, B. A.; Holm, C. Comput. Phys.
Commun. 2006, 174, 704−727.
(76) Berendsen, H. J. C.; Postma, J. P. M.; Vangunsteren, W. F.;
Dinola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684−3690.
(77) Kastenholz, M. A.; Hunenberger, P. H. J. Phys. Chem. B 2004,
108, 774−788.
(78) Patra, M.; Hyvonen, M. T.; Falck, E.; Sabouri-Ghomi, M.;
Vattulainen, I.; Karttunen, M. Comput. Phys. Commun. 2007, 176, 14−
22.
(79) Toxvaerd, S.; Dyre, J. C. J. Chem. Phys. 2011, 134, 081102.
(80) Baumketner, A. J. Chem. Phys. 2009, 130, 104106.
(81) Bussi, G.; Donadio, D.; Parrinello, M. J. Chem. Phys. 2007, 126,
014101.
(82) Bussi, G.; Zykova-Timan, T.; Parrinello, M. J. Chem. Phys. 2009,
130, 074101.
(83) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid,
E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kale, L.; Schulten, K. J. Comput.
Chem. 2005, 26, 1781−1802.
(84) Cozmuta, I.; O’Keeffe, J. T.; Bose, D.; Stolc, V. Mol. Simul.
2005, 31, 79−93.
(85) Tang, Y.; Yoo, J.; Yethiraj, A.; Cui, Q.; Chen, X. Biophys. J. 2008,
95, 581−596.
(86) Banerjee, S.; Murad, S.; Puri, I. K. Chem. Phys. Lett. 2007, 434,
292−296.
(87) Hennrich, F.; Arnold, K.; Lebedkin, S.; Quintilla, A.; Wenzel,
W.; Kappes, M. M. Phys. Status Solidi B 2007, 244, 3896−3900.
(88) Huang, B.; Xia, Y.; Zhao, M.; Li, F.; Liu, X.; Ji, Y.; Song, C. J.
Chem. Phys. 2005, 122, 84708.
(89) Koopman, E. A.; Lowe, C. P. J. Chem. Phys. 2006, 124, 204103.
(90) Nikunen, P.; Karttunen, M.; Vattulainen, I. Comput. Phys.
Commun. 2003, 153, 407−423.
(91) Christen, M.; Hunenberger, P. H.; Bakowies, D.; Baron, R.;
Burgi, R.; Geerke, D. P.; Heinz, T. N.; Kastenholz, M. A.; Krautler, V.;
Oostenbrink, C.; Peter, C.; Trzesniak, D.; van Gunsteren, W. F. J.
Comput. Chem. 2005, 26, 1719−1751.
(92) Plimpton, S. J. Comput. Phys. 1995, 117, 1−19.
(93) Oostenbrink, C.; Soares, T. A.; van der Vegt, N. F.; van
Gunsteren, W. F. Eur. Biophys. J. 2005, 34, 273−284.
(94) Cino, E. A.; Wong-ekkabut, J.; Karttunen, M.; Choy, W. Y. PLoS
One 2011, 6, 27371.
(95) Beauchamp, K. A.; Lin, Y. S.; Das, R.; Pande, V. S. J. Chem.
Theory Comput. 2012, 8, 1409−1414.
Journal of Chemical Theory and Computation Article
dx.doi.org/10.1021/ct3001359 | J. Chem. Theory Comput. 2012, 8, 2905−29112911