Biology Monte Carlo method

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

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.

It has become evident that the

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}}$

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.

Availability of high-resolution X-ray crystallographic measurements of complete

Ideally, the steric interactions between protein atoms and the ions in the aqueous medium are to use a repulsive potential like

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.

It is widely accepted that the ions and water molecules do not have the same mobility or diffusivity in confined regions as in bulk.

${\displaystyle \lambda }$

${\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.

In addition to having a diffusive effect on

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

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

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.

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

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

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.

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.

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

• Monte Carlo method
• Biology
• Computational biology