Electrostatic field solution

The electrostatic potential is computed at regular time intervals by solving the Poisson’s equation

${\displaystyle \nabla (\varepsilon (r)\nabla \phi (r,t))=-(\rho _{\text{ions}}(r,t)+\rho _{\text{perm}}(r))}$

where ${\displaystyle \rho _{\text{ions}}(r,t)}$ and ${\displaystyle \rho _{\text{perm}}(r)}$ are the charge density of ions and permanent charges on the protein, respectively. ${\displaystyle \epsilon (r)}$ is the local dielectric constant or permittivity, and ${\displaystyle \phi (r,t)}$ is the local electrostatic potential. Solving this equation provides a self-consistent way to include applied bias and the effects of image charges induced at dielectric boundaries.

The ion and partial charges on protein residues are assigned to a finite rectangular grid using the cloud-in-cell (CIC) scheme.^[3] Solving the Poisson equation on the grid counts for the particlemesh component of the P³M scheme. However, this discretization leads to an unavoidable truncation of the short-range component of electrostatic force, which can be corrected by computing the short-range charge-charge Coulombic interactions.

Dielectric coefficient

Assigning the appropriate values for dielectric permittivity of the protein, membrane, and aqueous regions is of great importance. The dielectric coefficient determines the strength of the interactions between charged particles and also the dielectric boundary forces (DBF) on ions approaching a boundary between two regions of different permittivity. However, in nano scales the task of assigning specific permittivity is problematic and not straightforward.

The protein or membrane environment could respond to an external field in a number of different ways.^{[1][4][5][6][7]} Field induced dipoles, reorientation of permanent dipoles, protonation and deprotonation of protein residues, larger scale reorganization of ionized side-chains and water molecules, both within the interior and on the surface of the protein, are all examples of how complicated the assignment of permittivity is. In MD simulations, where all the charges, dipoles, and field induced atomic dipoles are treated explicitly then it is suggested that a dielectric value of 1 is appropriate. However, in reduced-particle ion simulation programs, such as ours, where the protein, membrane, and water are continuum backgrounds and treated implicitly, and on top of that, the ion motion takes place on the same time-scale as the protein’s response to its presence, it is very difficult to assign the dielectric coefficients. In fact, changing the dielectric coefficients could easily alter the channel characteristics, such as ion permeation and selectivity The assignment of dielectric coefficient for water is another key issue. The water molecules inside ion channels could be very ordered due to tapered size of the pore, which is often lined with highly charged residues, or hydrogen bond formation between water molecules and protein.^[8] As a result, the dielectric constant of water inside an ion channel could be quite different from the value under bulk conditions. To make the matter even more complicated, the dielectric coefficients of water inside nanopores is not necessarily an isotropic scalar value, but an anisotropic tensor having different values in different directions.

Anisotropic permittivity

It has become evident that the macroscopic properties of a system do not necessarily extend to the molecular length scales. In a recent research study carried by Reza Toghraee, R. Jay Mashl, and Eric Jakobsson at the University of Illinois, Urbana-Champaign,^[4] they used Molecular Dynamics simulations to study the properties of water in featureless hydrophobic cylinders with diameters ranging from 1 to 12 nm. This study showed that water undergoes distinct transitions in structure, dielectric properties, and mobility as the tube diameter is varied. In particular they found that the dielectric properties in the range of 1 to 10 nm is quite different from bulk water and is in fact anisotropic in nature. Though, such featureless hydrophobic channels do not represent actual ion channels and more research has to be done in this area before one could use such data for ion channels, it is evident that water properties like permittivity inside an ion channel or nano-pore could be much more complicated that it has been thought before. While a high axial dielectric constant shields ion’s electrostatic charges in the axial direction (along the channel), low radial dielectric constant increases the interaction between the mobile ion and the partial charges, or the dielectric charge images on the channel, conveying stronger selectivity in ion channels.

Solving the Poisson equation based on an anisotropic permittivity has been incorporated into BioMOCA using the box integration discretization method,^[9] which has been briefly described below.

Box integration discretization

In order to use box integration for discretizing a D-dimensional Poisson equation

${\displaystyle \nabla (\varepsilon \nabla \varphi )=\rho }$

with ${\displaystyle \varepsilon }$ being a diagonal D × D tensor, this differential equation is reformulated as an integral equation. Integration the above equation over a D-dimensional region ${\displaystyle \Omega }$, and using Gauss theorem, then the integral formulation is obtained

${\displaystyle \oint _{\partial \Omega }{\hat {n}}(\varepsilon \nabla \varphi )=-\int _{\Omega }\rho }$

In this appendix it is assumed to be a two-dimensional case. Upgrading to a three-dimensional system would be straightforward and legitimate as the Gauss theorem is also valid for the one and three dimensions. ${\displaystyle \epsilon }$ is assumed to be given on the rectangular regions between nodes, while ${\displaystyle \varphi }$ is defined on the grid nodes (as illustrated on figure at the right).

Box integration for a two-dimensional tensor product grid. The integration region is indicated by the dashed rectangle. Charges are assumed to be given on the same nodes as potential

The integration regions ${\displaystyle \Omega }$ are then chosen as rectangles centered around node and extending to the 4 nearest neighbor nodes. The gradient ${\displaystyle \nabla \varphi }$ is then approximated using centered difference normal to the boundary of the integration region ${\displaystyle \Omega }$, and average ${\displaystyle \epsilon }$ over the integration surface ${\displaystyle \partial \Omega }$ . This approach allows us to approximate the left hand side of the Poisson equation above in first order as

${\displaystyle \oint _{\partial \Omega }{\hat {n}}(\varepsilon \nabla \varphi )={\frac {\varphi _{i+1,j}-\varphi _{i,j}}{h_{i}^{x}}}\left({\frac {h_{j}^{y}}{2}}\epsilon _{i,j}^{x}+{\frac {h_{j-1}^{y}}{2}}\varepsilon _{i,j-1}^{x}\right)}$

${\displaystyle {}-{\frac {\varphi _{i,j}-\varphi _{i-1,j}}{h_{i-1}^{x}}}\left({\frac {h_{j}^{y}}{2}}\epsilon _{i-1,j}^{x}+{\frac {h_{j-1}^{y}}{2}}\varepsilon _{i-1,j-1}^{x}\right)}$

${\displaystyle {}+{\frac {\varphi _{i,j+1}-\varphi _{i,j}}{h_{j}^{y}}}\left({\frac {h_{i}^{x}}{2}}\varepsilon _{i,j}^{y}+{\frac {h_{i-1}^{x}}{2}}\varepsilon _{i-1,j}^{y}\right)}$

${\displaystyle {}-{\frac {\varphi _{i,j}-\varphi _{i,j-1}}{h_{j-1}^{y}}}\left({\frac {h_{i}^{x}}{2}}\varepsilon _{i,j-1}^{y}+{\frac {h_{i-1}^{x}}{2}}\varepsilon _{i-1,j-1}^{y}\right)}$

where ${\displaystyle \varepsilon ^{x}}$ and ${\displaystyle \varepsilon ^{y}}$ are the two components of the diagonal of the tensor ${\displaystyle \epsilon }$. Discretizing the right-hand side of the Poisson equation is fairly simple. ${\displaystyle \rho }$ is discretized on the same grid nodes, as it's been done for ${\displaystyle \varphi }$.

${\displaystyle \int _{\Omega _{i}}\rho ={\text{Volume}}(\Omega _{i})\rho _{i}}$

Ion size

The finite size of ions is accounted for in BioMOCA using pairwise repulsive forces derived from the 6–12 Lennard-Jones potential. A truncated-shifted form of the Lennard-Jones potential is used in the simulator to mimic ionic core repulsion. The modified form of the Lennard-Jones pairwise potential that retains only the repulsive component is given by

${\displaystyle U_{LJ}(r_{ij})={\begin{cases}4\epsilon _{LJ}\left(\left({\frac {\sigma _{ij}}{r_{ij}}}\right)^{12}-\left({\frac {\sigma _{ij}}{r_{ij}}}\right)^{6}\right)+\epsilon _{LJ}&r_{ij}<2^{1/6}\sigma _{ij}\\0&r_{ij}>2^{1/6}\sigma _{ij}\end{cases}}}$

Here, ${\displaystyle \epsilon _{LJ}}$ is the Lennard-Jones energy parameter and ${\displaystyle \sigma _{ij}=(\sigma _{i}+\sigma _{j})/2}$ is the average of the individual Lennard-Jones distance parameters for particles i and j. Using a truncated form of the potential is computationally efficient while preventing the ions from overlapping or coalescing, something that would be clearly unphysical.

Ion-protein interaction

Availability of high-resolution X-ray crystallographic measurements of complete molecular structures provides information about the type and location of all atoms that forms the protein. In BioMOCA the protein atoms are modeled as static point charges embedded in a finite volume inaccessible to the ions and associated with a user-defined dielectric coefficient. Moreover, a number of force-field parameters are available that provide information about the charge and radii of atoms in different amino-acid groups. The conjunction of the molecular structure and force fields provide the coordinates, radii, and charge of each atom in the protein channel. BioMOCA uses such information in the standard PQR (Position-Charge-Radius) format to map the protein system onto a rectangular grid.

Ideally, the steric interactions between protein atoms and the ions in the aqueous medium are to use a repulsive potential like Lennard-Jones to prevent ions from penetrating the protein. As this approach could add a significant load to the amount of calculations, a simpler approach is chosen that treats the protein surfaces as predetermined hard wall boundaries. Many recent open source molecular biology packages have built-in facilities that determine the volume accessible to ions in a protein system. The Adaptive Poisson Boltzmann Solver (APBS) scheme^[10] has been incorporated to BioMOCA to obtain the accessible volume region and therefore partition the simulation domain into continuous regions.

Ions are deemed to have access to protein and lipid regions and if any point within the finite-size of ionic sphere crosses the protein or membrane boundary, a collision is assumed and the ion is reflected diffusively.

Ion-water interactions

As a reduced particle approach, BioMOCA replaces the explicit water molecules with continuum background and handles the ion-water interactions using BTMC method, in which, appropriate scattering rates should be chosen. In other words, ion trajectories are randomly interrupted by scattering events that account for the ions’ diffusive motion in water.^[1] In between these scattering events, ions follow the Newtonian forces. The free flight times, T_f, are generated statistically from the total scattering rate according to

${\displaystyle -\ln(r)=\int _{0}^{T_{f}}\lambda ({\vec {p}}(t))\,dt}$

where r is a random number uniformly distributed on the unit interval. ${\displaystyle \lambda }$, a function of momentum, is the total scattering rate for all collision mechanisms. At the end of each free flight, the ion’s velocity is reselected randomly from a Maxwellian distribution. As the correct scattering mechanism for ion-water interactions in nonbulk electrolyte solutions has yet to be developed, a position dependent scattering rate linked to the local diffusivity is used in our model. This dependency on position comes from the fact that water molecules can have different order of organization in different regions, which will affect the scattering rate.

Position-dependent diffusivity

It is widely accepted that the ions and water molecules do not have the same mobility or diffusivity in confined regions as in bulk.^[2][6] In fact, it is more likely to have a lessening in the effective mobility of ions in ion channels.^[5] In reduced particle methods where the channel water is assumed as implicit continuum background, a mean ion mobility is needed to reveal how ions could diffuse due to local electrostatic forces and random events. In Transport Monte Carlo simulations, the total scattering rate (${\displaystyle \lambda }$), is assumed to only result from ion-water interactions; it is related to ion diffusivity with the expression

${\displaystyle \lambda ={\frac {kT}{mD}}}$

where m is the mass of the ion and D is its diffusion constant. As the equation indicates, reduced diffusivity of ions inside the lumen of the channel renders to increased incidence of scattering events.

Hydration shells

In addition to having a diffusive effect on ion transport, water molecules also form hydration shells around individual ions due to their polar nature. The hydration shell not only shields the charge on ions from other ions but also modulates the ion radial distribution function causing the formation of peaks and troughs. The average minimum distance between two ions is increased as there is always at least one layer of water molecules present between them, acting as a physical deterrent preventing two ions from getting too close to each other, in a manner that is similar to the short-range repulsive component of the Lennard-Jones potential.

The theory of hydration shells is well developed in the physical chemistry literature however a simple model is required that captures the essential effects with as little computational overhead as possible. For this purpose the same pairwise potential discussed by Im and Roux^[11] is implemented to include the effect of hydration shells.

${\displaystyle U_{hy}=c_{0}\exp \left({\frac {c_{1}-r}{c_{2}}}\right)cos(c_{3}(c_{1}-r)\pi )+c_{4}\left({\frac {c_{1}}{r}}\right)^{6}}$

The coefficients c_i were determined empirically for a 1 M KCl solution, using MD simulations to benchmark the ion radial distribution functions against Equilibrium Monte Carlo simulations. The effect of hydration shells was found to be important in simulations at higher salt concentrations where the conductance of many ion channels, porin among them, is observed to saturate as the salt concentration in the electrolyte baths is further increased. Earlier simulations that did not include a model of hydration shells did not reproduce the conductance saturation behavior. This suggests an additional repulsive potential acting to prevent ion crowding, and hence limiting the concentration of ions and current density in the confined space of the pore even at high bath salt concentration. When the repulsive potential was included moderate channel conductance was observed.

Boundary conditions

The electrical and physiological properties of ion channels are experimentally measured by inserting the channel into a lipid membrane separating two baths containing solutions of specific concentrations. A constant electrostatic bias is applied across the channel by immersing the electrodes in the two baths. Formulating boundary conditions that accurately represent these contact regions may require enormously large bath regions and is a challenging task. Beyond a Debye length from the membrane the electrostatic potential and ion densities do not vary appreciably. This assumption has been supported by the results of continuum results presented earlier.^[12] For typical salt concentrations used in ion channel simulations, the Debye length is of the order of 10 Å. Using the assumption, Dirichlet boundary conditions are imposed on the potential at the two domain boundary planes that are transverse to the channel, taking care that these planes are sufficiently far from the membrane.

The other problem in duplicating the experimental conditions is the problem of maintaining fixed charge density in the two baths. This problem is treated by maintaining the specified density in two buffer regions extending from the boundary plane toward the membrane. The number of ions needed to maintain the density in the two buffer regions is calculated at the start of the simulations. The count of the ions in these buffers is sampled throughout the simulation and an ion is injected whenever a deficit is observed. The initial velocity of the injected particle is decided according to Maxwellian distribution. The ions can leave the system only by exiting through the two Dirichlet boundary planes and an ion is not removed artificially from these buffer regions. The reflections from the Neumann boundary planes are treated as elastic reflections.

Multi-grids and grid focusing method

In all most any of the methods in simulation of ion channels, the major computational cost comes from the calculation of electrostatic forces acting on the ions. In continuum models, for instance, where ionic density exist rather than explicit ions, the electrostatic potential is calculated in a self-consistent manner by solving the Poisson equation. In MD simulations, on the other hand, the electrostatic forces acting on the particles are calculated by explicit evaluation of the Coulombic force term, often splitting the short-range and long-range electrostatic forces so they could be computed with different methods. In a model such as a reduced particle method, the longrange electrostatic forces are evaluated by solving the Poisson equation and augmenting the forces so obtained with a short-range component. By solving the Poisson equation it is possible to self-consistently include the forces arising from the bias to the system, while this is a difficult issue to be addressed in MD simulations.

Currently there are two Poisson solvers implemented in BioMOCA based on the finite difference method. One uses the pre-conditioned Conjugate Gradient scheme (pCG) and is used by default. The later is borrowed from an APBS solver, which uses a V-multi-grid scheme. Other than the numerical approach to solve the Poisson equation, the main difference between the two solvers is on how they address the permittivity in the system. In the first solver, a dielectric value is assigned to each cell in the grid, while in the APBS solver the dielectric coefficients are defined on the grid nodes. As discussed earlier box integration method is used in the pCG solver, which treats the Poisson equation in the most accurate way. Even though a full multigrid solver based on box-integration method has been under development, there is a neat way to reuse the already exiting code and treat the ion channel systems.

Ion channel simulations require the presence of large bath regions for accurate treatment of screening.^[1] There being of such bath regions make the mesh domain of Poisson equation large and leads to either a large number of grid points with fine mesh resolution or a small number of grid points with very coarse discretization. From bulk simulations a coarse mesh is sufficient for describing the baths using the P³M scheme. However, a fine resolution is required in the channel domain because of the highly charged nature of these regions and the presence of spatially varying dielectric regions. Besides the ultimate interest is to study the channel behavior in terms of ion permeability, selectivity, gating, density, etc.... In other words, it is better off to put more computational resources in the channel region, and bare minimum in the baths to reduce the overall computational cost and speed up of simulations from weeks to perhaps days instead. A scheme based on the grid focusing method has been developed that makes it possible to satisfy the requirement of large bath region and a fine grid resolution in channel at the same time in a computationally effective way. This methodology is capable to have multiple fine mesh domains, which may be needed to describe multiple pore channels like OmpF porin, or an array of ion channels sharing the same bath regions or even having yet finer meshes inside a fine mesh for relatively large channels with narrow ion passages like Nicotine receptor channel.^[13]

The first grid is coarse mesh spanning the entire problem domain including the bath regions and the channel region. The second grid (and so on for any other grids, 3rd, 4th, etc.) is a relatively much finer mesh that spans a sub-domain of the system containing the region that requires fine resolution like the channel pore. The Poisson equation is first solved on the coarse mesh with all the Dirichlet and Neumann boundary conditions, taking into account the applied bias. Next the boundary conditions for the secondary meshes are obtained by interpolating from the first or previous solutions of the Poisson equation. The Poisson equation is solved again for the finer meshes using the new boundary conditions. In this way, electrostatic fields with different mesh discretization for different regions can be generated.

EMF and DBF

The electro-motive-force (EMF) is the measurement of the energy needed for a charged particle like ion to cross the ion channel embedded in a membrane. Part of this potential energy barrier is due to the interaction between the crossing ion and the permanent/partial charges on the protein residues. The other part comes from the induced dipoles in the protein/membrane dielectric medium, and is referred as dielectric-boundary-force (DBF). To compute the DBF alone, one may turn off all the static charges on the protein residues and drag the ion through the pore and compute the energy barrier using

${\displaystyle P_{DBF}=\int -d{\hat {z}}.{\vec {E}}}$

It is important to note that EMF or DBF measurements are just qualitative measurements, as an ion does not necessarily cross the channel through the center of its lumen in a straight line and it is often accompanied by other ions moving in the same or opposite directions, which dramatically changes the dynamics of the system. Moreover, unlike steered MD calculations where the protein residues dynamically reposition themselves as an ion or ions are bouncing across the channel, in our EMF or DBF calculations protein is modeled as a static continuum, which further affects the energy calculations in a more quantitative way. Another issue that additionally impacts the measurements is absence of water hydration molecules, which move with the ion and shield part of its charge. Having said all of above, still computing EMF or DBF is valuable to address channel selectivity or gating. Computing either of these two energy barriers is available as an option in BioMOCA.

Visualization using VMD

VMD visualization of Gramicidin 1MAG molecule along with the structure generated by BioMOCA, where green represents protein, red addresses the membrane (i.e. lipid), and purple is the channel and left and right baths

VMD^[14] was equipped with the option of loading BioMOCA structures. This is a very useful feature as one could load both the protein structure (i.e. PDB or PQR file) along with the structures generated by BioMOCA to make comparisons. Figure at the right shows how BioMOCA has generated a structure for Gramicidin channel with a membrane wrapped around it. Furthermore, BioMOCA also dumps the ion trajectories in standard formats so they could be later loaded to molecular visualization tools such as VMD and watched frame by frame in a movie format.

Recording trajectories in binary

Other than counting the number of ions crossing the channel, sometimes it is desirable to study their behavior at different regions of the channel. Such examples would be the average occupancy of ions or their average moving velocity inside the channel or a nanopore. BioMOCA has been equipped with the option of dumping every ions position, average and instantaneous velocities, potential and kinetic energies, average and instantaneous displacements and other info at every step (or few steps) of the simulations in ASCII format, so such trajectory information could be studied later on to gather further statistics. From a technical point of view however, dumping such information for tens of ions, even at every few hundreds of time steps, could slow down the simulations and end up with huge files accumulating to tens of gigabytes. Loading such files later on from disk storage is also a very time-consuming and computationally inefficient procedure. Over and above that, recoding the numerical information in ASCII format does not hold its machine precision and has loss of accuracy.

Solving such problems is actually an easy task and it is simply to avoid using ASCII format and use binary format instead. Not only it preserves the machine accuracy but also writing and reading to file system is a lot faster. The computational overhead to dump the trajectories becomes negligible and the trajectory files become about two orders of magnitude smaller in size. The downside might be that programming and decoding the data could become very tricky, but once it's done correctly and with care, the advantages of using binary format are well worth the extra effort. BioMOCA is now equipped with the tools to record the trajectory information in binary format.