.:: MCCE ::.
 
.:: Database ::.
 
 
.:: About Us ::.
 
 
.::  Login ::.
 
 
.:: Documentation ::.
 

2.2. Submitting and monitoring mcce job, a simple example

Here is an example of calculating pKa values for Hen egg lysozyme. In this example, mcce was installed under /usr/local/mcce-1.0/, and the directory /usr/local/mcce-1.0/bin/ is in my excutable searching path. First make a directory for this set of run. You get the structure file 4LZT from the rscb data site with command "getpdb" included in MCCE package:

jmao@metis jmao]$ mkdir 4lzt
[jmao@metis jmao]$ cd 4lzt
[jmao@metis 4lzt]$ getpdb 4lzt
Inquiring the remote file 4LZT.pdb ...
Saving as 4LZT.pdb ...
Download completed.
[jmao@metis 4lzt]$

Then copy a "run.prm" file to this directory and edit it:

[jmao@metis 4lzt]$ cp /usr/local/mcce/run.prm.quick run.prm

You will need to edit this "run.prm" file. In the edtitor, make first string of line (INPDB) "4LZT.pdb", and turn flags for all 4 steps to "t". You may want to run one step at a time and look at the output and run.log messages for each step.

The Most modified section will look like this:

======================================================================
Most modified entries
----------------------------------------------------------------------
Input and Output:
4LZT.pdb (INPDB)

Steps:
t step 1: pre-run, pdb->mcce pdb (DO_PREMCCE)
t step 2: make rotatmers (DO_ROTAMERS)
t step 3: do energy calculations (DO_ENERGY)
t step 4: monte carlo sampling (DO_MONTE)
======================================================================

MCCE can be invoked simply by running the executable file "mcce". As the program will take hours to days to finish, you may want to run the program at the background and redirect the output message to a log file using command "mcce > run.log &".

[jmao@metis 4lzt]$ mcce > run.log&
[1]8533
While running, MCCE writes out the progress, warning messages and error messages directed here into file "run.log". Check mcce run time messages and serious errors will appear in it. The "debug.log" and "progress.log" provide extra warning messages and progress report, which can be safely ignored most time.

In running 4LZT I noticed a line in the run.log file:

Error! premcce_confname(): Can't get conformer list of this residue NO3

This means a cofactor NO3 was not recognized by mcce. The program continues anyway and treats all 4 atoms of NO3 with 0 charge and radii of 1.7 Angstroms. You may want to stop the program and make a TPL (residue paramter file) file, but we will ignore the error. If NO3 is on the surface, it won't cause much error in pKa becuase of 1) the screening effect of the solvent and 2) Poisson-Boltzmann equation will take care of the salt effect if this surface NO3 is stripped off like surface water by the mcce program.

After step 1, the file "step1_out.pdb" is a mcce format pbd file with terminal residues split off to form an independently titratable NTR and CTR, groups renamed based in file "name.txt", and surface water and salt molecules with 10% solvent accessibility stripped off. Undefined cofactors will be identified as made up of non-charged atoms. A temporary parameter file "new.tpl" will be written out. This can help you star to make a tpl file for new groups. The file "acc.res" lists the solvent accessibility of residues and the file "acc.atm" lists the solvent accessibility of atoms. No conformers are made for groups with more than 30% surface accessibility in the configuration of "run.prm.default".

After step 2, you will get a multi-conformation pdb file "step2_out.pdb", which is the only file used by step 3. The rotamer making statistics are dynamically updated in file "rot_stat".

The DelPhi runs in step 3 are the most time consuming step in mcce. If you look at the bottom of run.log when DelPhi is running it will look like this:

[jmao@metis 4lzt]$ tail run.log

Creating connectivity table...
Done

Preparing DelPhi runs ...
3 focusing runs required for this protein.
Done

Computing pairwise from conformer 1 to 281 of 281 total conformers
see progress.log for progress...

So, we take a look at the additional progress file "progress.log":

[jmao@metis 4lzt]$ tail progress.log			
Doing pairwise 72 of 281 conformers. 5 seconds
Doing pairwise 73 of 281 conformers. 2 seconds
Doing pairwise 74 of 281 conformers. 4 seconds
Doing pairwise 75 of 281 conformers. 2 seconds
Doing pairwise 76 of 281 conformers. 3 seconds
Doing pairwise 77 of 281 conformers. 3 seconds
Doing pairwise 78 of 281 conformers. 4 seconds
Doing pairwise 79 of 281 conformers. 4 seconds
Doing pairwise 80 of 281 conformers. 5 seconds
Doing pairwise 81 of 281 conformers.[jmao@4lzt]$

This progress is reported in "progress.log" and it can be used to estimate how long step 3 will take. A pairwise calculation and a reaction field energy calculation will be performed for each conformer. So step 3 will take about 4 seconds*281 conformers*2 DelPhi's= 2248 seconds or 37 minutes. The output directory of step 3 is "energies" and the output file is "head3.lst".Both are needed as input of step 4.

Step 4 is Monte Carlo sampling, which simulates a pH or Eh titration. It writes the progress to file "mc_out". The occupancy of each conformer at each pH is in fort.38. The net charge of each ionizable residue is in "sum_crg.out". The pKa or Ems are in "pK.out". These are obtained by non-linear fitting of the titration curves in file "sum_crg.out".

One reason to run simulations is to be able to ask WHY. To understand what contributes to the ionization energy of a residue at a specific pH, use the tool "mfe.py". This program approximates the interactions between residues by Mean Field approach. To look at the free energy of the reaction GLU -> GLU- + H+ for Glu 58 at pH 5.078, the midpoint of the titration curve of this residue run.

[jmao@metis 4lzt]$ mfe.py GLU-_0035_ 5.078 0.2
Residue GLU-_0035_ pKa/Em=5.078
=================================
Terms pH meV Kcal
---------------------------------
vdw0 0.03 1.98 0.05
vdw1 0.08 4.91 0.12
tors -0.00 -0.00 -0.00
ebkb -1.39 -80.41 -1.89
dsol 2.72 157.67 3.71
offset -0.99 -57.45 -1.35
pH&pK0 -0.33 -19.04 -0.45
Eh&Em0 0.00 0.00 0.00
-TS 0.11 6.35 0.15
residues -0.15 -8.47 -0.20
*********************************
TOTAL 0.10 5.54 0.13
*********************************
ASP_0048_ 0.22 13.02 0.31
ASP_0052_ 0.95 54.90 1.29
ARG_0112_ -0.29 -16.64 -0.39
ARG_0114_ -0.24 -13.93 -0.33
=================================
The program mfe.py is bundled with the mcce package. The syntax of the first argument about how to refer to residue should be taken from pK.out. The second argument means the Boltzmann distribution of conformers for all residues at 5.078 will be used. The Last argument says only specific interactions with individual residues over 0.2 kcal/mol will be printed out individually. The result explains the free energy of the reaction GLU -> GLU+ at pH 5.078, the midpoint of the titration curve of this residue, where reactant and product are equally favored. So it is expected here that the ionization energy in line TOTAL is close to 0. The small error is caused by mean field approximation and this error may indicate the degree of coupled ionization or conformational changes with other residues. The breakdown of ionization energy helps identifying the factors controlling the pKa of the residue. Negative number favors ionization, so ASP52 in this protein shifts up the pKa of GLU35 by 0.95 pH unit.

Listed in the table are other energy terms. All these terms are the difference between the contributions from the disignated term to the ionized conformers and netural conformers. "vdw0" is van der Waals interaction of the side chain of this residue. "vdw1" is the van der Waals interaction of the side chain of this residue to the backbone atoms. tors is the torsion energy of this residue. "ebkb" is the electrostatic interaction of the side chain to the backbone. "dsolv" is the loss of solvation energy. "offset" is a correction for each residue type defined in "extra.tpl". "pH&pK0" is the pH effect, which equals solution pH minus solution pKa of that residue if the residue is an acid or solution pKa of that residue minus solution pH if the residue is a base. "Eh&Em0" is the redox potential effect, which is the current solution potential minus the standard redox potential of that cofactor. "-TS" is entropy effect. As the protonated and deprotonated residue have have different constraints in a protein, the entropy of the two ionization states could be different. "residues" refers to the mean field interaction with residues given their conformations at this pH. This is the term that is pH or Eh dependent.

A mcce web interface analysis tool is provided separately. If it is installed, just point your browser to "http://localhost/cgi-bin/mcce.py" (replace localhost by the host name of mcce server if you run mcce on a remote computer). For more information, refer to section 4, Understanding MCCE Results.

<Prev | Contents | Next >

 
 

:: Doc Index

 

:: Complete Doc