Printer Friendly - AutobusterLigandQMtutorialIntro2i0g

Direct use of weighted Quantum Chemical Energy for ligands in BUSTER refinement

QM introductory tutorial example:

benzopyran ligand bound to `2i0g`

Contents

Background
Starting files
(A) Preparation: add hydrogen atoms to the ligand(s)
(B) Refine with -qm I0G
(C) Control refine with dictionary
(D) Analysis comparing the QM and dictionary refine jobs
(E) Assessing ligand strain energy and motions for the QM result

Background

For an introduction to the ideas behind the weighted-QM energy as a restraint function for ligands see AutobusterLigandQMIntroduction.
Here we will use the method to treat the benzopyran ligand I0G in the estrogen receptor beta structure 2i0g.

image from ligand expo: http://ligand-expo.rcsb.org/reports/I/I0G.

2i0g is a structure of the estrogen receptor beta with the I0G ligand bound.

The structure was solved by Yong Wang (Eli Lilly) and diffracts to 2.5Å resolution. It has respectable geometry as gauged by MolProbity and its reciprocal space correlation coefficient plot looks very reasonable.
The receptor has 2-fold NCS with a ligand bound in both sites. The density for the ligands is clear for the resolution.

This tutorial shows that weighted-QM energy ligand restraints yield a conformation similar to that obtained from a GRADE dictionary. In addition the QM method indicates that the ligand is slightly strained on binding, with the phenolic ring being rotated from the unbound state.

Starting Files

2i0g.pdb	protein/ligand coordinates from the RCSB
2i0g.mtz	structure factors
I0G.grade_PDB_ligand.cif	`grade_PDB_ligand` dictionary

The PDB and MTZ files were obtained by using collect.pl, see AutoBusterExampleUsingPDB.
The GRADE dictionary for the ligand I0G, based on Mogul CSD information, was generated by:

grade_PDB_ligand I0G

See the GRADE wiki page for details.

(A) Preparation: add hydrogen atoms to the ligand(s)

To run the QM method it is essential that your ligand is chemically complete, in particular its hydrogen atoms must be present. In 2i0g.pdb the two copies of I0G do not have hydrogen atoms modelled so we must first add them.
To do this run the hydrogenate tool:

hydrogenate -p 2i0g.pdb \
           -l I0G.grade_PDB_ligand.cif \
           -ligonly -zero \
           -o 2i0g_addhI0G.pdb

hydrogenate

Uses Richardsons' Lab reduce program so this must be on your PATH. See http://kinemage.biochem.duke.edu/software/reduce.php for information (it is also included in MolProbity).
hydrogenate understands the hydrogen atoms specified in a Refmac/CCP4-style CIF restraint dictionary.
The -ligonly option means add hydrogen atoms only to the ligand and not the protein.
The -zero option means assign zero occupancy to the hydrogen atoms added. This is recommended because at this resolution explicit full occupancy hydrogen atoms are not recommended.
The -o option specifies the output filename. Use COOT to check that the hydrogen atoms have been added correctly:

coot -p 2i0g_addhI0G.pdb

You should see that both the A and B copies of the I0G ligand have hydrogen atoms added:

(B) Refine with `-qm I0G -l I0G.grade_PDB_ligand.cif`

First add hydrogen atoms to the ligand.
The refine option -qm I0G will use the QM method as restraints for every copy of I0G in the protein:

time refine                 \
-p 2i0g_addhI0G.pdb         \
-m 2i0g.mtz                 \
-autoncs                    \
-l I0G.grade_PDB_ligand.cif \
-qm I0G                     \
-d partB > partB.log &

As 2i0g has 2-fold NCS the -autoncs option is used. LSSR restraints will be used for both the protein and the I0G ligand.
Note that a CIF restraint dictionary must still be specified for the ligand. Although the QM method uses zero weights, the conventional stereochemical restraints between I0G atoms BUSTER still have to know the correct atom types to work out ideal-distance contacts and to set up normal BCORREL restraints for I0G.
To check that QM restraints are being used, use firefox (or other browser) to examine the initial file:

firefox partB/01-BUSTER/Cycle-1/LIST.html

Scroll down until you see a section:

first QM calculation of energy/gradient
.....
    QM_HELPER_LOG: Hamiltonian            =            RM1  Number of QM Atoms     =             39
    QM_HELPER_LOG: System Charge          =              0  System Multiplicity    =              1
    QM_HELPER_LOG: Number of Orbitals     =            102  Occupied Orbitals      =             54
    QM_HELPER_LOG: Number Alpha Electrons =             54  Number Beta Electrons  =             54
    QM_HELPER_LOG: Number of Shells       =             60  Number of Primitives   =            360
    QM_HELPER_LOG: Number of 2e Integrals =          24933  Energy Base Line       =     333346.698
    QM_HELPER_LOG: No. of Boundary Atoms  =              0  Matrix Size            =           5253

This shows that the RM1 semi-empirical QM method is used to assess the energy of the chain A copy of I0G. The logfile shows that I0G has 39 atoms and that 102 orbitals are used, with further details of the QM calculation.

Note that as there is a second copy of I0G in the B chain, a separate RM1 calculation is used to assess its energy.

While the refinement is running you can check how the QM energy for the two copies of I0G changes by:

graph_autobuster_QM -d partB -launch

This produces a plotmtv graph which looks like (here the run has only completed 60 small cycles):
The QM energy plotted is relative to the initial energy for the first set. It can be seen that during refinement the RM1 energy for both ligands rapidly drops and then plateaus. What does this mean? Basically because in the PDB deposition the IOG ligand has been refined using CNS restraints, the ideal geometry is slightly different. All that is happening is that the ligand is adjusting its conformation to that preferred by RM1. The adjustments are small.

While the job is running, graph_autobuster_R can be used to monitor the progress of Rwork and Rfree:

graph_autobuster_R -d partB -launch

In addition, after one big cycle is completed the refined structure and its maps can be quickly examined with:

visualise-geometry-coot partB

See the BusterOutputInterpretation page for more details.

Once partB has completed examine the IOG ligand geometry.

Start COOT by:

visualise-geometry-COOT partB

In COOT also load 2i0g.pdb for comparison.
Use COOT "Go to atom" to find IOG A 1.
Here the QM-refined ligand is shown in orange, and compared with the original PDB conformation in green. The 2Fo - Fc density is contoured at 1.3 RMSD and the difference density at 3 RMSD.
In terms of fitting the macromolecular data the changes are small and the fit to the map is the same.

However, if you use Mogul to examine the bond lengths and angles for the QM-refined I0G in comparison to the original I0G from the PDB file, many fewer deviations are found (see analysis below).
Looking at the Rwork/Rfree graph for the completed partB job, the refinement produces an nice drop in Rfree by about 1.2%. Is this due to the use of QM for the ligand? Of course it isn't! - the difference is the result of running a long refinement with BUSTER, particularly with -autoncs. We will show this by running a control refinement in part (C).

Advanced users You will notice that both Rwork and Rfree are still dropping at the end of the partB refinement. We also have not applied TLS. To get the most out of the model before additional building it is appropriate to do a further autoBUSTER refine run. Notice that this still uses -qm I0G but this time with -M TLSbasic as well as -autoncs.


  time refine                 \
  -p partB/refine.pdb         \
  -m 2i0g.mtz                 \
  -autoncs                    \
  -l I0G.grade_PDB_ligand.cif \
  -qm I0G                     \
  -M TLSbasic                 \
  -d partB2 > partB2.log      &

The further partB2 refinement with TLS results in a drop of Rwork/Rfree from 0.2153/0.2517 to 0.2053/0.2420 and further improvement in MolProbity scores. To keep things simple we will use the original run.

(C) Control refine with dictionary `-l I0G.grade_PDB_ligand.cif`

As a control do the same refine job but leave out the -qm so I0G geometry will be handled by the GRADE dictionary:

time refine                 \
-p 2i0g_addhI0G.pdb         \
-m 2i0g.mtz                 \
-autoncs                    \
-l I0G.grade_PDB_ligand.cif \
-d partC > partC.log &

(D) Analysis comparing the QM and dictionary `refine` jobs

Let's compare results, first concentrating on the overall refinement metrics:

What	original `2i0g.pdb`	partB `-qm`	partC GRADE dictionary
CPU secs (1)	-	1684	1638
Rwork (2)	0.2329	0.2157	0.2173
Rfree (2)	0.2677	0.2524	0.2512
MolProbity remark 40 clashscore	9.96	0.28	0.28
MolProbity remark 40 Ramachandran favoured	95.1%	97.9%	98.4%
RMS bond deviation in Å	0.0117	0.0092	0.0090
RMS angle deviation in degs.	1.35	1.06	1.04
Final X-ray weight	-	6.90	6.52

refine.log reported CPU time on a Intel(R) Core(TM)2 CPU 6700 @ 2.66GHz (4 year old basic workstation)
BUSTER values

Let's now look at the I0G ligand, concentrating on the A site (because the B site is really pretty much identical). First let's extract the A site ligand from the various coordinate files (together with the CRYST1 card so the resultant PDB file works with crystallographic tools).

egrep "CRYST1|I0G A" 2i0g.pdb > 2i0g.I0G_A.pdb
egrep "CRYST1|I0G A" partB/refine.pdb > partB.I0G_A.pdb
egrep "CRYST1|I0G A" partC/refine.pdb > partC.I0G_A.pdb

We can assess how well the different ligand conformations fit to the density by eye using COOT with the result that they appear very similar and fit the density about equally well. A correlation coefficient to the BUSTER 2Fo-Fc map can be found by running:

corr -p partC.I0G_A.pdb -m partB/refine.mtz

These confirm that the density fits are equally good.

What	original `2i0g.pdb`	partB `-qm`	partC GRADE dictionary
I0G ligand CC (to partB/refine.mtz)	0.937	0.938	0.937

What about the geometry of the ligand? This can be assessed by the CCDC Mogul program, compared to similar CSD structures. We have found that a useful metric is obtained by running Mogul in default mode and searching for "All fragments", "All bond fragments", sorting by Z and finding the number of bonds with Z above 2, as well as the average Z score (using a spreadsheet such as OpenOffice or gnumeric).

What	original `2i0g.pdb`	partB `-qm`	partC GRADE dictionary
I0G ligand site A `Mogul` number of bonds with Z score > 2.0	5/24	0/24	0/24
I0G ligand site A `Mogul` average Z score for bonds	1.3	0.49	0.19

So both the refinement using QM (RM1 method) for I0G and that using a GRADE dictionary improve the "geometric quality", as assessed using Mogul compared to the original PDB model. As GRADE uses Mogul data to derive many of its restraints the average Z-score value is lower for GRADE.

From these results we can conclude:

Using QM with the RM1 method only marginally slows a BUSTER refinement.
That both refinements show good improvements over the PDB entry,
and that looking at either the partB or partC maps shows that quite a bit of density for ordered solvent has appeared in particular next to atom O20 of the ligand.
Although the Rfree is marginally higher in the QM run, the differences are small and likely to be due to slightly different X-ray weight adjustment. In terms of the protein structure the two refinements are in essence identical (for instance compared to the effect of using TLS).
So the exact restraint method used for the ligand does not effect the protein refinement - as really would be expected.
That both QM (RM1) and refinement with a GRADE dictionary produce good ligand geometry as assessed by Mogul.

(E) Assessing ligand strain energy and motions for the QM result

In part (B) we saw how it is possible to monitor the QM energy for the ligand during a refinement relative to its initial position. However, what we really want to assess is the strain energy of the ligand. To do this we want to find what the quantum energy of the ligand in its refined position is with respect to the minimum energy conformation of the ligand. Normally this can be found by doing a "geometry optimization" of the ligand in isolation from its position in the refined complex.

First extract the ligand coordinates for the A site. This can be done with an editor or using the egrep command. N.B.: BUSTER requires that all PDB files have a CRYST1 line:


   egrep "I0G A|CRYST1" partB/refine.pdb > partB.IOGA.pdb

Then use gelly_refine to geometry optimize the ligand:

gelly_refine -p partB.IOGA.pdb \
    -l I0G.grade_PDB_ligand.cif -qm I0G \
   -d partE -o partE_out.pdb \
   -max 1000 -glim 0.01 -jiggle_xyz 0.01 > partE.log &

The run uses a dictionary and -qm IOG
It also increases the default maximum number of steps and reduces the convergence limit: this is useful to ensure that a minimum is actually reached. The -jiggle_xyz 0.01 argument jiggles the input coordinates by a small amount: this is useful to ensure the optimizer does not get stuck in a saddle point.
The command graph_autobuster_QM can be used to find the QM energy of the ligand in the optimization.

graph_autobuster_QM -d partE
...
# final raw QM energies   =   -405.27 kJ/mol

So relaxation drops the QM energy to -405.27 kJ/mol. This energy can be used as a reference point to find the strain energies for the partB result:

graph_autobuster_QM -d partB -qmzero -405.27
...
# final raw QM energies   =   -399.62  -399.61 kJ/mol
# -qmzero specified so report QM energies relative to -405.27
#     big   small   iters      QMener_01      QMener_02
#       5      33     346           5.65           5.66

So the strain energy of bound ligand is 6 kJ/mol relative to the isolated geometry-optimized form.
The optimized coordinates for the ligand partE_out.pdb (purple) can be compared to the refined ligand (orange) using COOT:
It can be seen the relaxation is principally a rotation of the phenyl ring. The 2Fo-Fc density in the refined protein defines the orientation of the ring. Relaxation of the ring would cause it to have a short contact with Leu 339.
The torsion angle C2-C14-C16-C17 measures the rotation of the phenyl ring. This torsion angle changes from -148.82 degs. (partB) to -103.10 degs.

It is interesting to compare these results with the Mogul distribution of similar torsion angles from CSD small molecule structures. In the refined structure the ligand phenyl torsion angle would be regarded as "unusual" on the edge of the expected torsion angle from related CSD structures :

But relaxation with RM1 allows shifts the torsion angle into the range to be expected for small molecule structures.

The relaxation relieves short contacts between the phenyl ring and rest of the structure.

Back to Ligand QM/force field top index page

Page by Oliver Smart and Ian Tickle from 14 Feb 2011, modified November 2013. Any questions regarding our software or this wiki should be directed to buster-develop@globalphasing.com