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) and particle-particle-particle-mesh (P³M). 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.

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 3-D 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.

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

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. 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.

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.