Biology Monte Carlo Suite

1. Biology Monte Carlo Suite An introduction to using the online BioMOCA Suite to run ion channel simulations. David Papke and Dr. Umberto Ravaioli Beckman Institute University of Illinois, Urbana-Champaign 3/30/2008

2. How does BioMOCA work? The BIoMOCA code was initially written by Trudi Van der Straaten, Gulzar Kathawala, and Dr. Umberto Ravaioli. NOTE: To learn more about the BioMOCA code and how it works, check out the following presentation by Reza Toghraee: • Toghraee, Reza (2007), "BioMOCA and Ion Channels", https://www.nanohub.org/resources/2751/, accessed on 2008-04-17 11:16:38. BioMOCA, or Biology Monte Carlo, is used to simulate ion flow through a channel. In these simulations, water and lipid are treated as a background dielectric materials, but ions and protein are treated discretely. Monte Carlo methods are used to determine when ion scatterings (from ion-water interactions) take place.

3. Overview • Step 0: Convert channel file format from .pdb to .pqr at http://pdb2pqr.sourceforge.net/ • Step 1: Use the map generator to load pqr structure and set up appropriate 3D simulation box. • Step 2: Use the lipid wrapper to set lipid boundaries and radius. • Step 3 (optional): Use the boundary force potential calculator to obtain the potential energy profile which faces a point charge traversing the channel. • Step 4: Run the Biology Monte Carlo (BioMOCA) Simulator to determine simulate ion flow through the channel.

4. Setting up the pqr file:How to obtain a pqr file from a pdb file • To begin, go to http://pdb2pqr.sourceforge.net/ and choose one of the three server mirrors. Note that currently the UCSD server provides the most lenient file size limits. • Enter all of the information requested by the server. For the force field field, choose CHARMM, and for the naming scheme, choose internal naming scheme (other naming schemes are sometimes incompatible with BioMOCA).

5. Setting up the pqr file:Aligning the pqr in VMD • Download and install VMD and the VMD orient script plug-in. These can be found at www.ks.uiuc.edu/Research/vmd • For BioMOCA, the pore of the channel must be aligned to the z-axis, with the intracellular opening pointing in the +z direction. Aligning the channel in this manner is important because in BioMOCA, positive ions flowing in the +z direction constitute positive currents (following electrophysiological convention). Extra-cellular end Intra-cellular end Here, the channel gramicidin is aligned with the intracellular domain pointing in the +z direction.

6. Setting up the pqr file:Rotating the channel • If your molecule is flipped (i.e. The intracellular end points in the –z direction), the following set of commands will flip it. set sel [atomselect 0 all] set A [transvecinv {0 0 -1}] $sel move $A set A [transaxis y -90] $sel move $A The TK Console can be accessed in the extensions menu. • These commands should be entered into the TK Console in VMD.

7. Setting up the pqr file:Orienting a channel • Occasionally, you will obtain a structure which isn’t aligned to any axis. To align it to the z-axis, use the following set of commands (note: you may still need to rotate the channel after aligning): package require Orient namespace import Orient::orient set sel [atomselect top "all"] set I [draw principalaxes $sel] set A [orient $sel [lindex $I 2] {0 0 1}] $sel move $A set I [draw principalaxes $sel] set A [orient $sel [lindex $I 1] {0 1 0}] $sel move $A set I [draw principalaxes $sel]

8. Setting up the pqr file:Aligning the pqr in VMD • After your pqr file has been properly aligned, open up the pqr file in a text editor to eliminate any headers or tails that the file may have. For our purposes, the first and last lines must looks as they do below An example of how the pqr file should look at its beginning and end

9. Using the map generator:Input • First, upload the pqr file in the pqr loader. Occasionally, you might get a forbidden access warning. Just try loading again if this happens, and your file should load within a few more attempts. • Next, set the Nx, Ny, Nz, dx, dy, and dz parameters to appropriate values. Note that Ni*di = box length along dimension i in Angstroms. It is important in the z direction to leave about 1/6 of the box on either end of the channel empty; this region will be occupied by bath ions during the simulation.

10. Using the map generator:Output • The visualization mode determines the number of x- and z-cross-sections you will see in the output. Mode 1 only gives two slices in each direction, but it renders much faster than mode 2, which gives all x- and z-cross-sections. Typically, mode 1 is sufficient; mode 2 is only really necessary for checking the lipid wrap. In addition to the cross-sections, the map generator also outputs a file called “Run Information,” which contains technical information about the program operation. The useful information for you are the warnings: if this files contains warnings about out of bounds atoms, it’s a sign that your box did not fully contain the protein. As a result, you need to define a bigger box before continuing to the next step. ~1/6 box length is free on both sides

11. Choosing an appropriate box size • To make sure your box is sufficiently large, check the “Run information” output file. If you see any “inside domain” warnings, than your box does not properly contain your protein.

12. Using the lipid wrapper After you have determined an appropriate 3D box size, the next step is to wrap your channel in lipid. In BioMOCA, the lipid is treated as a slab of material with a uniform dielectric constant. In the lipid wrapper, you enter the locations for the intra- and extra-cellular edges of this slab, as well as for the radius of the hole in which the protein exists. In general, the lipid should surround the protein as tightly as possible (i.e. without any dead space between the lipid and the protein) while never obstructing the pore. When picking an appropriate slab thickness, keep in mind that the plasma membrane of a cell is only around 30A thick. In the above picture, the protein (blue) has now been wrapped with lipid (green).

13. Boundary force potential calculator • In this sub-tool, which can only be executed after running the map generator and lipid wrapper, you can obtain the potential energy profile encountered by a point charge traversing your channel. • You can alter the dielectric constants to see how this affects the potential energy barrier. The potential energy profile of the gramicidin channel is shown above.

14. BioMOCA SimulatorOverview • Using BioMOCA, you can simulate particle flow through an ion channel. BioMOCA treats water and lipids as background dielectric materials, but treats ions and protein atoms as discrete particles. Every 100 time steps, the Poisson equation is solved for the entire system, and from it forces are calculated on each particle. Ion movement occurs every time step; to ensure that movements remain physical (i.e. no two ions overlap, or no ion overlaps with protein), close range calculations are performed every time step. • Prior to running BioMOCA, it is necessary to run the map generator and the lipid wrapper. However, it is not necessary to run the boundary force potential calculator.

15. BioMOCA Simulator:Inputs • Currently, the parameters you can adjust in BioMOCA are intra- and extra-cellular concentrations of Na+, Cl-, K+, Ca2+, and Mg2+, time step length, transmembrane potential, and simulation time.

16. BioMOCA Simulator:Output • For a small channel (< 1000 atoms), a 30ns simulation should take less than a day; for a large channel (>10,000 atoms), a 15-20ns simulation should take around 2 days. • Currently, the output of the BioMOCA Simulator consists of two files. The first file, “Run Information,” is a technical summary of program operation during the run.

17. BioMOCA Simulator:Output • All the useful information in this file pertinent to ion currents is also contained in the much easier to read “User Output” file. For each ion species, “User Output” relays the numbers of right to left and left to right crossings, as well as net numbers of crossings and currents. All of the information in this file conforms to electrophysiological conventions (for example, left to right, or inward to outward, flow of calcium ions constitutes a positive current, while a left to right flow of chloride ions would register negative).

18. Future plans In the future, we plan on adding ion concentration maps, potential distribution maps, and VMD movie files to the output of the BioMOCA Simulator. Additionally, we hope to give the user more options in BioMOCA inputs. Acknowledgements I’d like to thank Reza Toghraee and Dr. Umberto Ravaioli for their help in understanding and using BioMOCA.

19. Sources • H.M. Berman, J. Westbrook, Z. Feng, G. Gilliland, T.N. Bhat, H. Weissig, I.N. Shindyalov, P.E. Bourne: The Protein Data Bank. Nucleic Acids Research, 28 pp. 235-242 (2000). • Dolinsky TJ, Nielsen JE, MCCammon JA, Baker NA. PDB2PQR: an automated pipeline for the setup, execution, and analysis of Poisson-Boltzmann electrostatics calculations. Nucleic Acids Research, 32, W665-W667 (2004). • PDB2PQR software can be accessed at pbd2pqr.sourceforge.net. • Humphrey, W., Dalke, A. and Schulten, K., "VMD - Visual Molecular Dynamics", J. Molec. Graphics, 1996, vol. 14, pp. 33-38. • VMD and the orient script can be found at www.ks.uiuc.edu/Research/vmd. • Li H, Robertson AD, Jensen JH. Very Fast Empirical Prediction and Rationalization of Protein pKa Values. Proteins, 61, 704-721 (2005). • T.A. Van Der Straaten, G. Kathawala, A. Trellakis, R.S. Eisenberg, and U. Ravaioli, "BioMOCA - a Boltzmann transport Monte Carlo model for ion channel simulations," Molecular Simulatin, Vol. 31, No. 2-3, pp 151-171 (2005). David Papke email address: dpapke@uiuc.edu