Biology Monte Carlo methods (BioMOCA) have been developed at the University of Illinois at Urbana-Champaign to simulate ion transport in an electrolyte environment through ion channels or nano-pores embedded in membranes. It is a 3-D particle-based Monte Carlo simulator for analyzing and studying the ion transport problem in ion channel systems or similar nanopores in wet/biological environments. The system simulated consists of a protein forming an ion channel (or an artificial nanopores like a Carbon Nano Tube, CNT), with a membrane (i.e. lipid bilayer) that separates two ion baths on either side. BioMOCA is based on two methodologies, namely the Boltzmann transport Monte Carlo (BTMC)C. Jacoboni, P. Lugli, The Monte Carlo Method for Semiconductor Device Simulation, Springer Verlag, New York (1989) and particle-particle-particle-mesh (P3M).R. Hockney, J. Eastwood, Computer Simulation Using Particles, McGraw-Hill, New York (1981) The first one uses Monte Carlo method to solve the Boltzmann equation, while the later splits the electrostatic forces into short-range and long-range components.
Backgrounds In full-atomic molecular dynamics simulations of ion channels, most of the computational cost is for following the trajectory of water molecules in the system. However, in BioMOCA the water is treated as a continuum dielectric background media. In addition to that, the protein atoms of the ion channel are also modeled as static point charges embedded in a finite volume with a given dielectric coefficient. So is the lipid membrane, which is treated as a static dielectric region inaccessible to ions. In fact the only non-static particles in the system are ions. Their motion is assumed classical, interacting with other ions through electrostatic interactions and pairwise Lennard-Jones potential. They also interact with the water background media, which is modeled using a scattering mechanism.
The ensemble of ions in the simulation region, are propagated synchronously in time and 3D space by integrating the equations of motion using the second-order accurate leap-frog scheme. Ion positions r and forces F are defined at time steps t, and t + dt. The ion velocities are defined at t – dt/2, t + dt/2. The governing finite difference equations of motion are
where F is the sum of electrostatic and pairwise ion-ion interaction forces.
Electrostatic field solution The electrostatic potential is computed at regular time intervals by solving the Poisson’s equation
where and are the charge density of ions and permanent charges on the protein, respectively. is the local dielectric constant or permittivity, and 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. Solving the Poisson equation on the grid counts for the particlemesh component of the P3M 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. 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. 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, 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,S. Selberherr, Analysis and Simulation of Semiconductor Devices, New York, Springer-Verlag Wien, (1984). which has been briefly described below.
Calculations
Box integration discretization In order to use box integration for discretizing a D-dimensional Poisson equation
with being a diagonal D × D tensor, this differential equation is reformulated as an integral equation. Integration the above equation over a D-dimensional region , and using Gauss theorem, then the integral formulation is obtained
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. is assumed to be given on the rectangular regions between nodes, while is defined on the grid nodes (as illustrated on figure at the right).
The integration regions are then chosen as rectangles centered around node and extending to the 4 nearest neighbor nodes. The gradient is then approximated using centered difference normal to the boundary of the integration region , and average over the integration surface . This approach allows us to approximate the left hand side of the Poisson equation above in first order as
where and are the two components of the diagonal of the tensor . Discretizing the right-hand side of the Poisson equation is fairly simple. is discretized on the same grid nodes, as it's been done for .
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
Here, is the Lennard-Jones energy parameter and 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 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. In between these scattering events, ions follow the Newtonian forces. The free flight times, Tf, are generated statistically from the total scattering rate according to