A new semi-explicit atomistic molecular dynamics simulation method for membrane proteins

Abstract

The implementation and results of a new semi-explicit atomistic molecular dynamics simulation method applied to the potassium channel KcsA in a DPPC lipid bilayer and aqueous electrolyte environment is described in detail. For the non-bonded interactions a new and efficient implementation of the original fast multipole method is presented. To obtain the reaction field due to a large implicit lipid bilayer and aqueous electrolyte surrounding the explicit atomic region, a new implementation of the boundary element method valid for the three dielectric interfaces present is used. Rigorous comparisons with known systems that are exactly or numerically solvable are done to check the validity of the implementation method. With the help of detailed force field maps of the reaction field in the channel pore region due to the large implicit lipid bilayer and with the help of the structural changes seen in KcsA we show the importance of this technique. Due to the smaller explicit region now needed numerous runs on large systems are achievable in a reasonable time ensuring good statistics.

Keywords

Boundary element method fast multipole method reaction field Poisson-Boltzmann equation Debye constant implicit and explicit regions Helmholtz layer molecular dynamics dielectric interfaces membrane protein

1. Introduction

In a large class of simulational problems in biological physics the effects of a large dipolar lipid bilayer need to be incorporated accurately. These effects certainly cannot be neglected if there is an external potential applied across the membrane during the simulations. Except for the special case when the cell membrane is depolarized there always exists a potential across the membrane. In trying to simulate large enough lipid bilayers explicitly as the size of this explicit lipid bilayer increases the systems can become very large, nearing 10 ${}^{6}$ explicit atoms. It is a well-known fact [1] that in simulations of such large systems there are several problems. These include exponentially increasing equilibrium and non-equilibrium relaxation times, correlation times and lengths [1]. In very large systems there are problems like critical slowing down and other finite size effects that can take on special significance [1]. Further with any simulation the full characterization must have reliable error estimates. This can only be obtained after several runs from a number of different initial starting configurations. Thus if a method exists that can reduce the amount of explicit lipid and solvent needed by using implicit regions, more efficient simulations with good statistics and reproducibility are possible. The method can also be used to estimate the amount of explicit lipid/electrolyte needed before large scale simulations are undertaken. To demonstrate an example of this technique we embed the potassium channel protein (KcsA) in a pore within a lipid bilayer (DPPC) surrounded by aqueous electrolyte, and extend this explicit system into a continuum implicit one an order of magnitude bigger. Using the techniques described below we now find the reaction field due to this large implicit region and include its effects in the simulation.

2. Fast multipole method

Efficient calculation of the non-bonded interactions is extremely important in large atomic systems. An adaptive implementation of the fast multipole method (FMM) based on dynamic box sizing and summing of the moments on the fly is used for the explicit atomic level system. For the given system we find the optimal set of parameters to use and ensure that these parameters conserve the total energy of the system without a drift during the entire length of the simulation.

In FMM implementations [2, 3, 4, 5, 6, 7] there are certain basic similarities. These include the initial partitioning of the system into boxes, then three basic operators called A, B and C [2, 3] that translate mulipoles in space are used, and the final calculation of the total potential as a sum of the near and far field potentials is performed. The exact solution where the potential is calculated explicitly (complexity $O(N^{2})$ , where $N$ is the number of explicit atoms in the system) is also available. FMM implementations in [2, 3, 4, 5, 6, 7] also have subtle differences. In the original technique of Greengard and Rokhlin [2] that first gave a general formulation of this method, the multipole expansions were based on terms of $1/r^{k}$ , where $r$ is the point where the potential is desired and $k$ is the number of multipoles in the expansion. Based on this method [2], several other techniques were later developed that used associated Legendre polynomials, spherical harmonics or plane wave expansions for the mulipole moments [3, 4, 5, 6, 7]. The implementation in [3] is based on associated Legendre polynomials, and in that one [4] spherical harmonics are used, both giving $O(N)$ scaling. Although the FMM implementation in [4] is for continuous charge distributions, the modification is a dynamic adjustment of the number of multipoles that are used locally. This does show improvements for continuous charge distributions seen in quantum systems [4]. In the case of explicit atomic systems, where the total number of particles in the system is conserved (NVT or NPT ensemble), the number of multipoles to be used for acceptable accuracy needs to be optimized globally, not locally. For the FMM implementation in [4], the scaling is $O(N)$ . An adaptive FMM implementation based on spherical harmonics is used in [5]. The number of particles is not fixed in the adaptive FMM algorithm in the lowest level boxes. Here again $O(N)$ scaling is seen. An implementation for screened Coulomb interactions using plane wave expansions and spherical harmonics is presented in [6] again with $O(N)$ scaling. Thus, every method has its accuracy and speed, and all the FMM implementations give the same $O(N)$ scaling.

After the original representation of FMM in [2], the $1/r^{k}$ implementation received little attention. The reasons are twofold. First, when a large number of multipoles are needed and there are wide fluctuations in the particle positions, instabilities due to overflow or underflow errors can occur because of variation in the $1/r^{k}$ terms in the multipole expansions. Second, there is a translation of the multipole moments from both finer to coarser levels and vice versa. If dynamic boxes are used, there can be wide fluctuations and instabilities in the $1/r^{k}$ terms. While the concept of dynamic boxes is explored in [5], it is not relevant to the $1/r^{k}$ implementation and it does not consider the important case when the number of atoms in each of the finest boxes is kept constant. The reasons for introducing the dynamic boxes are also two. First, to optimize or maximize the number of finest boxes at the lowest level that are empty when the system is inhomogeneous, and second, to ensure that each of the finest boxes at the lowest level have either a fixed number of particles or that within some predetermined range. In the first case efficient partitioning is ensured, as empty boxes are skipped, and in the second case translation of the moments from the finest lowest level boxes is uniform in terms of the errors introduced due to the truncation of the multipole expansion to a fixed number. As compared to the other techniques [3, 4, 5, 6, 7] based on spherical harmonics, we find that for large explicit atomic systems this $1/r^{k}$ method is indeed faster by at least a factor of two. Our implementation based on the calculation of multipole moments using $1/r^{k}$ terms following the algorithm of [2] with certain modifications is called FMM-1. The implementation based on calculation of multipole moments using associated Legendre polynomials following the algorithm of [3] with the same modifications we call FMM-2. We compare these implementations and also define the optimal parameter set to use for our problem. The general implementations are the same in both the cases. The differences are in the actual calculation of the multipole moments.

2.1 Handling errors in the $1/r^{k}$ method

Our implementation of the $1/r^{k}$ method for explicit atomic systems required only a slight modification of the algorithm presented in [2]. In explicit atomic systems implementing the translational A, B and C operators [2] can cause widely varying values in the $1/r^{k}$ terms leading to overflow or underflow errors. The A, B and C operators presented in detail below contain the terms that depend on $r_{0}^{+n}$ and $r_{0}^{-n}$ , $r_{0}$ where is the distance between the parent and child boxes obtained after partitioning the system, and $n$ is the extent of the multipole moment. In some cases these terms occur twice in a summation and use coefficients from the preceding operators that also have such terms. Thus the overall variations are rather of the form $r_{0}^{+3n}$ and $r_{0}^{-3n}$ . However, the higher order $1/r^{k}$ terms in each of the operators can be easily found from the lower order terms and their repeated multiplications are very efficient. For example, if $k$ is equal to 3, then one needs only the term $1/r$ . The $1/r^{2}$ term is a multiplication of $1/r$ with $1/r$ , and $1/r^{3}$ is another multiplication of $1/r^{2}$ with $1/r$ and so on. Each moment and its coefficients is summed and saved on the fly. In the calculation of multipole moments in non- $1/r^{k}$ implementations [3, 4, 5, 6, 7] based on spherical harmonics this repeated multiplication is not possible. Thus, if these overflow and underflow errors are prevented, this $1/r^{k}$ original FMM implementation of Greengard and Rokhlin [2] is an excellent framework for large explicit atomic systems.

An automatic normalization precluding overflow and underflow errors is developed. The distances are normalized using the diagonals of the dynamic boxes, so that for higher order terms there is no overflow or underflow error. Suppose $p$ is the limit of the number of multipoles included in the expansion. Then for a large $p$ , if we have boxes of varying size in the effort to keep the number of atoms in the lowest level box approximately fixed, we can see that there may be overflow or underflow errors. For example, let the distances $r_{0}$ between a parent and its three child boxes be 2, 4 or 2 ${}^{2}$ and 8 or 2 ${}^{3}$ units. Then the multipole moments vary from $r_{0}^{+3p}$ to $r_{0}^{-3p}$ in the case of $r_{0}=$ 2, from $r_{0}^{+6p}$ to $r_{0}^{-6p}$ in the case of $r_{0}=$ 4 or 2 ${}^{2}$ , and from $r_{0}^{+9p}$ to $r_{0}^{-9p}$ in the case of $r_{0}=$ 8 or 2 ${}^{3}$ . If $p=$ 16, then for $r_{0}=$ 8 or 2 ${}^{3}$ we have the extent of the terms varying from 144 to $-$ 144. Thus, for every box at every level these distances need to be fixed by normalization, and they can later be easily multiplied out from the calculated coefficients. In order to determine the normalization constant, the general method would be to find the highest $r_{0}^{+n}$ and lowest $r_{0}^{-n}$ mulipole terms for each of the operators for the given box. The constant is then selected so that for any value of $p$ all the terms in the mulipole expansion are well-defined without any type of overflow or underflow errors. However, for up to $p=$ 16 we find in our case that the length of the box diagonal works as a good normalization constant. For the example considered above the distances now become $2/2\sqrt{3}$ , $4/4\sqrt{3}$ and $8/8\sqrt{3}$ , i.e. they are $1/\sqrt{3}$ . Thus, the system is equivalent to a system where at any given level all the boxes are of the same size, and a constant box partitioning is done as in [2, 3]. This is new in our implementation. The moments now vary from $r_{0}^{+3p}$ to $r_{0}^{-3p}$ , and for the case of $p=$ 16 we have the extent of the terms varying from 48 to $-$ 48. Float precision can now be used without any loss of accuracy. Once this is checked for the given system size and for $s$ and $p$ values, it can be excluded from the calculations and the length of the box diagonal can be used as the normalization constant. For the case of the FMM-2 implementation this normalization is not needed [3] as all the moments are expanded in terms of sines and cosines that associated Legendre polynomials are expressed as, and which are bounded between 0 and 1. The same is true for any other representation like spherical harmomics [6] or plane wave expansions [7].

2.2 Implementation

Our implementation is a four pass process with an initialization step. We give the local mulipole expansion and the A, B and C operators for the FMM-1 implementation. Those for FMM-2 are well described in [3] and hence are not reproduced here.

2.2.1 Initialization step

Let $p$ be the number of multipoles included in the expansion, $l$ the number of box levels from the coarsest level that involves the entire system to the finest level (lowest level). $s=l-1$ gives the number of child levels. Let there be $m$ explicit charges in each of the finest level boxes in the system. The value of $m$ can vary for each box at the finest level but it is usually fixed in some range. Assume these $m$ charges are located at points $r_{i},i=1,2,3,\ldots,m$ , with $|r_{i}|<r$ for the given lowest level child box. Here $r_{i}$ is measured from the center of each lowest level child box. Then for each lowest level child box the potential expansion at some point $r$ due to all the explicit atoms in that box is given by [2]

$\displaystyle\phi(r)=a_{0}\log(r)+\sum_{k=1}^{p}\frac{a_{k}}{r^{k}}$ (1)

where $a_{0}=\sum_{i=1}^{m}q_{i}$ , $q_{i}$ is the $i$ th charge, and $a_{k}=\sum_{k=1}^{m}\frac{-q_{i}r_{i}^{k}}{k}$ .

The first step in our method is to partition the system using dynamic box sizing, optimizing the number of atoms in each box. In three dimensions each parent will have typically 2 ${}^{3}$ child boxes. The dynamic box sizing allows box sizes to change but keeps the number of atoms in each box approximately constant. The initialization step thus involves the following sub-steps: 1) partition the system dynamically starting from the lowest (finest) level using dynamic boxes of variable size to keep the number of atoms in each box constant or within some predetermined acceptable range. Care must be taken in partitioning the boxes at the lowest level. The box sizes can be varied at this level so that there are a predetermined number of atoms in every box or the boxes are empty. The number of boxes here is a multiple of 2 ${}^{3}$ $\times$ 2 ${}^{3}$ . This will ensure that all parents at every level will have eight children even though the boxes may have variable sizes and that at least three levels are present; 2) for each box at every level find and store the box diagonal that will be used for the normalization (needed only for the FMM-1 implementation); 3) remove empty boxes from lowest level; 4) ignore parents that have all their child boxes empty; 5) using Eq. (1) form the lowest level multipole expansions.

2.2.2 First pass (Operator A)

The A operator forms a $p$ term multipole expansion of the potential of each of the child boxes of a parent shifted to the center of the parent box. The far-field potential at any point $r^{\prime}$ due to all particles contained in a parent box at a given level is

$\displaystyle\phi(r^{\prime})^{A-op}=a_{0}\log(r^{\prime})+\sum_{i=1}^{p}\frac% {b_{i}^{A}}{r^{\prime i}}$ (2)

where

$\displaystyle b_{i}^{A}=\left[\left(\sum_{k=1}^{i}a_{k}r_{A0}^{i-k}\big{(}% \begin{smallmatrix}i-1\\ k-1\\ \end{smallmatrix}\big{)}\right)-\frac{a_{0}r^{i}_{A0}}{i}\right],$

and $r_{A0}$ is the distance between the center of a child box and its parent center. For each level $a_{0}$ remains the same for each box as defined in Eq. (1). $a_{k}$ in Eq. (2) is the A operator coefficient that is the same as in Eq. (1) only at the lowest level but then it is the value of $b_{i}^{A}$ which is stored for the next level calculations. These $b_{i}^{A}$ ’s become the $a_{k}$ ’s in Eq. (2) for the next level calculations and so on. $\big{(}\begin{smallmatrix}i-1\\ k-1\\ \end{smallmatrix}\big{)}$ is the corresponding binomial coefficient. The A operator starts from the lowest level child boxes and proceeds upwards for all the $s$ (child) levels. Thus, after the A operator implementation stage, the parent box at each level has the multipole expansions representing the far-field potential due to all particles contained in the corresponding child boxes at the immediate lower level. The first pass consists of the following sub-steps: 1) starting from the lowest (finest) level use the A operator to translate the multipole expansion of each lowest level box to a parent box; 2) these moments with common origin are summed and stored in parent boxes on the fly; 3) save the $b_{i}^{A}$ coefficients at each level separately.

2.2.3 Second pass (Operator B)

The B operator used at each level between well-separated boxes obtains local expansions in each box due to all particles at that given level that are well-separated from and not contained in the given box. The far-field potential at any point $r^{\prime\prime}$ due to all the particles in a box at a given level is

$\displaystyle\phi(r^{\prime\prime})^{B-op}=\sum_{i=0}^{p}\frac{b_{i}^{B}}{r^{% \prime\prime i}}$ (3)

where

$\displaystyle b_{i}^{B}=\left[\frac{1}{r_{B0}^{i}}\left(\sum_{k=1}^{p}\frac{b_% {k}^{A}}{r_{A0}^{k}}\big{(}\begin{smallmatrix}i+k-1\\ k-1\\ \end{smallmatrix}\big{)}(-1)^{k}\right)-\frac{a_{0}}{ir^{i}_{A0}}\right]$

for $i>$ 0 and

$\displaystyle b_{0}^{B}=a_{0}\log(r_{B0})+\sum_{k=1}^{p}\frac{b_{k}^{A}}{r_{B0% }^{k}}(-1)^{k}.\big{(}\begin{smallmatrix}i+k-1\\ k-1\\ \end{smallmatrix}\big{)}$

is the corresponding binomial coefficient. $r_{B0}$ is the distance between the centers of two well-separated boxes. $b_{k}^{A}$ ’s are the coefficients for the given level that are stored in the first pass. The well separateness condition ( $W=$ 1) is when two child boxes are separated by one box in between and the parents of these boxes are nearest neighbors only. $W=$ 2 is when two child boxes are separated by two boxes in between and the parent boxes are either nearest neighbors or next nearest neighbor and so on. We used three values for $W=$ 1, 2, 3. The values of $b_{k}^{B}$ ’s are saved for each box as separate coefficients called the B operator coefficients. The B operator like the A operator starts from the lowest level child boxes, but unlike the A operator acts at the same level. A further speedup is obtained if one uses symmetry. For example, if a box is well-separated from another box, this other box is also well-separated from it. This means that at a given level for a box $X_{1}$ , only the mulipole expansions in boxes which are well-separated from $X_{1}$ , but are not well-separated from $X_{1}$ ’s parent are translated into local expansions in $X_{1}$ . Since the distances for the well-separated boxes of $X_{1}$ are the same, this can also be done for all of other boxes well-separated from $X_{1}$ . The second pass proceeds from bottom-up with all boxes at each level being considered. After the second pass, all particles in a given box are interacted with particles in other well-separated boxes at the given level in the tree. The errors for $W$ is known exactly [3]. For example, for $W=$ 2 it is given by

$\displaystyle\left|\Phi_{W=2}-\Phi_{\textit{exact}}\right|\leqslant\frac{\sum|% q_{i}|}{(5-\sqrt{3})d}\left(\frac{\sqrt{3}}{5}\right)^{p+1}$ (4)

where 2 $d$ is the grid size. From Eq. (4) it is clear that the errors depend on the number of atoms, the grid size and $p$ for a given value of $W$ . The best values of $W$ to use for the given problem need to be determined and then fixed. Thus, the higher the value of $W$ the greater the accuracy, as the local region where the explicit calculation or the near-field potential is done becomes bigger. The second pass consists of the following sub-steps: 1) use B operator to interact particles in a given box with multipole expansions of each well-separated boxes, summing on the fly; 2) save $b_{0}^{B}$ and $b_{i}^{B}$ coefficients at each level separately.

2.2.4 Third pass (Operator C)

The C operator starts from the coarsest level and expands the potential at the centers of each parent box to that about the center of each child box of that parent. This is done using the coefficients of the parents saved at each level when the B operator was used (B operator coefficients). Let $b_{i}^{B}$ be the B operator coefficients of the parent. Then the B operator coefficients of the child box at the next lower level is updated by adding the term

$\displaystyle c_{n}=\sum_{k=n}^{p}b_{k-n}^{B}\begin{pmatrix}k\\ n\\ \end{pmatrix}(-r_{C0})^{k-n}\text{ where }n=0,1,2,\ldots,p.$ (5)

$r_{C0}$ is the distance between the center of a child box and its parent. This transfers all neglected information from the parent to the child boxes and at the end of the C operation the lowest level child boxes all have local Taylor expansions that can be used to generate the far-field potential due to all particles outside the well-separated child boxes. The third pass consists of the following sub-steps: 1) in a top-down manner use the C operator to transfer each parent’s local Taylor expansion of the far-field potential from the center of the parent to child boxes, again summing the moments on the fly; 2) update the coefficients $b_{i}^{B}$ saved in the second pass with this operation.

2.2.5 Fourth pass (Calculate far- and near-field potentials)

The lowest level boxes at this stage contain a Taylor expansion representing the effects of all particles well-separated from it. Now calculate the far-field potential for each particle using the coefficients $c_{n}$ of the lowest level child boxes from Eq. (5) in Eq. (3). At this stage as the far-field potential is being calculated it is possible to directly calculate the interaction of each particle with the particles in non-well-separated boxes including itself (near-field potential), which is then added to the far-field potential calculated from Eq. (3). This gives the total potential.

Figure 1.

Comparison of FMM-1 and FMM-2 techniques. (a) Plot showing the % deviation from the exact energy of the total system calculated explicitly with that obtained using FMM-1 and FMM-2 with simulation time. (b) Variation of time taken to calculate the non-bonded energy with explicit system size for the exact case, FMM-1 and FMM-2. (c) Variation of the average error of the total system potential energy with explicit system size using FMM-1 and FMM-2. (d) Variation of the average error per particle of the total system potential energy with explicit system size using FMM-1 and FMM-2.

2.3 Optimization and tests of FMM parameters

A first test was to check that there was no drift in the total energy of the system with time during the simulations. Dynamic box sizing was done for a binning range of 0–5, 5–10, 10–15 and 15–20 atoms respectively. The values for $s$ varied accordingly. Three values for $W=$ 1, 2 and 3 were used. The values of $p$ tried were 4, 8, 16 and 20. It was found that for $p=$ 8 or higher and with $W=$ 2 or higher, and up to 10–15 atoms per box, the average total (potential) energy fluctuations with simulation time for the equilibrated system (explained in more detail below) were within 0.1% with no systematic shifts or trends either up or down. Figure 1a gives the case of $p=$ 16, $W=$ 2, a 10–15 atom binning in the lowest level box and $N=$ 60327 explicit atoms for the test system. The test molecular dynamic runs (described later) were for 2 ns in 1 fs steps. As the box size increased so did the number of atoms in each box. A coarser box resulted in slight increases in the total system error given by

$\displaystyle\delta\phi_{\textit{error system}}=\frac{\sum_{i=1}^{N}\phi_{i-% \textit{Exact}}-\sum_{i=1}^{N}\phi_{i-\textit{FMM}}}{\sum_{i=1}^{N}\phi_{i-% \textit{Exact}}}\times 100\%$ (6a) with a corresponding improvement in time. The error per particle given by $\displaystyle\delta\phi_{\textit{error/particle}}=\frac{\sum_{i=1}^{N}|\phi_{i% -\textit{Exact}}-\phi_{i-\textit{FMM}}|}{\sum_{i=1}^{N}\phi_{i-\textit{Exact}}% }\times 100\%$ (6b)

remained largely constant.

In the present problem it was found that above 10–15 atoms per box for $W=$ 2 and $p=$ 16, there was a sudden jump in the average fluctuations of the total system energy from 0.1% to 4%. Above 10–15 atoms per box both the errors given by Eqs (6a) and (6b) showed a jump, as did the average fluctuations of the total system energy, by nearly an order of magnitude. Thus the best value for the number of atoms in a lowest level box for this problem was 10–15, giving the best balance between box coarseness, time improvement, fluctuations in total system energy and errors (Eqs (6a) and (6b)). $W=$ 2, $p=$ 16, and a box density of 10–15 explicit atoms were the values used in all subsequent runs with the FMM implementation. What was more important was the improvement in time by a factor of 2 as compared to FMM-2. Figure 1b compares results for FMM-1 with FMM-2. Improvement for large systems ( $>$ 130 K atoms) is better than 2-fold. The relative errors in the total potential energy of the system having $N$ explicit atoms as defined by Eq. (6a) are less than 2% and potential energy per particle defined by Eq. (6b) are less than 10 ${}^{-5}$ %. Figure 1c gives the variations of the average total system errors Eq. (6a) and Fig. 1d the variations of the average per particle errors Eq. (6b) in the total potential energies of the given system with system size for both FMM-1 and FMM-2. Both FMM-1 and FMM-2 implementation give similar results as far as the errors were concerned although the former is faster in the higher atom limit (Fig. 1b). For the given values of $W$ , $p$ and with 10–15 atoms in the lowest level box we approach the region where the errors defined by Eq. (6a) become less than 2% and by Eq. (6b) less than 10 ${}^{-5}$ % when $N$ (number of explicit atoms) is greater than 50,000. This is seen in Fig. 1c and d. For the case of $N<$ 50,000 another set of optimizations of $W$ , $p$ and lowest level box atom densities is needed. For $N\sim$ 25000 atoms for 10–15 atoms in the lowest level box the errors of the order equivalent to the case for $N>$ 50,000 as shown in Fig. 1c and d can be obtained but $W$ must be increased to 3. This slows the computation mainly because the number of explicit non-bonded calculations increases and this can be considerable for large systems. If W remains 2 then the number of atoms in the lowest level box must be reduced to 0–5 atoms per box to get the same error. This increases the number of boxes and again slows the computation although the latter technique is preferred. However after this as $N$ is reduced still further ( $<$ 15000) the value of $W$ is the only variable that can be increased, and at this stage there is only a slight improvement when using the FMM technique. This further stresses on the fact that the FMM technique yields appreciable improvement in computation time only in the high atom limit ( $>$ 100 K atoms).

While the relative error in the total potential energy of the system appears high it remains the same between any two points in time during the simulation. This error can be reduced to less than 0.1% by decreasing the box size and having 1–3 atoms in the box. However the error being a systematic one is not of concern here. For the current model system ( $N=$ 60327 explicit atoms) and after optimization with $W=$ 2, $p=$ 16, and a lowest level box density of 10–15 explicit atoms, if the differences in the exact total system energies at two points in the simulations are compared with the differences in total system energies calculated using the FMM-1 or the FMM-2 technique at the same two points in the simulations, it is found that these are within 10 ${}^{-3}$ % of each other. Thus in cases where changes in energy at two points of the simulation are needed as in Monte-Carlo or Molecular Dynamics this is acceptable. In situations where the absolute system energy is needed then for a given problem the offset or systematic error needs to be determined for the given values of $W$ , $p$ , and lowest level box occupancy. This can then be used to accurately correct the absolute total system (potential) energy bringing the error down to less than 10 ${}^{-3}$ %. For the case $W=$ 2, $p=$ 16, and a lowest level box density of 10–15 explicit atoms (Fig. 1a) this correction term would typically be $\sim$ $-$ 0.062% for the total system energy obtained by either FMM-1 or FMM-2 techniques.

The FMM-1 method is the faster and the cross-over point is at a system size of $\sim$ 10,000 atoms (Fig. 1b), while both the methods have similar accuracy for the given values of $W$ , $p$ , and box occupancy. Comparison of the current implementation (FMM-1) with that in NAMD or DPMTA in NAMD also showed that in the high atom limit ( $>$ 130 K) FMM-1 gives a more than 2-fold time improvement.

If particles enter or leave the system or if the entire system size changes this technique can also be used. However depending on the rate at which the number of particles or the system size changes the optimal parameters described above may need to be changed on the fly and tests to check the drifts in the total system energy may also be needed. In our case though this does not apply and the optimal parameters were found to be $p=$ 16, $W=$ 2, and a 10–15 box occupancy with the FMM-1 method selected over the FMM-2 one.

3. Boundary element method

In order to incorporate the effects of an extended lipid bilayer and electrolyte region around the explicit system we modify the boundary element method (BEM) [8]. In our application the extended (implicit) region is nearly an order of magnitude larger than the explicit region comprising $\sim$ 60,000 atoms. We extend the BEM method to treat three dielectric cases. Previous work has been limited to two dielectric regions. Figure 2a illustrates the geometry. The explicit domain of dielectric constant $\varepsilon_{1}$ is surrounded by two continuums, one with dielectric constant $\varepsilon_{2}$ is an extension of the explicit lipid bilayer, and the other with dielectric constant $\varepsilon_{3}$ describes the aqueous electrolyte. The red interface depicts the Helmholtz layer introduced to eliminate singularities if the charges come too close to the boundary elements. The implicit region is chosen symmetric about the z-axis, with the $\varepsilon_{2}/\varepsilon_{3}$ interface parallel to the xy-plane. The boundary elements are located at the interfaces between the dielectric regions on the outer side of the Helmholtz layer while explicit system is constrained within the inner side of the Helmholtz layer with reentrant boundary conditions (Fig. 2a). The Helmholtz layer serves two purposes. It can be scaled to minimize errors due to atoms approaching the boundary elements too closely [8]. For sufficiently small boundary elements a well-chosen Helmholtz layer permits direct use of the elemental area of each boundary element. As shown in Fig. 2b the rectangular boundary elements must be carefully sited to treat a three dielectric interface. Each element should abet only two dielectric regions, one on each side. When approximating curved surfaces boundary element normals are usually not parallel to the surface normal at the point of intersection of the element with the surface. In the case of flat rectangular boundary elements at rectilinear interfaces both the normals are parallel greatly reducing the computational overhead. With this geometry solving the integral equations (Appendix A) and integrating over the resulting induced surface charge and dipole densities to obtain the reaction field is trivial, not requiring the use of any complex Gauss curvature area integration techniques. By contrast the method as originally applied to model solvent around molecules [8] can be very expensive for surfaces that are time dependent [9]. Our approach avoids these problems. Also optimization of the size of the Helmoltz layer and of the boundary elements needs to be done only once. This works because of our choice of fixed rectilinear interfaces, and because all explicit molecules are constrained in a single region of dielectric constant also of fixed size. Implementation details are given in Appendix A. To find the optimal size of the Helmholtz layer to be used with the size of the boundary element an exact three dielectric solution (Appendix B) for a similar geometry is used. It is found that for a Helmhotz layer at least twice the size of the boundary element length, the differences of the reaction field energies computed using both the methods are minimal. Far from the explicit region larger boundary elements can be used without sacrificing accuracy. We did not do this. Table 1 provides calculation times needed to determine the reaction field for our model system as a function of the size of the implicit region.

Figure 2.

(a) The entire simulation system after equilibriation. The red region is the Helmholtz layer. The boundary between the implicit regions is shown for Case I1 and extends the average location of the explicit water-lipid interfaces. For Case I2 it is closer and extends the average location of the phosphorous atom in the lipid head groups. In both the cases the planes of the interfaces are symmetric about the z-axis and parallel to it. The lipid bilayer is the center region and the electrolyte regions are the extremes. The x-axis is into the plane of the paper and the z-axis is upwards in the plane of the paper. The y-axis is in the plane of the paper and is from right to left. The origin is at the center of the explicit region. The implicit regions are symmetric about the xy-plane and the z-axis for both Case I1 and Case I2. (b) Incorrect and correct ways of assigning the boundary elements at the three dielectric interfaces. The black lines are the dielectric boundary interfaces and the red lines are the boundary elements.

4. Comparison of BEM with exact solutions

We checked our BEM implementation by comparison with four systems for which either analytical or numerical solutions can be found. Unless otherwise indicated reaction field energies (RFE) ( $q^{\prime}\phi^{\prime}$ ) are computed for a reaction field potential (RFP) $\phi^{\prime}$ and a unit positive test charge $q^{\prime}$ at the origin.

4.1 Case 1

A charge $q$ at the center of a cube of dielectric constant $\varepsilon_{1}$ surrounded by a medium of dielectric constant $\varepsilon_{2}$ (Fig. 3).

Table 1
Variation of the calculation time of the reaction field using BEM technique with the size of the implicit lipid region

Size of implicit lipid region (Å)	BEM (minutes)
100	48
500	184
1000	1273
1500	3741

Table 2

Comparison of the reaction field for Case 1 with that obtained using the BEM technique

Reaction field energy (kcal/mole)	Monopole $q=$ 1
	$q$ at origin	$q$ at $x=y=z=$ 40 Å
$E_{\text{exact}}$	$-$ 6.557	$-$ 6.557
$E_{\text{BEM}}$	$-$ 6.753	$-$ 6.753
% Error	2.900	2.900

Figure 3.

Case 1 with a monopole in a cube surrounded by a dielectric.

RFP in the cube is

$\displaystyle\phi^{\prime}=\frac{-(\varepsilon-1)q}{\varepsilon L/2}\text{ % with }\left(\varepsilon=\frac{\varepsilon_{2}}{\varepsilon_{1}}\right).$ (7)

Table 2 compares the exact RFE (Eexact) from Eq. (7) with the BEM approximation. Unless otherwise indicated, the Helmhotz layer is 2 Å, the grid size is 1 Å and $L=$ 100 Å. From Eq. (7) it is clear that reaction field is independent of the location of the charge $q$ . For this case the BEM method gave similar results when the charge was off-center. This shows that the reaction field depends only on the strength of the source in the central dielectric, in agreement with Gauss’ Law. BEM errors can be reduced (Table 2) by using a smaller grid size.

4.2 Case 2

An infinite slab (dielectric constant $\varepsilon_{1}$ ) surrounded by a medium of dielectric constant $\varepsilon_{2}$ (Fig. 4a).

The RFP is

$\displaystyle\phi^{\prime}=q(1-\varepsilon)\int_{0}^{\infty}\frac{e^{k|z|}J_{0% }(k\rho)dk}{1+\varepsilon e^{kL}}.$ (8)

Table 3 compares the exact reaction field with that obtained using the BEM method choosing $\varepsilon_{1}=$ 1 and $\varepsilon_{2}=$ 2. For the BEM technique the slab must be truncated. We fixed its length at 2000 Å while $L$ is an order of magnitude smaller at 100 Å. We consider centered and off center charge source $q$ . For the off center charge the symmetry is broken and the error in the RFE is sensitive to grid size, as illustrated in Fig. 4b. Optimal grid size is thus critical. In this case the optimal grid size would be $\sim$ 2 Å. By choosing the maximal grid that keeps errors fairly small the reaction field computation can be accelerated without loss of accuracy.

Table 3

Comparison of the reaction field for Case 2 and Case 3 with that obtained using the BEM technique

Reaction field energy (kcal/mole)	$q=$ 1 (placed at origin)		$q=$ 1 (placed at $z=$ 40 Å and $\rho=$ 40 Å)
	Case 2	Case 3	Case 2	Case 3
$E_{\text{exact}}$	$-$ 3.19719	$-$ 7.41452	$-$ 32.23	$-$ 44.0682
$E_{\text{BEM}}$	$-$ 3.2008	$-$ 7.49231	$-$ 32.78	$-$ 48.1245
% Error	$-$ 1.395	$-$ 5.042	$-$ 1.706	$-$ 9.21

Figure 4.

(a) Case 2 with a charge in the center region surrounded by two similar and parallel infinite dielectric slabs. (b) Plot showing the variation of the reaction field energy with BEM size for Case 2 as the charge is fixed off-center.

Figure 5.

Case 3 showing a three dielectric problem. Similar to Case 2 but two different infinite parallel dielectric slabs.

Figure 6.

(a) Case 4 showing a cylinder symmetrically placed with respect to two infinite dielectric regions. The charge is at the origin. The problem is symmetric about the z-axis. (b) Variation of error compared with the exact solution in Case 4 as a function of a finite implicit region size for BEM.

4.3 Case 3

Similar to Case 2 with one important difference. The slab is bounded by two regions of different $\varepsilon$ as illustrated (Fig. 5).

The RFP is then

$\displaystyle\phi^{\prime}=q\left[\int_{0}^{\infty}\xi_{1}(k)e^{+k|z|}J_{0}(k% \rho)dk+\int_{0}^{\infty}\xi_{2}(k)e^{-k|z|}J_{0}(k\rho)dk\right],$ (9)

with

$\displaystyle\xi_{1}(k)=\frac{(1-\varepsilon^{\prime})}{(1+\varepsilon^{\prime% })}\left(\frac{e^{-kL}-1}{Ae^{-kL}-e^{kL}}\right),\xi_{2}(k)=\left(\frac{1-Ae^% {-kL}}{Ae^{-kL}-e^{kL}}\right),\varepsilon^{\prime}=\frac{\varepsilon_{2}}{% \varepsilon_{1}},\varepsilon^{\prime\prime}=\frac{\varepsilon_{3}}{\varepsilon% _{1}}$

and

$\displaystyle A=\frac{(1-\varepsilon^{\prime})(1-\varepsilon^{\prime\prime})}{% (1+\varepsilon^{\prime})(1+\varepsilon^{\prime\prime})}.$

This three dielectric problem has a unique, closed form solution. With $\varepsilon_{1}=$ 1, $\varepsilon_{2}=$ 2 and $\varepsilon_{3}=$ 80, Table 3 presents the errors for RFP for a charge at the origin and off-center. As in Case 2 the errors are sensitive to the location of the charge. Smaller grids will reduce the errors significantly. Large number of charges will cause screening of the inner charges (and far away reaction field induced charges) near a given interface and localize the reaction field. If equal numbers of positive and negative charges are present there will also be a weakening of the source. These effects also help to reduce error. Also for Cases 2 and 3 increasing the length of the slab used in the BEM implementation from 2000 Å will lower the errors as the computational domain will approach the exact infinite slab solution.

4.4 Case 4

A finite cylinder of dielectric constant $\varepsilon_{1}$ in a slab of dielectric constant $\varepsilon_{2}$ surrounded by a dielectric constant $\varepsilon_{3}$ (Fig. 6a).

The RFP is

$\displaystyle\phi^{\prime}=q\left[\int_{0}^{\infty}\left.\xi_{1}(k)e^{+k|z|}J_% {0}(k\rho)dk\right|_{|z|\leqslant L/2}+\int_{0}^{\infty}\left.\xi_{1}(k)e^{+k|% z|}J_{0}(k\rho)dk\right|_{L/2<|z|\leqslant W/2}\right],$ (10)

with

$\displaystyle\xi_{1}(k)|_{|z|\leqslant L/2}=\frac{(1-\varepsilon^{\prime})(e^{% -kL/2}+e^{kL/2})}{2\varepsilon^{\prime}}\left(Ae^{-k(W-L)}-1\right),$ $\displaystyle\xi_{1}(k)|_{L/2<|z|\leqslant W/2}=Ae^{-k(W-L)},\varepsilon^{% \prime}=\frac{\varepsilon_{2}}{\varepsilon_{1}},\varepsilon^{\prime\prime}=% \frac{\varepsilon_{3}}{\varepsilon_{1}}\text{ and }A=\frac{(1-\varepsilon^{% \prime\prime})}{(1+\varepsilon^{\prime\prime})}.$

With $q$ at the origin the solution is symmetric and is found using the boundary conditions at $\rho=R$ , $|z|\leqslant L/2$ and at $\rho=R$ , $L/2<|z|\leqslant W/2$ , which yield equations for $\xi_{1}(k)$ valid in the two regions. This geometry corresponds to a planar lipid bilayer with a symmetric pore. For the case of $L=$ 100 Å, $R=$ 50 Å, $W=$ 200 Å, $\varepsilon_{1}=$ 1, $\varepsilon_{2}=$ 2 and $\varepsilon_{3}=$ 80 the exact solution for the reaction field energy (Eexact $=$ $-$ 4.2748 kcal/mole) is for an infinite dielectric slab. As with Cases 2 and 3 BEM calculations are sensitive to the length of the computational slab, as shown in Fig. 6b for a square planar implicit lipid (dielectric $\varepsilon_{2}$ ) of variable half width, from 100 Å to 1400 Å. The results converge slowly to the exact value given above. This problem is not exactly solvable for an off-center charge as symmetry is broken.

5. Simulational system

H atoms and incomplete residues for KcsA (Protein Data Bank (PDB) entry 1BL8) were generated. The parameter and topographic files of CHARMM force field was used for this. For each residue in KcsA a lookup was done iteratively and at the correct bond lengths the missing residue atoms and the missing H atoms were added. Initially there were overlaps after the first pass between residues. In the second pass each residue with an overlap is rotated keeping bond length fixed so that there is no overlap. Several runs are done in parallel with the structure from the first pass starting from a different chain and chain end to give non-overlapping structures. The energy is calculated for each structure after the second pass (bonded and non-bonded interactions) and the structure with the lowest energy is selected for constructing the system. DPPC bilayer is constructed at correct density for fully hydrated case [10, 11] and equilibrated at 300 ${}^{\circ}$ K. Water molecules are assembled corresponding to density at 300 ${}^{\circ}$ K and equilibrated (300 ${}^{\circ}$ K). Before constructing the system the conjugate gradient (CG) method is used on KcsA and DPPC systems above until change in energy is less than 0.1 kcal/mol (RMSD $<$ 0.1 in both cases.). Then KcsA/Ions, DPPC systems are combined and overlapped DPPC is removed. CG is again done until change is less than 0.1 kcal/mol (RMSD changes being within $\sim$ 0.2.). Then a 2 ns molecular dynamics (MD) run (described in detail below) is done keeping the protein fixed. The DPPC core now has the asymmetric structure of KcsA and clamps it (Fig. 7). $+$ z is out of the plane of the paper. Finally water is added and overlapping waters removed to fully solvate the system. K ${}^{+}$ and Cl ${}^{-}$ atoms are added randomly to give 4 mM extra- and 140 mM intracellular concentrations (35x difference). The Helmholtz layer is 2 Å thick and Boundary Elements are 1 Å square. The entire system (60327 atoms) is illustrated (Fig. 2a) with the dielectric continuum regions and the explicit atomic regions. $+$ x is out of the plane of the paper and origin is at the center of the system. The Helmholtz layer is in red, the dielectric constant in the explicit region is one, with $\varepsilon_{1}=$ 2 and $\varepsilon_{2}=$ 80. The $\varepsilon_{1}/\varepsilon_{2}$ interface can be chosen at the average position of the lipid-water interface (system I1) or shifted to the average position of the phosphorous atoms of the bilayer headgroups (system I2). The boundary separations are 55.8 Å and 40.53 Å respectively. We do this since the dielectric constant of the lipid bilayer varies from 10–15 in the head group region to 2 in the aliphatic core. Thus there is arbitrariness in the siting of the boundary of the implicit bilayer. Macroscopic effects can be studied exploiting this flexibility, and the influence on the reaction field in the protein quantified.

Figure 7.

Top or extracelluar view of the equilibriated simulation system shown in Fig. 2 with water, ions and counter-ions removed. KcsA is in the center clamped by DPPC. The z-axis is out the plane of the paper.

MD simulations were performed in a Canonical Ensemble (NVT) with the velocity Verlet algorithm [12] and a Nosé-Hoover thermostat [13] at 300 ${}^{\circ}$ K. All runs unless indicated were for 2 ns with 1 fs time steps and re-entrant boundary conditions (R-BC) with the implicit region using the BEM technique. Constant external membrane potentials (V) of 0, $+$ 100 mV and $-$ 100 mV were applied across the system. The electric field is perpendicular to the membrane and exerts a small additional force on the atoms in the explicit region. Using $V=E\cdot d$ and $F=qE$ ( $q$ is the partial charge on the atoms, and $d$ is the separation between the upper and lower implicit lipid-water interface) the force ( $F$ ) due to the external potential can be found. This is included for all the atoms in the explicit region located within the upper and lower lipid-water boundary separation of the implicit region. Thus for system I2 (as compared to system I1) as this boundary is nearer a fewer number of explicit atoms are affected by this external potential. For the implicit region a surface charge density of $k\varepsilon_{2}V/d$ is added to the lipid-water interface surface charge and dipole densities $\times$ (A20) from the BEM method for the source of the reaction field. Test runs using a wide combination of external membrane potentials, intra- and extra-cellular K ${}^{+}$ ion concentrations, and implicit region size were conducted. Within the time frame of the current simulations for membrane potentials of $\pm$ 100 mV, 4 mM extra- and 140 mM intra-cellular K ${}^{+}$ ion concentration, and implicit region size of 1000 Å we were able to see structural and large energetic changes in KcsA. The total system size (implicit $+$ explicit) is 1100 Å square. The explicit region ( $\sim$ 100 Å in length) has approximately 4 lipid shells while the implicit region represented approximately 20 lipid shells. Thus an equivalent explicit system would entail at least 24 lipid shells. Note this did not take into account the huge amount of explicit water and ion/counter-ion molecules that would also be required. In every case unless indicated the final KcsA structures obtained were checked for stability by 5 ns runs. The first shell of the lipid bilayer around KcsA (Fig. 7) was updated every 10 fs, and the second and third shells every 100 ps. The fourth shell was fixed. Before each run the total momentum of the system was set to zero. The starting system had three K ${}^{+}$ ions with the waters in the selectivity filter and inner pore regions. There are two waters in the selectivity filter with two K ${}^{+}$ ions. Twenty-seven waters surround the third K ${}^{+}$ in the inner pore and there is one more water below this making a total of 28 waters. This formed the initial starting simulation system. Counter ions do not enter this region. The system was electrically neutral by randomly placing four Cl ${}^{-}$ ions. Nosé-Hoover Thermostat relaxation time was 0.01 ps. The CHARMM force field was used. Averages and error bars where are from 8 different runs each starting from a new CG minimized system. All runs were done on a server with dual E5-2670 v3 Xeon Intel processors @2.3 Ghz with 1 TB RAM and dual Intel Xeon Phi 7110P coprocessors. Code compilation and optimization was done using Intel C++ Compiler 16.0. A 2 ns simulation in 1 fs time steps with the BEM method and with the lipid update time interval as indicated above took between 9.6–10.4 hours. This was because the reaction field using BEM was recalculated only when the RMSD of the protein changed at least 0.1 Å or on an average every 10 ps. Recalculating the reaction field every MD step was not needed due to the small motions in the atoms. The detailed implementation of the optimizations used for the Intel Xeon Phi Coprocessor are beyond the scope of this paper and will be published separately.

Figure 8.

Sample views showing KcsA from top or extracellular side ( $+z$ implies that the z axis is out of the plane of the paper) and from bottom or intracellular side ( $-z$ implies that the z-axis is into the plane of the paper). S is the standard starting configuration, I1 is the system where the dielectric interface is at the average position of the lipid-water interface, and I2 is the system where the dielectric interface is at the average position of the phosphorous atoms of the bilayer headgroups. O-I1 and O-I2 correspond to MD runs on systems I1 and I2 with a $+$ 100 mV membrane potential respectively. C-I1 and C-I2 correspond to MD runs on systems I1 and I2 with a $-$ 100 mV membrane potential respectively.

Table 4

Summary of the results for structures S, O and C for both I1 and I2 cases

Structure	Av. H ${}_{2}$ O in inner	Av. Sel. filter	Av. Inner pore	Energy Change	RMSD change
	pore region	dia. (Å)	dia. (Å)	w.r.t S (%)	w.r.t S (Å)
S	28(2)	4.23	4.49	0	0
I1
O-I1	37(6)	4.45	5.16	$-$ 9.34	3.749
C-I1	24(3)	4.36	2.37	$-$ 11.25	4.258
I2
O-I2	36(5)	4.27	5.07	$-$ 10.73	5.324
C-I2	22(2)	4.34	2.29	$-$ 11.82	3.265

Table 5

Average energy and RMSD changes of C with respect to O for cases I1 and I2

System	Av. Energy change w.r.t O (%)	Av. RMSD change w.r.t O (Å)
I1	$-$ 10.27	2.79
I2	$-$ 8.76	2.28

Figure 9.

Plot showing the characteristic changes in RMSD with respect to C-I1 from one run for Case I1. I1 is the system where the dielectric interface is at the average position of the lipid-water interface, O-I1 correspond to MD runs on I1 with a $+$ 100 mV membrane potential and C-I1 correspond to MD runs on I1 with a $-$ 100 mV membrane potential.

6. Results and discussions

Using this technique we find that KcsA shows appreciable structural mobility and reversible behavior, energetically and dynamically favoring one state to the other. Averaging of the various runs did show the structural and energetic changes were present consistently. Figure 8 shows a representative view of one CG minimized KcsA from a 2 ns MD simulation. These starting configurations were denoted by S. On S further MD simulations with a membrane potential of $+$ 100 mV were then done for the I1 and I2 case. We refer to these structures as O-I1 and O-I2 respectively and a sample view of these are shown in Fig. 8. This can be compared with a similar region in a KcsA CG minimized starting structure S. On structures O-I1 and O-I2, 5 ns runs were done to check for stability. The changes in RMSD were less than 0.01 Å and the changes in the total energy of KcsA (bonded $+$ non-bonded) was less than 0.1 kcal/mole indicating a stable state. After this 2 ns MD runs were then done with a $-$ 100 mV membrane potential. A sample configuration of the conformations obtained in also shown in Fig. 8. We refer to these as C-I1 and C-I2 respectively. Similar 5 ns runs were also done on C-I1 and C-I2. In this case the RMSD were less than 0.007 Å and energy changes less than 0.04 kcal/mole. Indicating this structure could be the more stable one. With the external potential reversed ( $-$ 100 mV) runs from C-I1 gave structures similar to O-I1 but the change was more gradual. Further as can be seen from Fig. 9 the time for the structural change from O-I1 to C-I1 was quicker and the slope is steeper. For I2 results are similar. If further runs were done on the structures C-I1 or C-I2 with the potential reversed ( $+$ 100 mV) structures and energies similar to O-I1 and O-I2 were obtained. For the given set of parameters by changing the potential it was possible to get the system to oscillate between these two structural states. Table 4 gives the energy change of the structures as compared to the starting structure (as a percentage) which is the CG minimized structure of KcsA (S). The average value of this (Fig. 2a shows the top view of one case) was $\sim$ 6.5 $\times$ 105 kcal/mole. Table 4 also gives the average inner pore diameter and the average selectivity filter diameter, as well as the number of H ${}_{2}$ O molecules in the inner pore region below the selectivity filter of KcsA. The number in parenthesis is the number of water in the selectivity region. The sum of the two gives the total number of water molecules in KcsA. Table 4 also gives the corresponding values when system I2 is used. O-I2 and C-I2 are analogous counterparts as that when system I1 is used. In all cases counterions do not enter KcsA and the number of K ${}^{+}$ ions in KcsA remained the same as in S with one in the inner pore and three in the selectivity filter region. Table 5 gives the changes in RMSD and energies between O-I1 and C-I1 and similarly between O-I2 and C-I2. It is now interesting to compare the forces generated by the reaction field in the entire explicit region and in the region of the pore that has KcsA for the two cases of the implicit interface locations. To understand the effects of the changes in the reaction field with the location of the dielectric interface we use force diagrams. A general map in this region is complicated with a myriad of vortices where the lines of force intersect. Plots depicting the full three dimensional reaction force field (RFF) are very complicated due to the innumerable vortices present. This is because the RFF is poloidal [15, 16] with the lines of force twisting over many times in a spiral as one moves along the z-direction of the membrane pore. An easier approach is to study the RFF vectors in important xy planes and along the z direction with respect to the xy plane. Figure 10a shows all the tryptophan residues in KcsA. The upper plane ( $+$ z upwards) is at $z1=$ $+$ 13.5 Å while the lower plane is at $z2=$ $-$ 14.5 Å. Similarly Fig. 10b shows all the glutamate residues in KcsA. The topmost plane here at $z3=$ $+$ 21 Å is the region corresponding to the outer pore slightly above the selectivity filter near the extra-cellular side. The lowest plane at $z4=$ $-$ 24.5 Å is near the inner pore entrance at the intra-cellular side. Finally Fig. 10c shows all the alanine residues in KcsA. The lower plane at $z5=$ 1.5 Å passes through the center of KcsA.

Figure 10.

(a) Location of all the tryptophan residues in KcsA. The upper plane ( $+z$ upwards and in the plane of the paper) is at $z1=$ $+$ 13.5 Å while the lower plane is at $z2=$ $-$ 14.5 Å. (b) Location of all the glutamate residues in KcsA. The topmost plane ( $+z$ upwards and in the plane of the paper) is at $z3=$ $+$ 21 Å and the lowest plane is at $z4=$ $-$ 24.5 Å. (c) Location of all the alanine residues in KcsA. The lower plane at $z_{5}=$ 1.5 Å passes through the center of KcsA.

Figure 11 shows the comparison in the forces due to reaction field in the region of the pore containing KcsA for I1 and I2 cases. For the force vectors shown in Fig. 11 the entire plane has been divided into grids and the force is found for a unit charge at the center of the grid. The grid size is 10 Å $\times$ 10 Å. The force vectors for $+$ 100 mV and $-$ 100 mV external membrane potential is shown for both Case I1 (red arrows) and Case I2 (blue arrows). If the three dimensional force vectors are viewed perpendicular to the xy-plane, then this is the projection of the force vector onto the xy-plane. These force vectors are due to the reaction field boundary charges. Note there is attenuation of the reaction field in the membrane pore due to the explicit region. Within the time-scale of the present MD simulation as the relaxations in KcsA described above are seen only when the reaction field is included it appears that the rather weak forces present due to the reaction field could influence this behavior in KcsA. From Fig. 11 (I1 case) the changes due to reversing the membrane voltage are seen in all planes shown ( $z1$ – $z5$ ). Several local regions of order can be seen. These may be important in causing the structural changes in KcsA seen. In the I2 case (Fig. 11) the difference from I1 is apparent. In the planes $z1$ and $z3$ the forces in the xy plane are very different. The local vertices are fewer and the major forces appear to be vertical with greater degree of orientation among the force vectors. There appears to be a greater focusing of the forces as compared to I1. These effects may translate to the present structural differences seen in KcsA in I1 and I2 mainly in the selectivity filter region (plane $z1$ ). In order to quantitatively compare the changes in the orientations of the force vectors for the two cases (I1 and I2) an order parameter defined by the Maier-Saupe [1] mean field theory

$\displaystyle\text{$<$$m$$>$}=\frac{1}{2}(3<\cos^{2}\theta_{z}>-1)$ (11)

is used. In Eq. (11) $<>$ indicates ensemble averages and $\theta_{z}$ is the angle made by the force vector with the z-azis, and 0 $<$ $<$ $m$ $>$ $<$ 1. A value of 1 indicates all vectors are along the z-axis and a value of 0 is when all vectors are either perpendicular to the z-axis or perfectly random.

Table 6

Average value of $<$ $m$ $>$ for cases I1 and I2 in the pore region, and for the full explicit region. All five planes $z1$ – $z5$ are considered

System	Pore region ( $<$ $m$ $>$ )										Explicit region ( $<$ $m$ $>$ )
	$+$ 100 mV					$-$ 100 mV					$+$ 100 mV					$-$ 100 mV
	z1	z2	z3	z4	z5	z1	z2	z3	z4	z5	z1	z2	z3	z4	z5	z1	z2	z3	z4	z5
I1	0.28	0.29	0.35	0.33	0.24	0.42	0.47	0.35	0.38	0.20	0.16	0.18	0.27	0.29	0.16	0.14	0.17	0.19	0.23	0.12
	(.08)	(.06)	(.09)	(.13)	(.04)	(.05)	(.08)	(.03)	(.06)	(.02)	(.03)	(.05)	(.08)	(.03)	(.02)	(.02)	(.01)	(.04)	(.04)	(.01)
I2	0.65	0.73	0.62	0.71	0.63	0.59	0.62	0.71	0.73	0.64	0.68	0.65	0.67	0.62	0.56	0.62	0.66	0.68	0.61	0.53
	(.10)	(.11)	(.14)	(.12)	(.09)	(.09)	(.10)	(.12)	(.10)	(.08)	(.14)	(.12)	(.10)	(.10)	(.08)	(.14)	(.13)	(.11)	(.13)	(.08)

Numbers in brackets gives the error bars.

Figure 11.

Sample plots of the reaction field forces for I1 (red arrows) and I2 (green arrows) in the five planes $z1$ – $z5$ for membrane potentials $+$ 100 mV (top) and $-$ 100 mV (bottom). $+z$ axis is out of the plane of the paper and the xy-plane is in the plane of the paper. The circular region is a rough location of the pore containing the channel protein KcsA, while the square region shows the exact boundaries to the Helmholtz layer in the explicit region. I1 is the system where the dielectric interface is at the average position of the lipid-water interface, and I2 is the system where the dielectric interface is at the average position of the phosphorous atoms of the bilayer head. The upper plane ( $+z$ upwards) is at $z1=$ $+$ 13.5 Å while the lower plane is at $z2=$ $-$ 14.5 Å (tryptophan residues region – refer to Fig. 10a). $z3=$ $+$ 21 Å is the region corresponding to the outer pore slightly above the selectivity filter near the extra-cellular side (glutamate residues region – refer to Fig. 10b). The lowest plane at $z4=$ $-$ 24.5 Å is near the inner pore entrance at the intra-cellular side (glutamate residues region – refer to Fig. 10b). The lower plane at $z5=$ 1.5 Å passes through the center of KcsA (alanine residues region – refer to Fig. 10c).

Figure 12.

Overlapped conformations O (blue) and C (red) for I1 and I2 cases for comparison. Inset is the S conformation showing the four chains of KcsA. Water, ions, counter-ions and lipids are not shown. Straight green arrows show the relative translational motions between the chains from C to O that occur with the rotations of each of the four chains in KcsA shown with the curved arrows. These changes are reversible as described in the text. I1 is the system where the dielectric interface is at the average position of the lipid-water interface, and I2 is the system where the dielectric interface is at the average position of the phosphorous atoms of the bilayer head. O correspond to MD runs on systems I1 and I2 with $+$ 100 mV membrane potential while C correspond to MD runs on systems I1 and I2 with $-$ 100 mV membrane potential.

Table 6 summarizes the results for both the rough circular region of the pore that contains the channel protein KcsA and for the exact explicit region ending at the Helmholtz layer. The error bars are computed from each of the eight final configurations obtained for the O and C conformations for both the I1 and I2 cases. Several interesting points are worth noting here. For the I1 case all the $<$ $m$ $>$ values were typically less than 0.5 while for the I2 case these were greater than 0.5. Thus there is a greater ordering of the reaction field with respect to the z-axis in the I2 case. A comparison of $<$ $m$ $>$ values between the regions for each case also showed that the numbers for $<$ $m$ $>$ were consistently higher in the pore region than in entire region suggesting further ordering in the inner pore cavity. The focusing is clearly greater for the I2 case. However the other interesting point is the absence of a global symmetry in the field. While close examination of Fig. 11 did show local cyclic ordering of the vectors there was no general symmetry. Thus it appears that on equilibration the lipid bilayer and the channel protein lose symmetry at least about the xy-plane. An inclined planar cut with respect to the z-axis or a more complex surface cut may reveal more global symmetry in these force vectors. In this case the search space is large. The iso-potential maps in the planes $z1$ – $z5$ (not shown) for the two cases I1 and I2 are very complex and yield no further insight. There may exist some symmetry on other complex surface cuts but this is difficult to find. However it is now established that the continuum location from the explicit region to the implicit region of the lipid bilayer is important and must be selected carefully.

Figure 13.

Sample views of KcsA selectivity filter region. A and C chains (top) and D and B chains (bottom) with the three K ${}^{+}$ ions in the selectivity filter. Water in this region is not shown. Counter-ions do not enter this region.

Figure 14.

Close up of the selectivity filter region view from extra-cellular side ( $+z$ ) and intra-cellular side ( $-z$ ) for the I2 case. Water and lipids are removed. Counter-ions do not enter this region.

Figure 15.

Schematic showing the boundaries for the three-element dielectric case. The closed Gauss surface encloses the implicit lipid region and the explicit region. The system is symmetric about the $x y$ plane and the $z$ -axis. The $z$ -axis is upwards and in the plane of the paper, and the $x$ -axis into the plane of the paper. The $y$ -axis is from right to left in the plane of the paper. The origin is at center of the explicit region.

Figure 12 shows the qualitative differences in the location of this region on the response of KcsA. For the I1 case (left) the blue is the O-I1 overlapped with the red which is the C-I1 structure. The overall relaxations can be considered to be one of rotations of the inner trans-membrane (TM) helices (remember DPPC has clamped the outer regions of KcsA preventing large scale motions of the outer TM helices) and horizontal movement of the regions of the inner pore and the selectivity filter. This moving away in effect opens the selectivity filter and the three K ${}^{+}$ ions in S now move away and selectively coordinate with some of the original binding sites. This is shown in Fig. 13. In S two of the three K ${}^{+}$ ions are coordinated with the oxygen in the carbonyls of Y78 and G77 and V76 and T75 respectively. In the case of O-I1 there is a largely horizontal movement coupled with the unscrewing motion of the chains, all three K ${}^{+}$ enter the selectivity filter region and a typical partial coordination is seen. For the third K ${}^{+}$ ions it appears that the hydroxyl oxygen of T75 is also available. In case of C-I2 the change is again largely horizontal with a reversing of the rotational motions and the region in the selectivity filter narrows. There is no change in the coordination of the K ${}^{+}$ ions but as can be seen from Fig. 8 the inner pore in the intra-cellular mouth closes considerably. In the I2 case as can be seen in Fig. 13 the circular movements of the inner TM helices remain but now the selectivity filter moves largely vertically. The selectivity filter is preserved and all the three K ${}^{+}$ ions are tightly bound. The intriguing change is in the inner pore diameter. This is further seen in Fig. 15. For both O-I2 and C-I2 the oxygens of the carbonyls of Y78, G77, V76 and T75 are well presented as prospective binding sites. The hydroxyl oxygen of T75 also is available. Thus in both configurations of the I2 case with the exception of the inner pore diameter the overall structure of the selectivity filter is conserved. To compensate in the closing of the inner pore (C-I2) the selectivity filter moves down in this case. The number of water in the selectivity filter and inner pore for this case is the smallest (Table 4) and this is a compact configuration.

These results demonstrate that the location or the extension of the explicit lipid bilayer to the implicit one is important at least in the case of the present DPPC/KcsA test case. There is a change in the nature of the reaction field when this is done, evident both from the force field maps and the structural changes seen in KcsA. When the external potential is removed and internal and external electrolyte concentrations are equalized both the force field maps (not shown) and KcsA (not shown) do not show appreciable changes and the system is not sensitive to the location of the implicit extension of the explicit lipid bilayer. These results also show that there may be a need for including the effects of a large lipid bilayer that has a voltage and ion concentration gradient across it in studies for similar types of ion channel systems. This inclusion may be via the vastly more computationally expensive explicit formulation of the problem or via a semi-explicit approach presented here. The extent of the lipid bilayer needed in such a case will need to be optimized. Before large scale tests on fully explicit systems are undertaken our technique can provide a quick way to estimate the size of the explicit region that may eventually be adequate.

7. Conclusions

We developed two techniques applicable to both MD and Monte-Carlo (MC) simulations. In the case of the FMM method an almost 2-fold time improvement as opposed to the standard implementations is obtained. The BEM (semi-explicit method) reduces the amount of explicit lipid and water needed to build simulational models that can incorporate the effects of a large lipid bilayer in the pore. If a fully explicit system were used the number of atoms would exceed 10 ${}^{7}$ . This would be equivalent to at least 24 lipid shells and electrolyte around the channel protein in addition to the explicit electrolyte and lipid environment. The problem would also be greatly exacerbated by the large equilibration and relaxation times and correlation time and length scales. As long-range electrostatic force effects are from large areas of the channel’s surroundings the resulting effects on the channel protein should not be ignored. We find that radial forces generated by the reaction field can perturb this system but only if an external potential and electrolyte concentration gradient is present and the BEM technique is used with a large implicit region. What is unclear are the precise effects attributable to the surrounding’s reaction field with an external potential present. With no specific electro-chemical gradient protein conformations were unchanged. Longer simulations and similar studies based on other channel structures may be informative. The most important advantage of the BEM method is the ability to include effects of large lipid-solvent regions around the explicit region without increasing the explicit system size. The optimal size of this inclusion will ultimately depend on the specific systems being studied. Further by comparing the changes in the charge and dipole density distributions during a simulation it is possible to implement this method more efficiently by staggering the times between which the surface charge densities are calculated due to the reaction field and by using an increasing grid size as one moves farther away from the explicit region. Finally the BEM technique is suitable for large scale parallelization as it is a grid based method and uses rectilinear interfaces making the partitioning trivial.

Appendix B

In this three element dielectric method the discrete form of Gauss’ Theorem is used and the symmetries are exploited to obtain the reaction field. This is analogous to the BEM case when the electrolyte concentration is set to zero ( $\kappa=$ 0 in Eq. (A13)). Figure 15 gives a simplified picture of this concept. The Gauss surface is around the lipid/protein implicit/explicit region as shown. The discrete form of Gauss’ law for multiple dielectric media is given by

$\displaystyle\sum_{i=1}^{n}\vec{D}^{i}\bullet\vec{A}^{i}=\sum_{i=1}^{N}q_{i}$ (B1)

$D^{i}$ are the $x$ , $y$ and $z$ components of the electric displacement and $A^{i}$ are the $x$ , $y$ and $z$ components of the area of the surface enclosing charges $q_{i}$ . To clarify the notion assign $i=$ 1, 2 and 3 and these are the $x$ , $y$ and $z$ components of the vector of along $\hat{k}$ , $\hat{i}$ and $\hat{j}$ respectively, where $i=$ 1 is the explicit region of dielectric constant 1, $i=$ 2 is the lipid bilayer continuum region and $i=$ 3 is the water continuum region. Let $\hat{n}_{ij}$ be a unit normal to the surface directed from region $j$ to region $i$ . Then the equations of the electric displacement ( $D$ ), electric field ( $E$ ) and electric polarization ( $P$ ) across the dielectric interface are $(\vec{D}_{i}-\vec{D}_{j})\bullet\vec{n}_{ij}=\sigma_{\text{sc}}$ , $(\vec{E}_{i}-\vec{E}_{j})\times\vec{n}_{ij}=0$ , and $-(\vec{P}_{i}-\vec{P}_{j})\bullet\vec{n}_{ij}=\sigma_{\text{pol}}$ respectively. Thus the tangential component of the electric field is continuous across the interface and the normal component of the electric displacement is also continuous if the bound surface charge density ( $\sigma_{\text{sc}}$ ) is zero. In this case the total surface charge density will equal the polarization surface charge density ( $\sigma_{\text{pol}}$ ). This can be written as

$\displaystyle{}^{n}D^{i}={}^{n}D^{j}\text{ (normal component of electric % displacement ${}^{n}D$ is continuous)},$ (B2)

and

$\displaystyle{}^{t}E^{i}={}^{t}E^{j}\text{ (tangential component of electric % field ${}^{t}E$ is continuous)},$ (B3)

Note that there is no charge deposited on the dielectric interfaces and the only charges are the induced surface charges here due to the reaction field. The Helmholtz layer used in this case prevents any explicit charge to get bound to the dielectric surfaces (keeping $\sigma_{\text{sc}}=$ 0). Also since $\vec{D}=\varepsilon\vec{E},(\varepsilon=\frac{\varepsilon_{i}}{\varepsilon_{0}})$ and $\vec{P}=\vec{D}-\varepsilon_{0}\vec{E}$ the problem reduces to finding all $D^{i}$ satisfying Eq. (B1). The polarization can then be calculated at every point across the interface. It is easy to see from the symmetry of the problem at any point on the interfaces the $x$ , $y$ and $z$ components of the electric displacement, electric field and electric polarization vectors are either parallel or perpendicular to the $x$ , $y$ and $z$ components of $\vec{A}^{i}$ . The normal and tangential components of the electric displacement, electric field and electric polarization are the also the respective $x$ , $y$ and $z$ components. Now if the boundary conditions of the continuity of the normal component of the electric displacement and the tangential component of the electric field are imposed, and using $\vec{D}=\varepsilon\vec{E},(\varepsilon=\frac{\varepsilon_{i}}{\varepsilon_{0}})$ , the solution to obtain $\sigma_{\text{pol}}$ from the reaction field involves solving a set of algebraic equations for the components of the electric displacement on either side of each interface. The simultaneous solution of this set of algebraic equations from the impositions of the boundary conditions with Eq. (B1) will give the components of the electric displacement on either side of each interface. From this the components of the electric polarization are found, which in turn gives the polarization surface charge density. Note also that the volume bound charge given by $\rho_{vb}=-\vec{\nabla}\bullet\vec{P}$ is zero as there is a homogeneous isotropic and uniform medium. Now the calculation of the reaction field from the polarization surface charge density is trivial.

Footnotes

Acknowledgments

The author wishes to thank Intel Corporation and in particular Michael Moretti of Intel Corporation for loaning a large memory dual CPU Intel Xeon Server with dual Intel Xeon Phi Coprocessors on which this work was done.

Conflict of interest

The author declares no conflict of interest.

Appendix A

The BEM method involves transforming a set of differential equations with boundary values into a set of integral equations thereby reducing the dimensionality in the problem by 1. Let $q_{i}$ be the $i^{\text{th}}$ charge with position vector $\vec{r_{i}}$ . Let $N$ be the total number of charges (explicit atoms). In the explicit region Poisson’s equation given by

(A1) $\displaystyle\nabla^{2}\phi(\vec{r})=-\sum_{i=1}^{N}q_{i}\delta(\vec{r}-\vec{r% _{i}})/\varepsilon^{\prime},$

holds, where $\varepsilon^{\prime}=$ 1 for scaling of charges. In the continuum region the linearized form which is one approximation of the Poisson-Boltzmann (PB) equation given by

(A2) $\displaystyle\nabla^{2}\Psi(\vec{r})=\kappa_{i}^{2}\Psi(\vec{r})$

holds, where

(A3) $\displaystyle\kappa_{i}^{2}=\frac{2\textit{IF}^{2}}{\varepsilon_{i}\textit{RT}% },\text{(square of the inverse Debye screening length).}$ (A4) $\displaystyle I=\frac{1}{2}\sum_{i}c_{i}z_{i}^{2}\text{ ($i=$ 1, 2, or 3 % depending on the continuum region)},$

$\varepsilon_{i}=\bar{\varepsilon_{i}}\varepsilon_{0}$ with $\bar{\varepsilon_{i}}$ being the dielectric constant in the region. $c_{i}$ is the concentration of the $i^{\text{th}}$ ionic species, $z_{i}$ its charge, $F$ is Faraday’s constant, $R$ is the Universal Gas constant and $T$ is the absolute temperature. The position vector $\vec{r}$ can be in any of the three continuum dielectric regions in this problem $\varepsilon_{1}$ , $\varepsilon_{2}$ , or $\varepsilon_{3}$ (Fig. ).

$\Psi$ satisfies the regularity condition at infinity. $|\vec{r}|\Psi(\vec{r})$ and $|\vec{r}|^{2}\vec{\nabla}\Psi(\vec{r})$ are bounded for $\vec{r}\rightarrow\infty$ . Hence the problem has a unique solution. The boundary now can be defined as:

(A5) $\displaystyle\phi(\vec{r})=\Psi(\vec{r})$

on the interface/boundary and,

(A6) $\displaystyle\partial\phi(\vec{r})/\partial n=\hat{\varepsilon}\partial\Psi(% \vec{r})/\partial n$

where $\vec{n}$ is the outward unit normal to the surface at $\vec{r}$ on the interface, and $\hat{\varepsilon}=\frac{\varepsilon_{i}}{\varepsilon^{\prime}}$ . In this formulation the difficulty when the source charges approach too close to the interface causing numerically inaccuracies due to strong variations in the near and far field potentials is solved by incorporating a Helmhotz layer and the physical and mathematical basis of this is given in []. Let the boundary/interface be $\Sigma$ . If $\Sigma$ is chosen in such a way so that one component of the outward unit normal in Cartesian coordinates ( $n_{x}$ , $n_{y}$ , or $n_{z}$ ) is always parallel to one component of the boundary element area vector ( $A_{x}$ , $A_{y}$ or $A_{z}$ ), and if the boundary elements now are sufficiently small and rectilinear the use of complicated Gauss curvature integration techniques to approximate a curved surface is avoided. The interface surface normals and the boundary element surface normals are parallel and this creates a simple geometry. The computational overheads to keep track of and correct for the case when these two normals are not parallel is not present.

For points inside the explicit region the position vector is $\vec{r_{0}}^{-}$ , and for points in the continuum region it is $\vec{r_{0}}^{+}$ . Points on the interface/boundary are represented by $\vec{r_{0}}$ . To simplify the representation let $h(\vec{r})=\partial\phi(\vec{r})/\partial n$ . Then the transformation into two coupled integral equations is straightforward. Let

(A7) $\displaystyle F(\vec{r_{i}},\vec{r_{0}})=\frac{1}{4\pi\left|\vec{r_{i}}-\vec{r% _{0}}\right|}$

be the fundamental solution to Eq. () and let

(A8) $\displaystyle P(\vec{r_{i}},\vec{r_{0}})=e^{-K\left|\vec{r_{i}}-\vec{r_{0}}% \right|}F(\vec{r_{i}},\vec{r_{0}})$

be the fundamental solution to Eq. (). Equations (A1) and (A2) also satisfy

(A9) $\displaystyle\nabla^{2}F(\vec{r_{i}},\vec{r_{0}})=-\delta(\vec{r_{i}}-\vec{r_{% 0}})$

and

(A10) $\displaystyle\nabla^{2}P(\vec{r_{i}},\vec{r_{0}})=\kappa^{2}-\delta(\vec{r_{i}% }-\vec{r_{0}})$

respectively. Then multiplying $\phi(\vec{r})$ with Eq. (A9), subtracting $F(\vec{r_{i}},\vec{r_{0}})$ times Eq. (A1) and using Greens second theorem on the volume inside $\Sigma$ and multiplying $\phi(\vec{r})$ with Eq. (A10), $P(\vec{r_{i}},\vec{r_{0}})$ times Eq. (A2) and using Greens second theorem with Eq. () for the volume outside $\Sigma$ gives two equations. These two equations can now be rearranged to

(A11) $\displaystyle\frac{1}{2}(1+\hat{\varepsilon})\phi(\vec{r_{0}})-\iint\limits_{% \Sigma}\left[L_{1}(\vec{r},\vec{r_{0}})\phi(\vec{r})+L_{2}(\vec{r},\vec{r_{0}}% )h(\vec{r})\right]d\vec{r}=\sum_{i}q_{i}F(\vec{r_{i}},\vec{r_{0}})$

and,

(A12) $\displaystyle\frac{1}{2}(1+\frac{1}{\hat{\varepsilon}})h(\vec{r_{0}})-\iint% \limits_{\Sigma}\left[L_{3}(\vec{r},\vec{r_{0}})\phi(\vec{r})+L_{4}(\vec{r},% \vec{r_{0}})h(\vec{r})\right]d\vec{r}=\sum_{i}q_{i}\frac{\partial F(\vec{r_{i}% },\vec{r_{0}})}{\partial n_{0}}$

with

(A13) $\displaystyle L_{1}(\vec{r},\vec{r_{0}})=\hat{\varepsilon}\frac{\partial P(% \vec{r},\vec{r_{0}})}{\partial n}-\frac{\partial F(\vec{r},\vec{r_{0}})}{% \partial n}$ (A14) $\displaystyle L_{2}(\vec{r},\vec{r_{0}})=F(\vec{r},\vec{r_{0}})-P(\vec{r},\vec% {r_{0}})$ (A15) $\displaystyle L_{3}(\vec{r},\vec{r_{0}})=g(\vec{r},\vec{r_{0}})$ (A16) $\displaystyle L_{4}(\vec{r},\vec{r_{0}})=\frac{\partial F(\vec{r},\vec{r_{0}})% }{\partial n_{0}}-\frac{1}{\hat{\varepsilon}}\frac{\partial P(\vec{r},\vec{r_{% 0}})}{\partial n_{0}}$

$\vec{n_{0}}$ is the outward unit normal at point $\vec{r_{0}}$ on $\Sigma$ .

$\displaystyle g(\vec{r},\vec{r_{0}})=\frac{\partial^{2}P(\vec{r},\vec{r_{0}})}% {\partial n_{0}\partial n}-\frac{\partial^{2}F(\vec{r},\vec{r_{0}})}{\partial n% _{0}\partial n},$

and the relation

$\displaystyle\lim\limits_{\vec{r}^{\pm}\rightarrow\vec{r}}\iint\limits_{\Sigma% }g(\vec{r},\vec{r_{0}}^{\pm})d\vec{r}=\iint\limits_{\Sigma}g(\vec{r},\vec{r_{0% }})d\vec{r}$

holds and defines the continuity due to the boundary conditions. Thus $g(\vec{r},\vec{r_{0}})$ which contains second derivatives is also integrable over $\Sigma$ . The potentials inside and outside the boundaries are given by:

(A17) $\displaystyle\phi(\vec{r_{0}}^{-})=\iint\limits_{\Sigma}\left[L_{1}(\vec{r},% \vec{r_{0}})\phi(\vec{r})+L_{2}(\vec{r},\vec{r_{0}}^{-})h(\vec{r})\right]d\vec% {r}+\sum_{i}q_{i}F(\vec{r_{i}},\vec{r_{0}}^{-})$ (A18) $\displaystyle\hat{\varepsilon}\Psi(\vec{r_{0}}^{+})=\iint\limits_{\Sigma}\left% [L_{1}(\vec{r},\vec{r_{0}}^{+})\phi(\vec{r})+L_{2}(\vec{r},\vec{r_{0}}^{+})h(% \vec{r})\right]d\vec{r}+\sum_{i}q_{i}F(\vec{r_{i}},\vec{r_{0}}^{+})$

To solve Eqs (A11) and (A12) first partition the surface into $w$ rectangular elements. Each rectangular element satisfies the condition of having its surface area normal parallel to the interface normal. In addition when using this implementation for the three dielectric case each rectangular element must not be positioned so that three or more dielectric surfaces make contact with it. If this happens then the rectangular element needs to be split into two smaller ones, so that each side of it has only one dielectric region as shown in Fig. 2b. With this constraint the method can be applied to any number of dielectric boundaries with the limitation being that the entire boundary element surface has every boundary element having only two continuous and uniform dielectric regions in contact with it, one on each side. Not only does this make the selection of $\phi$ and $h$ for the case of rectangular elements imperative but also makes the integration easy as there are no curved surfaces to be approximated. $\phi$ and $h$ thus approximated by continuous trial functions $\hat{\phi}$ and $\hat{h}$ , and these are polynomials uniquely determined by their values at the nodes. For rectangular elements there are 4 nodes ( $n=$ 4). So $\hat{\phi}(\vec{r})=\sum_{k=1}^{n}a_{k}S_{k}(\vec{r})$ and $\hat{h}(\vec{r})=\sum_{k=1}^{n}b_{k}S_{k}(\vec{r})$ and $S_{k}(\vec{r})$ is trial basis function that is 1 at the $k^{\text{th}}$ node and 0 elsewhere. To determine the coefficients $a$ and $b$ the collation method is used []. Here the coefficients are selected so that the trial functions satisfy the integral equations exactly at the collation points ( $w$ points). The integral equations can now be written as sum of integrals over elements giving a set of linear equations:

(A19) $\displaystyle\begin{pmatrix}a\\ b\\ \end{pmatrix}-\begin{pmatrix}A_{1}&A_{2}\\ A_{3}&A_{4}\\ \end{pmatrix}\begin{pmatrix}a\\ b\\ \end{pmatrix}=\begin{pmatrix}c\\ d\\ \end{pmatrix}$

$A_{1}$ – $A_{4}$ are matrices of local surface integrals involving kernels $L_{1}$ to $L_{4}$ and the RHS of Eq. () are the source terms in the integral equations. This equation can be rewritten in a single matrix equation

(A20) $\displaystyle(I-A)x=p$

of order 2 $w$ which can be solved by LU decomposition.

The vector $x$ in Eq. () represents the charge and dipole surface densities and serves as the source of the reaction field. The internal potential is now given by

(A21) $\displaystyle\phi(\vec{r_{0}}^{-})=\sum_{j=1}^{2n}t_{j}x_{j}+\phi^{C}(\vec{r_{% 0}}^{-})$

where the $t_{j}$ ’s are the local surface integrals representing the potential of the local surface charge and dipole densities and $\phi^{C}$ is the direct Coulomb interaction. This gives the full potential where the first term in Eq. () is that due to the reaction field and the second the explicit Coulomb interaction to which the other bonded and non-bonded terms can be added. A similar equation holds for the external potential but is not of interest here.

References

Stanley

H.E.

, Introduction of Phase Transitions and Critical Phenomena, Oxford University Press, New York, 1971.

Greengard

and Rokhlin

, A fast algorithm for particle simulations, J Comput Phys 73 (1987), 325.

White

C.A.

and Head-Gordon

, Derivation and efficient implementation of the fast multipole method, J Chem Phys 101 (1994), 6593.

Schmidt

K.E.

and Lee

M.A.

, Implementing the fast multipole method in three dimensions, J Stat Phys 63 (1991), 1223.

Rudberg

and Salek

, Efficient implementation of the fast multipole method, J Chem Phys 125 (2006), 084106.

Cheng

Greengard

and Rokhlin

, A fast adaptive multipole algorithm in three dimensions, J Comput Phys 155 (1999), 468.

Greengard

and Huang

, A new version of the fast multipole method for screened coulomb interactions in three dimensions, J Comput Phys 180 (2002), 642.

Juffer

A.H.

Botta

E.F.F.

van Keulen

B.A.M.

van der Ploeg

and Berendsen

H.J.C.

, The electric potential of a macromolecule in a solvent: A fundamental approach, J Comput Phys 97 (1991), 144.

Boschitsch

Fenley

M.O.

and Zhou

H.X.

, Fast boundary element method for the linear Poisson-Boltzmann Equation, J Phys Chem B 106 (2002), 2741.

10.

Antonov

V.F.

Anosov

A.A.

Norik

V.P.

and Smirnnova

E.Y.

, Soft perforation of planar bilayer lipid membranes of dipalmitoylphosphatidylcholine at the temperature of the phase transition from the liquid crystalline to the gel state, Eur Biophys J 34 (2005), 155.

11.

Negrete

H.O.

Rivers

R.L.

Gough

A.H.

Colombini

and Zeidel

M.L.

, Individual leaflets of a membrane bilayer can independently regulate permeability, J Biol Chem 271 (1996), 11627.

12.

Scope

W.C.

Anderson

H.C.

Berens

P.H.

and Wilson

K.R.

, A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters, J Chem Phys 76 (1982), 637.

13.

Hoover

W.G.

, Canonical dynamics: Equilibrium phase-space distributions, Phys Rev A 31 (1985), 1695.

14.

Gao

Sotomayor

Villa

Lee

E.H.

and Schulten

, Molecular mechanisms of cellular mechanics, Phys Chem Chem Phys 8 (2006), 3692.

15.

Groves

J.T.

Boxer

S.G.

and McConnell

H.M.

, Electric field-induced critical demixing in lipid bilayer membranes, Proc Natl Acad Sci 95 (1998), 935.

16.

McWhirter

J.L.

Ayton

and Voth

G.A.

, Coupling field theory with mesoscopic dynamical simulations of multicomponent lipid bilayers, Biophysical Journal 87 (2004), 3242.

A new semi-explicit atomistic molecular dynamics simulation method for membrane proteins

Abstract

Keywords

1. Introduction

2. Fast multipole method

2.1 Handling errors in the 1 / r k method

2.2 Implementation

2.2.1 Initialization step

4.1 Case 1

Table 1 Variation of the calculation time of the reaction field using BEM technique with the size of the implicit lipid region

Appendix B

Footnotes

Acknowledgments

Conflict of interest

Appendix A

References

2.1 Handling errors in the $1/r^{k}$ method

Table 1
Variation of the calculation time of the reaction field using BEM technique with the size of the implicit lipid region