Letter pubs.acs.org/JCTC
PACO: PArticle COunting Method To Enforce Concentrations in Dynamic Simulations Claudio Berti,*,† Simone Furini,‡ and Dirk Gillespie† †
Department of Molecular Biophysics and Physiology, Rush University Medical Center, Chicago, Illinois 60612, United States Department of Medical Biotechnologies, University of Siena, viale Mario Bracci 16, I-53100, Siena, Italy
‡
S Supporting Information *
ABSTRACT: We present PACO, a computationally efficient method for concentration boundary conditions in nonequilibrium particle simulations. Because it requires only particle counting, its computational effort is significantly smaller than other methods. PACO enables Brownian dynamics simulations of micromolar electrolytes (3 orders of magnitude lower than previously simulated). PACO for Brownian dynamics is integrated in the BROWNIES package (www.phys.rush.edu/ BROWNIES). We also introduce a molecular dynamics PACO implementation that allows for very accurate control of concentration gradients.
B
edges of the discrete region a (possibly large) fraction of the inserted ions quickly diffuses into the continuum region and are deleted. Here we present a simple and computationally efficient PArticle COunting method (PACO) to provide the desired time-averaged ion concentrations in nonequilibrium particle simulations. With respect to other methods, PACO’s computational time saving is twofold. First, ions are inserted into and deleted from the CCs only when needed, in contrast to the repeated insertion/deletion attempts in DVC and GCMC. Second, PACO needs only ion counting inside the CCsthe first step in the DVC and GCMC methodwithout the subsequent computationally expensive energy calculation (Figure S2 in the Supporting Information). On the other hand, PACO can do exactly what other methods can, including providing very low concentrations. PACO is based on two main ideas: ion insertion/deletion to achieve the desired concentrations and handling of ions that reenter the CCs. We first present the basic ion insertion/deletion idea to introduce the method and to illustrate issues related to concentration boundary conditions in particle simulations. Then we present PACOBD, the implementation for BD that includes the handling of re-entrant ions. Because it is simple and lightweight, PACOBD enables accurate BD simulations with micromolar ion concentrations (3−4 orders of magnitude smaller than previously simulated) for hundreds of microseconds (time scales that are rarely done). Simulations of millimolar and lower ion concentrations are important because they are common in many applications and in biology. PACOBD is integrated in the BROWNIES package5 and is
rownian dynamics (BD) and molecular dynamics (MD) simulations are used to study a large number of particle transport phenomena including ions in bulk electrolytes and their permeation through biological ion channels. Compared to methods based on continuous charge distributions (e.g., Nernst−Planck (NP)), BD and MD preserve the discrete nature of electrolytes, providing a detailed description of the microscopic interactions. For all the advantages of BD and MD, the direct simulation of particles (typically ions) introduces computational complications in the fulfillment of concentration boundary conditions. The most widely used methods to prescribe reservoir concentrations use control cells (CCs). Dual Volume Control (DVC)1 and Grand Canonical Monte Carlo (GCMC),2 for example, are techniques that maintain the desired ion concentrations on a time-averaged basis with repeated ion insertions/deletions during the simulation. The insertions/ deletions are accepted based on the energy of the system and the chemical potential of the ion species. At each time step, a number of ion insertion/deletion are attempted. These techniques are effective in providing the right concentrations (even for low concentrations) and thus have been widely adopted.2,3 However, there are two aspects that are computationally expensive. First, a (possibly large) fraction of the attempted insertion/deletion are rejected, resulting in wasted time. Second, the energy calculation is computationally expensive, comparable to the BD or MD step itself. Another approach to enforce concentrations is a hybrid continuum-particle technique in which the core of the system is modeled in detail with BD or MD and the CCs with continuum approaches like NP theory.4 However, this method can also be inefficient. Specifically, when ions are inserted very close to the © XXXX American Chemical Society
A
DOI: 10.1021/acs.jctc.5b01044 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Letter
Journal of Chemical Theory and Computation
averaged number of ions through step T. (The floor operator ⌊x ⌋ gives the largest integer not greater than x, and the ceiling operator ⌈x ⌉ gives the smallest integer not less than x). The summations in eq 1 are over the entire course of the simulation (i.e., they are not continually recalculated over some smaller block of time steps). To exemplify this basic ion insertion/deletion algorithm, we consider a BD simulation of bulk NaCl under a 40 to 0.4 mM concentration gradient. Figure 2A shows the desired number of Na+ ions ηNa (black line) and the time-averaged number of Na+
available for free download. Finally, we introduce PACOMD, the implementation for MD, that enables more accurate control of ion concentrations of electrolytes in MD simulations. The typical BD/MD simulation domain is a box divided in a transport region (TR) and two CCs (CCL and CCR) (Figure 1A). The CCs provide the desired ion concentrations. In our
Figure 1. (A) Simulation domain of a typical BD and MD simulation of bulk electrolyte and an ion channel. (B) Simulation domain for PACOBD. Within PACOMD each OR (pink regions) is subdivided into reservoir and buffer regions (Figure S6 in the Supporting Information). The ORs are linked with periodic boundary conditions. In both (A) and (B), the membrane region is in gray, and white circles represent ion channel amino acids. The membrane and the channel are not included in the case of bulk electrolytes simulations. The size of the different regions is not to scale.
setup, we extend the simulation box beyond the CCs and introduce two outer regions (ORs, called ORL and ORR) (Figure 1B). Their importance is discussed later. For ion channel simulations, a membrane with a pore is set in the middle of the TR. For simplicity, we limit our discussion to a single CC and a single ion species. The ideal number of ions of species k in a CC is ηk = ckVCC, where ck is the desired concentration of species k and VCC is the CC’s volume. ηk is typically not an integer. In the simulation, at each time step i, one has an integer number of ions in the CC, denoted Nik. This differs from ηk by xik = ηk − Nik. At each time step one would ideally insert (if xik > 0) or delete (if xik < 0) xik ions, but this is not generally possible since xik is a fraction. Moreover, approximating xik as the closest integer is not the solution; if ηk < 0.5, then an ion will never be inserted. Instead of forcing the desired number of ions at each time step to be into the CC, we relax this requirement and consider the sums of ηk and Nik until step T: T
XkT =
T
Figure 2. (A, B, C) BD simulation of a bulk NaCl solution gradient. The left CC (CCL) and right CC (CCR) (white regions) contain 40 mM and 0.4 mM NaCl, respectively. Along z, the CCs, the ORs, and the TR are 1, 5, and 40 Å wide, respectively. Along x and y, the simulation box is 100 Å wide. The simulations were 4 μs long with a time step of 5 fs. Ion parameters are in Table S1 of the Supporting Information, and details of the BD implementation can be found in ref 3. Ion concentrations are the time-averaged number of ions per volume in 1 pm slabs. Ion currents are the number of ions crossing the center plane. Blue lines refer to the basic ion insertion/deletion (i/d) method; red lines refer to the PACOBD method. (A) Evolution of ηNa (black line) and ∑Ti=1NiNa/T during the course of the simulation. (B) Concentration of Na+ ions within CCL and (C) the whole simulation domain. (D) MD simulation of a bulk KCl solution for a 112 to 11.2 mM concentration gradient. In (B), (C), and (D), the percentages are errors with respect to NP theory. Panel B shows only CCL and is a magnified form of the left-hand-side of (C).
T
∑ ηk −
∑ Nki = ηkT −
∑ Nki
i=1
i=1
i=1
(1)
The quantity T T ⎧ ⎪⌊Xk ⌋ , if X k > 0 YkT = ⎨ ⎪ T T ⎩⌈Xk ⌉, if X k ⩽ 0
(2)
XTk
XTk
is the number of ions to insert (for > 0) or delete (for < 0) at time step T such that the CC has the desired timeB
DOI: 10.1021/acs.jctc.5b01044 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Letter
Journal of Chemical Theory and Computation ions in the CC (∑Ti=1NiNa/T, blue line) during the simulation. ∑Ti=1NiNa/T approaches ηNa quickly (∼200 time steps), and the two curves become indistinguishable after just a few thousand time steps (typical BD and MD simulations consist of millions of time steps). This shows that the method provides the right average concentration from the beginning of the simulation with high accuracy (the error is 0.9%). However, it is important to note that it introduces concentration gradients within the CC; the concentration within the CC is not flat (Figure 2B, blue line). The deviation from a constant concentration is quantified by the high (32.1) percent relative standard deviation in concentration (%RSD) within the CC. This produces the wrong ion concentrations at the edges of the CC (Figure 2B, blue line), which in turn produces wrong ion concentrations throughout the simulation domain (Figure 2C, blue line). The resulting error in the current, with respect to NP theory, is 42.2%. In our studies we found that this error varies unpredictably with simulation parameters like the time step length. For example, a 50 fs (versus 5 fs in Figure 2A-C) time step produced a 12.2% error in the current with 11.2%RSD. The low concentrations at the edges of the CCs occur because ions are inserted uniformly within the CC but are typically deleted at the edges of the CC. This happens because when the CC needs to delete an ion of species k, the ion is removed as soon as it crosses the CC edge after diffusing from the TR or the ORs. Many ions that leave the CC quickly reenter and are typically deleted because in the meantime a new ion has been added to replace it. Since ions are inserted uniformly within the CC, but are preferentially deleted at the edges, ions accumulate at the center to give the correct average concentration in the CC (Figure 2B, blue line). For the simulation, what counts is the concentration at the plane facing the TR and having an incorrect concentration there produces large errors (Figure 2C, blue line). This is an issue that affects any method (including DVC and GCMC) that ensures only the average concentration in the CC. Because the concentration gradient is largely caused by a CC deleting an ion it has recently inserted, a simple solution is to prevent each CC from deleting the ions it inserted; ions inserted in the left (right) CC can only be deleted after they have diffused through the simulation domain into the right (left) CC. This is physically justified because it is more like what actually occurs; ions can diffuse anywhere, but a bias is created when ions are deleted on one side of an arbitrary interface but not on the other. By having only the other CC delete an ion, we are minimizing this bias, while at the same time maintaining accurate concentrations. This ion deletion adjustment defines PACOBD (Figure S1 in the Supporting Information). Using only CCL as an example, at each time step T, one must
Unlike basic ion insertion/deletion, when YTk ≤ −1 it is necessary to wait until an ion of species k that was inserted in the right CC gets into the left CC to delete it. Similarly, the right CC can only delete ions inserted in the left CC. The treatment of ions that re-enter the CC is crucial for assuring the right concentration at the CC/TR interface. How ions interact with the opposite side of the CC is also important. For example, if there were a hard wall at the outer edge of the CC with no ion recycling along z (Figure 1A), the ion concentrations would not be flat there. This would throw off the CC/TR concentration when the average concentration in the CC is maintained. For this reason we have placed an OR next to each CC. The OR allows the ions to equilibrate at the CC side opposite the CC/TR interface. Because re-entrant ions are not deleted as they cross either the CC/OR or the CC/TR interfaces, concentration gradients at both ends of the CC are minimized, leading to very accurate CC/TR interface concentrations. In order to keep a continuous current of ions, ORL and ORR are connected with periodic boundary conditions; an ion leaving one OR on the opposite side of the CC/OR interface enters the other OR. Together, the ORs form a transport region analogous to the TR at the center of the simulation domain, but whose purpose is to act as a buffer region between the two CCs to keep them as thermodynamically independent as possible. The concentration gradients in the ORs are not a problem because these ions are far from the TR where the “real” simulation is taking place. PACOBD significantly reduces the concentration gradients within the CCs (%RSD = 1.18, Figure 2B, red line) without affecting the desired ion concentrations (the error is
