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.) Another part of the tutorial deals with a bifunctional enzyme and substrate channeling.
BrownDye can be easily installed on Linux or Mac OS X. If you have Windows, you can install Cygwin , which provides a free Linux emulator, and then to install BrownDye using the Cygwin terminal. Go here for more details. Other software required will be APBS , Ocaml, and VMD for the visualization.
sudo apt-get install ocamlYou will also need APBS; this can also quickly be installed on Ubuntu by typing
sudo apt-get install apbsThe following steps will also work on the Mac terminal and Windows Cygwin terminal. 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
PATH=$PATH:fullpath/browndye/binwhere
fullpath
is the full path name of the directory containing the
browndye distribution.
If you now want to run the thrombin-thrombomodulin tutorial, go to
browndye/thrombin-example
.
The atomic coordinates for thrombin and thrombomodulin are in
files t.pqr
and 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
t.in
and m.in
are provided. The grids are
output to t.dx
and m.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. (Some older versions of APBS will name the output files
t-PE0.dx
and m-PE0.dx
; if that happens, be sure to
rename the files.)
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 t-atoms.xml
and m-atoms.xml
, and a file,
protein-protein-contacts.xml
, which describes which pairs of atom
types can define a contact
(It is actually stored as protein-protein-contacts.xml.bak
; you
need to copy from this file.)
The program 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.xml -mol1 m-atoms.xml -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 4.9 -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 4.9 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 input.xml
.
(This must be copied from the file input.xml.bak
.)
It contains, among other things, information on the solvent,
information on each molecule, and parameters governing the simulation
itself.
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
next step:
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_top
run 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 t-charges.xml
and m-charges.xml
, and run bd_top
again.
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 input.xml.bak
and
protein-protein-contacts.xml.bak
.
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.xml
to
see the simulation progress. (This is called nam_simulation
after
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 compute_rate_constant
to
analyze results.xml
and obtain an estimate of and 95% confidence bounds
on
the reaction rate constant in units of M/s. The file t-m-solvent.xml
must also be given to this program:
compute_rate_constant < results.xmlThis will put the rate constant results to standard output.
To visualize a trajectory leading to a reaction, you must include the line
<trajectory-file> trajectory </trajectory-file>in the input file, run
bd_top
and nam_simulation
again. The value trajectory
gives its names to the resulting trajectory
files, and can be what you want.
Because the files can get quite large, you might want
to also include a line such as
< n-steps-per-output > 10 < /n-steps-per-output >which causes only every 10th configuration to be output. At the end, the files
trajectory0.xml
and
trajectory0.index.xml
will be generated. The first file contains the actual trajectories in a
compressed format, and the second file functions as an index to the
first file.
If you are running with more than one thread, more files
like this will also be generated, with names trajectory1.xm1
,
trajectory1.index.xml
, etc., with one pair of files per thread.
After the initial trajectory files are generated, they must be processed further. For example, if you wish to view a trajectory that results in a reaction, you
must find such a trajectory by running
process_trajectories -traj trajectory0.xml -index trajectory0.index.xml -srxn associationwhich prints out the numbers of reactive trajectories. Then, pick one of the numbers, say 8, and run
process_trajectories -traj trajectory0.xml -index trajectory0.index.xml -n 8 > trajectory.xmlThis uncompresses the contents of
trajectory0.xml
and outputs
the desired trajectory into trajectory.xml
.
At this point, you can cut down even further on the number of configurations by
adding to the above to get
process_trajectories -traj trajectory0.xml -index trajectory0.index.xml -n 8 -nstride 10 > trajectory.xmlwhich outputs only every tenth configuration in
trajectory0.xml
.
The final step is to generate an xyz-format file for visualization by
VMD, which is obtained by running
xyz_trajectory -mol0 t-atoms.xml -mol1 m-atoms.xml -trajf trajectory.xml > trajectory.xyzThis file can be very large, since all of the atom positions for every configuration are output to
traj.xyz
, so it can be important to
reduce the number of trajectories upstream. This file then can be passed
to VMD for visualization. If you are worried about the size of the file,
you can run
xyz_trajectory -mol0 t-atoms.xml -mol1 m-atoms.xml -trajf trajectory.xml -trialwhich will output nothing except an estimate of the file size. If it looks like it might be too large, you can go back and increase the stride argument to
process_trajectories
.
At last, you can start up VMD, load in the xyz file, and watch the animation.
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
input.xml
in the n-copies
tag.
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 it converges, the bin information is placed 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_we
will
simply refuse to provide an answer.
5,10-methylenetetrahydrofolate + deoxyuridine monophosphate ↔ dihydrofolate + thymidine monophosphate
The dihydrofolate is then reduced by the DHFR. In some organisms, the two enzymes are joined together, and there is some evidence that a positivly charged channel helps convey the negatively charged dihydrofolate from one active site to another. A Brownian dynamics simulation by Elcock et al. some years ago on enzyme complex from the protozoan that causes a nasty tropical disease, leishmaniasis, showed a very pronounced channeling effect. We unfortunately cannot reproduce that here since the coordinates of that enzyme complex are not publicly available, but I have set up a simple system using the complex from the agent of African sleeping sickness, another nasty disease.
The enzyme complex is 3QGT from the protein data bank. It is missing many residues in the region between the two different enzymes, so we may not have a very complete channel in this model. Furthermore, work would need to be done to find the binding mode of the dihydrofolate to the DHFR active site. So, for now, we will need to settle for a small sphere with a -2 charge to assess the hypothesis of channeling. The simulation starts the "substrate" in the TS active site. The complex itself is a dimer, with two TS and two DHFR units. The subtrate is advanced until it either reacts or escapes. You can explore the effect of electrostatic channeling by changing the charge on the substrate sphere.
You will find that the effect here is noticeable, but not dramatic; this may because of the incompleteness of the model or actual reality. There are other structures of the DHFR-TS complex in the protein data bank from medically interesting organisms. Perhaps some refining of the model along with upcoming improvements in Browndye such as flexible molecules will lead to some interesting results.
The data for the channeling simulation is found in the channeling-example
directory in the Browndye distribution.
You can change the random number generator seed, under the seed
tag.
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 4.9, 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 t.in
and m.in
and
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 input.xml
and
see it run under several processors. So far, it runs only on shared-memory
machines using pthreads
.
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
program -help
t.pqr
and t.in
found in the
input files,
generate the APBS grid t.dx
.
m.pqr
and m.in
found in the input files,
generate the APBS grid t.dx
.
Browndye: pqr2xml
, select the option Input URL
,
using another browser tab go to the input files, and copy and paste the URL of t.pqr
into the Input URL
box. 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
.
Submit
button
t-atoms.xml
in new browser tab, and do not delete the tab.
m.pqr
to generate the file m-atoms.xml
in its own tab.
The following steps are used to generate the reaction pairs file:
Browndye: make_rxn_pairs
.
t-atoms.xml
from Browndye: pqr2xml
service into the box for Molecule 0 Input URL
.
m-atoms.xml
from Browndye: pqr2xml
service into the box for Molecule 1 Input URL
.
protein-protein-contacts.xml
in the input files into the
box for Contacts Type URL
.
Search Distance
box
t-m-rxn-pairs.xml
into the Output File
box.
t-m-rxn-pairs.xml
file.
Do not delete the tab.
The following steps are used to generate the reaction description file:
Browndye: make_rxn_file
.
t-m-rxn-pairs.xml
from the output of
make_rxn_pairs
to the Pairs Input URL
box.
Reaction Distance
box.
Number of Required Contacts
box.
t-m-rxns.xml
into the Output File
box.
The following steps are used to run bd_top
and nam_simulation
bd_top
.
t-atoms.xml
and m-atoms.xml
into
the boxes for Molecule 0 URL
and Molecule 1 URL
.
t-m-rxns.xml
into the Reaction File URL
box.
t.dx
and m.dx
into the APBS boxes
for Molecule 0 and Molecule 1.
input.xml
in the input files into the
Input File URL
box.
results.xml
finally appears in the Output Base URL
,
keep checking it if you want until 1000 trajectories have been run and
the job stops.
results.xml
into the input URL box
for the Browndye: compute_rate_constant
service, and check the file
stdout.txt
for the rate constant information.
setup.exe
file. The only part that requires some
effort is finding the necessary packages to install. If you do nothing,
Cygwin installs a minimal subset, but to run Browndye you need gcc,
make, autoconf, and ocaml. You could install all optional packages, but
it makes for a huge and lengthy download. Once Cygwin is installed,
you can pop open a terminal window, giving you access to a directory
tree just like any other Unix system. If you are missing packages,
you can run setup.exe
again, which will again present you
the choice of packages.
Another option for Windows is to install Ubuntu Linux using Wubi. From Windows, you download and run the Wubi installer, which in turn downloads and installs Ubuntu Linux and the boot-up tools necessary to boot into either Windows or Ubuntu. No disk partioning is necessary, making it very easy. The download usually takes several minutes depending on the speed of your connection. Then, when you restart your computer, you are given a choice between booting into Windows or Ubuntu. Afterwards, if you want, Ubuntu is easily deleted by running the Wubi Uninstall program from within Windows.
In order to obtain gcc for OS X on the Mac, you need to go to Apple's website and register as a developer. If you are doing your work in a university, this is free. Then, you install the XCode Developer tools, which include gcc. Ocaml is available for the Mac from the Ocaml website.
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)