Many neuronal diseases such as Alzheimer's dementia are related to a loss of inter-neuron communication. The nicotinic acetylcholine receptor (nAChR) plays a crucial role in this process and is severely affected upon disease progression. Successful therapy approaches rely on modulation of response signals, initiated by the flux of ions through the receptor integrated ion channel at the post-synaptic membrane. We here present the comparison of two methods for calculating the potential of mean force of nAChR mediated ion permeation, in terms of accuracy and performance.
The nicotinic acetylcholine receptor (nAChR) is a cation-selective
channel made up from five homologous subunits symmetrically
arranged around a central pore, thus forming a structure with
fivefold pseudo-symmetry
The nAChR plays a crucial role in neuronal signaling and therefore
holds a key function for learning and memory processes. It is
severely affected by neuronal diseases such as Alzheimer’s dementia
A homology model for the transmembrane part of human α7 nAChR was
constructed on the basis of an electron microscopy structure of torpedo
marmorata (PDB code 1OED)
The topology for the lipid bilayer (POPC (16:0-18:1 Diester PC,
1Palmitoyl-2-Oleoyl-sn-Glycero-3-Phosphocholine) are the same used in
earlier studies
Each pore model was embedded into a hydrated POPC bilayer system, by
removing overlapping lipid and water molecules. After minimization the
systems were equilibrated for 3x1 ns while initial position constraints were
stepwise reduced (protein, backbone, C-alpha). In order to reduce stress
induced by lateral pressure fluctuations each pore model was simulated
during the equilibration phase using surface-tension pressure coupling with
a tension of 37.5 mN/m
US simulations follow the procedure described by Hub and de Groot
The PMFs were constructed with the WHAM procedure of Hub and de
Groot
During pulling SMD semi-isotropic pressure coupling was used applying 1
bar in z-direction while keeping the xy-plane fixed. The ion was placed at
the top or bottom of the simulation system on the central pore axis. The
virtual spring attached to the ion had a constant of 100 kJ/mol/nm2 and was
moved at 0.00375 nm/ns. Two times twelve independent simulations were
carried out, pulling the ion through the pore from either side.
The construction of PMFs follows the method of Anishkin and Sukharev
The error estimate for the SMD PMFs was calculated using block
averaging considering them to be correlated fluctuating quantity. The sets are
divided in a number of blocks and averages are calculated for each block.
The error for the total average is calculated from the variance between
averages of the m blocks Bi as follows: stderr2 = ∑(Bi-<B>)2/(m·(m-1)). The
complete derivation is given in the literature
The simulations were run on a DELL Studio XPS (8 cores, i7 920) and on facilities of the Paderborn Center for Parallel Computing PC2 (http://wwwcs.uni-paderborn.de/pc2/, Arminius (400 cores, Xeon 3.2 GHz EM64T, Infiniband) and Bisgrid (64 cores, dual-core Opteron 2.8 GHz, Infiniband)). Plots and pictures were made with Xmgrace-5.1.22, VMD-1.8.7 and MOE2008.10. The sodium ion permeation through the nAChR was studied with two different methods for the construction of PMFs. The PMFs represent the energy barrier a passing ion has to overcome. This barrier is directly correlated to conductivities accessible by experiments. Therefore PMFs are invaluable tools to access scientific problems such as ion selectivity, channel open-/closing, mutations inside the pore or allosteric modulation.
Basically a PMF describes the probability to encounter the ion at a certain position along the reaction coordinate, compared to the bulk phase (equation 1).
PMF (z) = −kBT ln(P(z) / P0 )
(1)
For US, which can be applied to a broad variety of other problems,
the probability is derived from umbrella histograms of the
restrained atom
The construction of PMFs for pulling SMD follows the Langevin
equation (equation 2) considering the forces exerted on the ion
along the reaction coordinate as irreversible work against the
friction force
Nz PMF (z) = ∑ ( ΔFn − γ Δzn / Δtn )Δzn n=1 (2) Both approaches usually assume that the end points lie within the bulk phase of the solvent. As the ion does not experience any potential at these points, they have to be equal and zero (equation 3).
PMF (zstart ) = PMF (zend ) = 0
(3)
-1) 40
l
o
(Jm30
k
()z 20
G
70
60
50
40
30
20
10
0
The conductive properties and time dependent characteristics of an
ion channel are determined by various features of the pore, as well
as the influence of the extra- and intra-cellular parts. Beside
general aspects like the length of the pore and its diameter, the
orientation and polarity of side chains pointing into the lumen of the pore
are highly relevant. The energy barrier of ion permeation not only
has an enthalpic part but usually has also a large entropic part
US simulations yield robust and well reproducible PMFs as long as a high sampling density can be achieved. This can be a challenging endeavor if the probability to encounter the ion at the gorge portion of a pore is low. Related to this problem is the artificial (partial) loss of the hydration shell of the ion during system preparation, when it is placed in hydrophobic areas of the pore.
Another commonly observed effect for pore US is due to small
sampling errors, that the ends of the PMF are not equal unlike
being proposed in equation 3.. This effect can be overcome with
cyclic correction as described by Hub and de Groot
In this study the nAChR channel was prepared in three slightly different ways, yielding three identical systems in terms of sequence, secondary structure, lipid environment and solvent content. The difference lies within the small conformational diversity. For each of them at least 300 restrained US were carried out. The single 300 ps simulations took 7.8 hours each using 4 cores. The height of the PMFs was determined to be 51.4, 63.1 and 49.5 kJ/mol. The standard deviation for US PMFs determined by bootstrap analysis is in most cases around 2.5 kJ/mol, but in some cases can reach as high as 7.6 kJ/mol (Figure 2). This is largely influenced by the convergence of each simulation. Especially for the gorge portion of the channel larger deviations have to be expected.
An earlier study using US on restrained nAChR-M2 helices
embedded into a bilayer-mimetic slab made from CH4 molecules
reports a barrier of 10.5 kbT for the sodium ion permeation
When an ion is pulled through a pore like the transmembrane part
of the nAChR it follows the natural path the ion would take in a
living organism. As shown in Figure 3 this path may vary
significantly for each individual permeation event. Hence different
distances were covered in different times, yielding slightly different
velocities. Therefore the friction coefficient has to be recalculated
for each simulation separately in order to fulfill equation 3.
The shape of pulling PMFs is slightly deformed as a constant
friction coefficient was used. It was suggested by Gullingrud et. al. to
use a velocity autocorrelation approach in order to apply a
modulated friction coefficient along the reaction coordinate
According to equation 2 different pulling velocities should not
affect the final PMF. In order to evaluate the dependence of the
friction coefficient on different pulling velocities a sodium ion was
pulled through the nAChR pore with 4 different velocities, each
with 6 repetitions. As shown in Figure 4 the absolute height of the
PMFs is not affected within the margin of error. The increasing
minimum around -2 nm has to be contributed to the assumption of
a constant friction coefficient. The higher the velocity of the ion
the less this assumption holds true. As stated above the
deformation of the PMF is related to this limitation. For the highest pulling
velocities of 0.03 nm/ns a flattening of the energy barrier can be
observed. This effect is best described as ‘rupture’. The ion
develops such a high speed that it cannot develop all atomic interactions
with surrounding residues. This can be compared to an
experimental condition under strong electric field called ‘electroporation’
showing analogy to some extend
The single 2100 ps pulling simulations took 18 hours each using 32 cores. In order to achieve a reasonable error estimate 12 repetitions for the influx and 12 for the efflux were prepared. 3.3
The two methods US and SMD are carried out on identical simulation systems, consisting of exactly the same number of atoms with exactly the same atomic interactions. Although the principals for PMF creation are fundamentally different, in terms of the simulation it is the same to constraint the ion to a certain position along the reaction coordinate or the pull it along the same.
As stated above 300 separate simulations of 300 ps length each are
required for the construction of one US PMF. With three
repetitions a total of 9360 CPU hours are required. To construct a
comparable SMD pulling PMF twelve 2100 ps simulations with
6490 CPU hours are needed. On the first glance the SMD method
is one third less expensive than the US method. This conclusion
may only be partially true, strongly depending on the architecture
of available compute resources. As shown in Figure 6 the
simulation setup used in this study scales superlinear on up to 8 cores.
Both HPC clusters are equipped with low latency infiniband
interconnects, while Arminius has two single-core CPUs (Intel XeonDP
3.2 GHz) on each node, Bisgrid has four dual-core CPUs on each
node (AMD Opteron 8220, 2.8 GHz). Although the type of CPU
used in both clusters differs by three generations of technological
development the effective performance in ns/day differs only 30 %
on 8 cores. The sharing of the same system bus by multiple cores
on one node is presumably the main bottle neck hindering better
performance on more nodes. This effect is overlaid by the
influence of the inter node communication via Infiniband thus leading
to the quasi parallel shift observed for large number of cores in
Figure 6. Arminus and Bisgird show a linear speedup on up to 32
and 40 cores, respectively. Using more cores leads to a fall-off of
performance. It has to be stated that the simulation system used in
this studies was not optimized for highly parallel computing.
Aspects like the domain decomposition and grid spacing for PME
electrostatics have a considerable impact, especially when using
larger numbers of cores
To achieve an optimal investment of computational time a high performance per core was anticipated for US simulating. For SMD it was considered more important to minimize the simulated time (wall clock). The relative short US simulations (300 ps) can be finished using just 4 cores within 7.8 hours, while the longer SMD simulations (2100 ps) need 18 hours on 32 cores, both measured on Arminius. Considering this discrepancy and comparing to Figure 6 (lower panel) the relative efficiency (ns per day and core) is the highest on 2 cores for both clusters. The level remains high on up 70 60 50 40 l)o30 m / ()Jk 20 z (G10
0 -10 -20 -30 70 60 50 40 l)o30 m / ()Jk 20 z (G10
0 -10 -20 -30-4 -2 to 8 cores and then drops sharply. The benefits of low latency infiniband interconnects are not used to their full extend. This illustrates clearly that such kind of simulations with GROMACS would be best carried out using distributed computing on multicore machines. 24 cores 24 cores 12 10 y /ad 6 s n 8 4 2 0 50 45 40 35 cores Fig. 6. Performance (upper panel), speedup (middle panel) and relative performance (lower panel) of the transmembrane nAChR system presented in this study. The molecular setup consisting of 71994 atoms has been simulated with GROMACS-4.0.4 on two clusters with low latency Infiniband interconnects (blue – Bisgrid: 4 dual-core AMD Opteron 8220, 2.8 GHz per node; red - Arminius: 2 single-core Intel XeonDP, 3.2 GHz per node). On up to 8 cores superlinear scaling could be observed for both systems. Up to 40 cores the speedup increases nearly linear on Arminius, while the performance breaks in at 32 cores on Bisgrid. In order to achieve an optimal investment of computational time a high performance per core was used for US simulating on 4 cores. For SMD it was considered more important to minimize the simulated time (wall clock) and 32 cores were used. 0 4 8 12 16 20 28 32 36 40 44 48 0 4
US and SMD based PMFs have been investigated to help understanding cation flux of the nAChR ion channel. Both procedures yield comparable energy barriers for the permeation of Na+ through the nAChR. The errors estimated for both methods do not differ considerably. From a technical point of view US tends to be more robust, while SMD enables the differentiation of in- and efflux of ions. Both methods differ by 30 % with respect to the computational cost. Due to the higher number of short simulations required for US this computational procedure is more recommended for low latency distributed computing than SMD.
We thank the PC2 University of Paderborn for providing computer time. JK gratefully acknowledges a research scholarship from the Alexander von Humboldt-Foundation. Support of the e-Science Institute Edinburgh is acknowledged.