Chemistry · Analytical and physical chemistry · Year 2

Computational chemistry and docking

How is a new drug “tried out” on a computer before it is synthesised in a flask?

The situation

At the Institute of Bioorganic Chemistry of the Academy of Sciences of Uzbekistan in Tashkent, substances isolated from plants and their derivatives are studied as drugs. Before synthesising a new compound, a chemist “builds” it on the computer: draws the molecule, obtains its 3D model, minimises its energy, calculates its vibrational spectrum and “inserts” it into the pocket of an enzyme. This work is done in programs such as Avogadro, HyperChem and Gaussian. How accurately does the computer predict experiment, which method suits which problem, and what in the results can be trusted?

Topic

Computer modelling in chemistry. Computer room of the Institute of Bioorganic Chemistry of the Academy of Sciences of Uzbekistan: before a new drug candidate reaches the laboratory it is “built” on a computer — a molecule is drawn, converted to 3D,…

Key ideas

Force field and optimisation
Quantum-chemical calculation
Docking and QSAR

Formulas

E = \sum k_{b} {(r - r_{0})}^{2} + \sum k_{θ} {(θ - θ_{0})}^{2} + \sum \frac{V_{n}}{2} (1 \pm \cos n φ) + \sum ε [{(\frac{r^{*}}{r})}^{12} - 2 {(\frac{r^{*}}{r})}^{6}] + \sum \frac{q_{i} q_{j}}{4 π ε_{0} ε_{r} r}

Force field (molecular mechanics): bond, angle, torsion, van der Waals and electrostatic terms.

𝐫_{k + 1} = 𝐫_{k} - s \nabla E

Geometry optimisation: steepest descent (in the conjugate gradient method the direction also “remembers” the previous step).

\hat{H} Ψ = E Ψ, 𝐅 𝐂 = 𝐒 𝐂 𝜺

The Schrödinger equation and the Hartree–Fock–Roothaan equations: orbitals as a sum of basis functions.

\tilde{ν} = \frac{1}{2 π c} \sqrt{\frac{k}{μ}}

Wavenumber of a harmonic vibration: force constant $k$ and reduced mass $μ$ (in a polyatomic molecule — the eigenvalues of the Hessian matrix).

Δ G_{bind} = R T \ln K_{d}

Docking: binding energy and dissociation constant ( $K_{d}$ small — the ligand binds strongly), $T$ = 298 K.

Theory

Computational chemistry: what is calculated and which programs are used

Computational chemistry calculates from theory the quantities that are measured in experiment: bond lengths, bond and torsion angles, the energies of conformers, bond energy and bond order, heat of formation, ionisation potential (Koopmans' theorem: $I \approx - ε_{HOMO}$ ), electron affinity, the energy gap (HOMO–LUMO), dipole moment and spectra. A calculated result is always compared with experiment — these are the “calculated” and “experiment” columns in the tables on the stage. There are two kinds of program: graphical (ISIS Draw, ChemWin, ChemDraw) — drawing formulae, naming, 2D pictures; and computational (Avogadro, HyperChem, Gaussian, ChemOffice/Chem3D) — 3D geometry, energy, spectra. A calculation starts from an initial geometry: the 2D drawing is “lifted” into 3D (VSEPR templates, a z-matrix — bond length, angle and torsion) and then optimised.

Molecular mechanics (empirical method) and optimisation

A molecule is treated as balls joined by springs: the force field is the sum of bond-stretching, angle-bending, torsion, van der Waals and electrostatic terms. Their parameters ( $k$ , $r_{0}$ , $θ_{0}$ , $V_{n}$ ) are fitted to experiment (IR force constants, electron-diffraction geometries, conformer energies) — this is parametrisation. $r_{0}$ is the ideal (strain-free) length; the real length in an optimised molecule differs from it slightly because of the influence of neighbouring atoms. Geometry optimisation is the search for the point where the energy gradient is zero: steepest descent is simple but slow (“zigzag”), while conjugate gradients and quasi-Newton (BFGS) are faster. Optimisation finds only a local minimum; for the global minimum a conformer search is needed (a systematic scan or many random starts). The force field in this experiment has been fitted so that it gives an ethane barrier of 12.0 kJ/mol and, for butane, gauche − anti = 3.8 kJ/mol. In Avogadro it is the same: “Auto Optimize” (MMFF94, UFF force fields).

Quantum-chemical methods: ab initio, basis sets, choosing a method

Ab initio (“from the beginning”) methods have no experimental parameters: in the Hartree–Fock method each electron moves in the average field of the others, and the orbitals are sought as a sum of basis functions (Gaussian functions $e^{- α r^{2}}$ ). As the basis grows, the energy decreases from above, by the variational principle, and approaches a limit (on the stage: the H atom, n = 1…8). But for helium even the Hartree–Fock limit lies 1.14 eV above the exact energy — this is the electron correlation energy; it is taken into account by MP2, CCSD(T) and density functional methods (DFT, for example B3LYP). Rules for choosing a basis: minimal (STO-3G) — fast, low quality; 6-31G* — with polarising d-functions, for many organic problems; diffuse (+) functions for anions and hydrogen bonds; cc-pVTZ for accurate energies. The cost of a calculation grows quickly: HF ~ $N^{4}$ , MP2 ~ $N^{5}$ , CCSD(T) ~ $N^{7}$ ( $N$ is the number of basis functions). Choosing a method: thousands of atoms — molecular mechanics; hundreds — semi-empirical or DFT; accurate thermochemistry — CCSD(T). Gaussian is the most widely used software package for ab initio and DFT calculations (geometry, frequencies, NMR, UV, thermochemistry).

Semi-empirical methods, charges and the UV spectrum

Semi-empirical methods simplify quantum mechanics and replace some of the integrals with parameters taken from experiment: the simplest is the Hückel π-electron method (parameters $α$ and $β$ ), then extended Hückel, CNDO, MNDO, AM1, PM3, PM7 (HyperChem, MOPAC). They are fast and give energetic quantities (heat of formation, ionisation potential), the charge distribution, the dipole, and a graphical display of the orbitals. The partial charges on atoms depend on the method: Gasteiger (electronegativity equalisation), Mulliken, ESP fitting; for water the Gasteiger charges give a dipole of 1.2 D, while liquid-water force fields (TIP3P, SPC/E) give 2.3–2.4 D (experiment, gas — 1.85 D): if you change the force field, the charges and results change. UV spectrum: in the Hückel method the HOMO→LUMO transition energy is $Δ E = (x_{HOMO} - x_{LUMO}) | β |$ ; when $| β |$ is fitted to butadiene (217 nm), the wavelength is overestimated for long polyenes — the method ignores bond alternation (more accurate: ZINDO, TD-DFT). Hückel π-orbitals and aromaticity are covered in the “The molecule: orbitals and symmetry” experiment.

Calculating and interpreting spectra

IR spectrum: the matrix of second derivatives of the energy (the Hessian matrix) is divided by the masses and diagonalised — the eigenvalues give the frequencies of the normal modes and the eigenvectors give the motion of the atoms; if the dipole changes during the vibration it is active in the IR (intensity $\propto {(\partial μ / \partial Q)}^{2}$ ), and if not, it appears only in the combination-scattering (Raman) spectrum. The harmonic approximation overestimates frequencies, so calculated frequencies are scaled (≈ 0.9 for HF, ≈ 0.96 for B3LYP) or the force field is fitted directly to experimental frequencies (this is how it was done for water in this experiment). UV — electronic transitions (TD-DFT, ZINDO); NMR — chemical shifts (the GIAO method; ChemOffice estimates ¹H and ¹³C shifts using increment — additivity — tables); mass spectrum — fragmentation pathways. A calculated spectrum helps in interpreting an experimental one: which line belongs to which vibration (in detail — in the “Spectra: determining structure” experiment).

Biological activity: QSAR, PASS, docking; limits of the model

Structure–activity relationships (SAR/QSAR): activity is predicted by regression from molecular descriptors (lgP, molecular mass, H-bond donors/acceptors, electronic effects) (the Hansch equation); for oral drugs, Lipinski's “rule of five”: M ≤ 500, lgP ≤ 5, donors ≤ 5, acceptors ≤ 10. PASS (V. Poroikov and colleagues, Moscow) gives, by a Bayesian method, the probabilities of hundreds of kinds of activity from the structural fragments of a molecule: if Pa (active) > Pi (inactive), the activity is possible, and Pa > 0.7 means high. Molecular docking is placing a ligand in the pocket of an enzyme: a search (Monte Carlo, genetic algorithm) plus a scoring function (steric fit, H-bonds, desolvation, entropy); programs such as AutoDock Vina. Virtual screening sorts thousands of compounds this way and chooses the best for synthesis. Limits of the model: the force field in this experiment is simplified for teaching and fitted to a few molecules; the ab initio calculation is only for atoms with s-type Gaussian functions; the enzyme pocket is not a real protein but a schematic model; the docking energies and $K_{d}$ are model values, not experimental ones. A calculated result is a hypothesis, and experiment confirms it.

Did you know?

In 2013 the Nobel Prize in Chemistry was awarded to M. Karplus, M. Levitt and A. Warshel “for the development of multiscale models of complex chemical systems” — in which molecular mechanics and quantum chemistry are combined.
The 1998 Nobel Prize in Chemistry was shared between Walter Kohn (density functional theory) and John Pople (computational methods of quantum chemistry; he was one of the founders of the Gaussian program).

Inside the experiment

The experiment starts from a real-life situation: you predict the result, change the controls, run the experiment and draw a conclusion. Inside are missions, theory, a quiz and a guided tour in which the camera moves and explains the key points.

More experiments in this section

Spectra: determining structure — Year 2
Alloys: phase diagrams — Year 2
Chemical kinetics and Arrhenius — Year 2
Catalysis and enzymes — Year 2
Colloids and adsorption — Year 2