This tutorial will walk you through two different types of BD simulations on the Thrombin-Thrombomodulin system, which is an important component of the blood-clotting cascade. Thrombin consists of 295 amino acids, and thrombomodulin consists of 117 amino acids. (My thanks to Adam Van Wynsberghe for the necessary data on these two molecules.)
This tutorial assumes that you have Linux on your computer. If you have only Windows, a quick and easy way to get Ubuntu Linux onto your machine is Wubi (Windows Ubuntu Installer). Other necessary software is the Ocaml compiler, which is free and can be obtained from here . If you are using Ubuntu or a similar Linux based on Debian, you can just type
sudo apt-get install ocamlYou will also need APBS; this can also quickly be installed on Ubuntu by typing
sudo apt-get install apbsFinally, you can install BrownDye on your machine by going to the website , downloading the file
browndye.tar.gz, moving the file to a directory of your choice, and typing
tar xvfz browndye.tar.gz cd browndye make allAll of the executables are now in the directory
browndye/bin; you need to add this directory to your path. If you use a bash shell, this can be done by typing
fullpathis the full path name of the directory containing the browndye distribution. If you now want to run the thrombin-thrombomodulin tutorial, go to
The atomic coordinates for thrombin and thrombomodulin are in
m.pqr. Because BrownDye works primarily with
XML files, you must convert this two files to an equivalent XML format:
pqr2xml < t.pqr > t-atoms.xml pqr2xml < m.pqr > m-atoms.xmlNext, you must generate the electrostatic grids, in dx format, using APBS:
apbs t.in apbs m.inwhere the input files
m.inare provided. The grids are output to
m-PE0.dx. Throughout this session, thrombin will be denoted by the prefix "t" while thrombomodulin will be denoted by prefix "m". Be sure to take note of the Debye length in the APBS output if you don't feel like calculating it by hand; it will be needed later.
In addition to atomic coordinates and grids, the third key input
is the set of reaction criteria. This can be generated from the
two coordinate files
m-atoms.pqrxml, and a file,
protein-protein-contacts.xml, which describes which pairs of atom
types can define a contact
(It is actually stored as
need to copy from this file.)
make_rxn_pairs takes these
three files and a search distance to generate a file of reaction
pairs. This assumes that the coordinates of the two molecules
are consistent with the bound state.
make_rxn_pairs -mol0 t-atoms.pqrxml -mol1 m-atoms.pqrxml -ctypes protein-protein-contacts.xml -dist 6.0 > t-m-pairs.xmlThe resulting file still is not suitable for input into the simulation programs, however. I've made the program general enough to have more than one reaction in a simulation, so one could envision having several such reaction pair files that would need to combined into a final reaction description file. For now, you can use the program
make_rxn_file, which generates an input file for the case of one reaction:
make_rxn_file -pairs t-m-pairs.xml -distance 5.5 -nneeded 3 > t-m-rxns.xmlThis generates a reaction description file which tells the simulation programs that if any 3 of the atom pairs approach within 5.5 Angstroms, a reaction occurs.
A note: if you don't feel like typing the above commands in, especially if you make changes and need to do it repeatedly, I have included a Makefile in the example directory. Typing
make allshould run the above commands.
The remaining pieces of information are contained in the file
(This must be copied from the file
It contains, among other things, information on the solvent,
information on each molecule, and parameters governing the simulation
For the sake of efficiency, the larger molecule
should be "molecule0" if there is a large size difference.
Also, each molecule is assigned a prefix as mentioned above; these
are used in naming the intermediate files that are generated in the
bd_top input.xmlThe bd_top program is written in Ocaml using a Unix "make"-like utility that I wrote to help orchestrate the creation of the files. Like "make", if an intermediate file is changed or replaced, running bd_top (which is analogous to a Makefile) will run only those commands necessary to re-generate files that depend on the updated file. Unlike "make", this utility, which I call the "Orchestrator", can also read information from xml files and have chains of dependent calculations (eventually I want to write a version in Python so it will look more familiar to most people). So, when the command is executed, the following files are generated for thrombin:
The following files are generated for both molecules:
bd_toprun again to update everything. For example, right now I'm using a simple test-charge approximation by default, but one could easily generate effective charges using another program such as SDA, convert the output into the appropriate XML format, replace
m-charges.xml, and run
Another useability note: if you want to clean things up, you
can delete all *.dx and *.xml files; the two xml files that you
need to get started are also available as
At this point, you can choose to do a simulation of one trajectory at a time, or you can do a weighted-ensemble simulation. In general, the weighted-ensemble method is not as efficient at the single-trajectory method unless the probability of a reaction event is very low.
To perform the single-trajectory simulation, the following is executed:
nam_simulation t-m-simulation.xmlThe results end up in
results.xml, as designated in
input.xml. As the simulation proceeds, you can look at
results.xmlto see the simulation progress. (This is called
nam_simulationafter Northrup, Allison, and McCammon, who came up with the first algorithm of this type. This code uses a fancy variation on the orginal algorithm.).
At any point, you can use
results.xml and obtain an estimate of and 95% confidence bounds
the reaction rate constant in units of M/s. The file
must also be given to this program:
compute_rate_constant < results.xmlThis will put the rate constant results to standard output.
To perform the weighted-ensemble method, you must first generate the bins for the system copies:
build_bins t-m-simulation.xmlThe number of system copies used in the bin-building process are given in
n-copiestag. As it runs, reaction coordinate numbers will go scrolling past; they should keep getting smaller and eventually stop. If that does not happen, i.e., the numbers keep on going, you might need to increase the number of system copies, or it might be that your reaction criterion is unattainable. Assuming is converges, the bin information is place in
t-m-bins.xml. The actual weighted-ensemble simulation is then run:
we_simulation t-m-simulation.xmlAs before, the results are output to
results.xml. In each row of output numbers, the right-most number is the flux of system copies that escaped without reaction, while the other ones are reactive fluxes. So, even for a rare reaction event, you should at least see small numbers for the reactive fluxes after the system has reacted steady-state. This can be visually examined at any point, and can also be analyzed as above, but with a different program:
compute_rate_constant_we -sim results.xml -solvent t-m-solvent.xmlBecause the streams of numbers are autocorrelated, a more sophisticated approach for computing confidence intervals is used, and if there are not enough data points, the program
compute_rate_constant_wewill simply refuse to provide an answer.
You can change the random number generator seed, under the
Good to do if you're bored but don't have the energy to do anything else.
One parameter to play with is the reaction criterion distance, which is the
-distance input to the program
make_rxn_file. The number given in
the tutorial and the Makefile is 5.5, but you can change that by re-running
make_rxn_file or by changing it in the Makefile.
You can also change the ionic strength in files
generate new APBS grids. Note: you must take note of the new Debye length
and put that value in the file
input.xml. So far, it is not possible to
automatically get the Debye length from the output DX file of APBS.
If your machine has several processors, you
can change the value under the tag
n-threads in file
see it run under several processors. So far, it runs only on shared-memory
Even on a single-processor
machine, you can still run several threads, but it does not make the
programs go any faster.
A final useability note: Most of the programs will output a description of themselves and their options if you type in
t.infound in the input files, generate the APBS grid
m.infound in the input files, generate the APBS grid
Browndye: pqr2xml, select the option
Input URL, using another browser tab go to the input files, and copy and paste the URL of
Input URLbox. Copy and pasting URL's can easily be done by moving the curser over the link, pressing the right mouse button and selecting
Copy Link Location.
t-atoms.pqrxmlin new browser tab, and do not delete the tab.
m.pqrto generate the file
m-atoms.pqrxmlin its own tab.
The following steps are used to generate the reaction pairs file:
Browndye: pqr2xmlservice into the box for
Molecule 0 Input URL.
Browndye: pqr2xmlservice into the box for
Molecule 1 Input URL.
protein-protein-contacts.xmlin the input files into the box for
Contacts Type URL.
t-m-rxn-pairs.xmlfile. Do not delete the tab.
The following steps are used to generate the reaction description file:
t-m-rxn-pairs.xmlfrom the output of
Pairs Input URLbox.
- Enter 5.5 into the
- Enter 3 into the
Number of Required Contactsbox.
- Submit the job, and inspect the results when finished.
The following steps are used to run
m-atoms.pqrxmlinto the boxes for
Molecule 0 URLand
Molecule 1 URL.
Reaction File URLbox.
m.dxinto the APBS boxes for Molecule 0 and Molecule 1.
input.xmlin the input files into the
Input File URLbox.
results.xmlfinally appears in the
Output Base URL, keep checking it if you want until 1000 trajectories have been run and the job stops.
results.xmlinto the input URL box for the
Browndye: compute_rate_constantservice, and check the file
stdout.txtfor the rate constant information.
Huber, GA and McCammon, JA. Browndye: A Software Package for Brownian Dynamics. Computer Physics Communications 181, 1896-1905 (2010)
Ermak, DL and McCammon, JA. Brownian Dynamics with Hydrodynamic Interactions, J. Chem. Phys. 69, 1352-1360 (1978)
Northrup SH, Allison SA and McCammon JA. Brownian Dynamics Simulation of Diffusion-Influenced Bimolecular Reactions J. Chem. Phys. 80, 1517-1526
Luty BA, McCammon JA and Zhou HX. Diffusive Reaction-Rates From Brownian Dynamics Simulations - Replacing the Outer Cutoff Surface by an Analytical Treatment, J. Chem. Phys. 97, 5682-5686 (1992)
Huber GA and Kim S. Weighted-ensemble Brownian dynamics simulations for protein association reactions, Biophys. J 70, 97-110 (1996)
Gabdoulline RR and Wade RC. Effective charges for macromolecules in solvent, J. Phys. Chem 100, 3868-3878 (1996)