Lit. Sem. 2010/6/30
Tatsuhiko Yoshino (M2)

QSAR in Catalysis ~ In Silico Catalyst Design ~Contents

1. Introduction : What is QSAR ?
2. Pioneer Work : Predicting Model for Tsuji-Trost Allylation
3. CoMFA : Analysis of Asymmetric Diels-Alder Reaction
4. QM-QSAR : Work of Prof. Mariza C. Kozlowski
5. Neural Network : Non-Linear Regression Methodology
6. Summary

1. Introduction : What is QSAR ?
1.1 Concept
QSAR = Quantitative Structure Activity Relationships
Method to predict "activity" of target compound "quantitatively" from caluculated parameters.
Major concept in drug design. (Activity = IC50 etc..)

QSAR = Quantitative Structure Asymmetry Relationships
QSSR = Quantitative Structure Selectivity Relationships
1.2 Descriptor
Descriptor : Variables (parameteres) which describe the feature of molecules.
Descriptors should be obtained by experiment or calculation.
Examples : Melting point, Log P, Dipole-moment, Bond-length, Dihedral-angle, number of functional group...
QSAR model correlates Descriptor and Activity in quantitative manner.

1.3 Training/Prediction
Fitting from Experimental Results (Training)
Compound 1
Compound 2
Compound 3D11005075D2102D30.150.60
1.00IC5030
30.5

Linear Correlation Modellog (1/IC50) = c0 + c1D1 + c2D2 + c3D3 + ...

Coefficients (cm) can be obtained by mathematical method.
(Multi-regression by least square, PLS analysis and so on)
Validation
Cross validation by LOO (leave one out) or LSO (leave several out) is usually employed for validation.
LOO (Leave One Out)
Correlation model is re-calculated with traning set where one compound is exclueded.
Then, activity of excluded compound is predicted using the new model and compared with real value.
LSO (Leave Sereral Out)
Almost same as LOO. Several compounds were excluded in this case.
Prediction

Activity of other compounds can be calculated using model equation.

1/15

1.4 What is Important ?
The most important point is to select the appropriate descriptors to describe molecular structre.
For drug discovery, CoMFA is one of the most general method.
Application of QSAR methodology to predict enantio-selectivity is main topic in this seminar.

1.5 Advantages and Disadvantages
Compared with Ab initio Transition State Calculation
Advantages : Rapid calculation for prediction.
Easy to understand what is important for selectivity.
Disadvantages : There is no theoretical, chemical guarantee. It is only statistical result.

2. Pioneer Work : Predicting Model for Tsuji-Trost Allylation
Norrby, P-O. et. al. "Steric Influences on the Selectivity in Palladium-Catalyzed Allylation"
Organometallics, 1997, 16, 3015.

2.1 Regio-Selectivity in Tsuji-Trost Allylation
EtOOC

COOEt
R1

R1

R2

R1

R2
orNu

R2
Nu

Pd

2.2 Reaction Conditions
Sodium diethyl malonate + E-allylic acetate in DMF

AllylsLigandsCatalytic and stoichiometric (using isolated -allyls)
readtions are performed.

E,Z-Isomerization
In catalytic reaction : fast (Boltzmann distribution by calculation)
stoichiometric reaction : slow

Products2/15

2.3 QSAR Model with Molecular Mechanics (MM2)
Final Selected Discriptors of Steric Feature of -allyl Pd"Selectivity" was conveted to Gibbs energy.
(Reaction Rate = K*exp(-Ea/RT) )Linear model(Structures of -allyl Pd were generated by MM2 method)

Discriptors
1) Breaking Pd-C bond length2) Dihedrals N1-N2-Pd-Ct
3) Dihedrals Du-Pd-Cc-Ct4) Steric Interaction with Nu

Steric Interaction with Nu : Set Ar probe atom to Nu position (figure) and calculate the increased energyRegression was performed with Levenberg-Marquardt algorithm (normal multi-regression)

2.4 Results and Discussion

Cross Validation Value : LOO Q2 = 0.86
LSO Q2 = 0.87

·The most important factor (descriptorw that coefficient has largest absolute value) is Pd-C bond.
·About result of cross validation of 1f which has unsymmetrical allyls and chiral ligand
1) Error occured because the fact that crossover between enantiometric path cannot take place.
2) Error occured because asymmetrical electronic effect was neglected.
Modification for these problem did not give better models.
Just a possible Error in MM2 system ?
3/15

3. CoMFA : Analysis of Asymmetric Diels-Alder Reaction

Lipkowitz, K. B. et. al.
"Computational Studies of Chiral Catalysts: A Comparative Molecular Field Analysis of an Asymmetric Diels-Alder
Reaction with Catalysts Containing Bisoxazoline or Phosphinooxazoline Ligands"
J. Org. Chem. 2003, 68, 4648.

3.1 Method
Comparative Molecular Field Analysis (= CoMFA) is now widely used method for drug design.
This method was first reported in 1988 in JACS.
Ref) Cramer, R. D. III et. al. "Comparative Molecular Field Analysis (CoMFA). 1. Effect of Shape on Binding ofSteroids to Carrier Proteins" J. Am. Chem. Soc. 1988, 110, 5959.

Key for CoMFA
1) All the analyzed compounds are set in grid space in appropriate manner.
Then interaction energy between "probe" and each coumpound at all grid points.
These energy values are used as descriptors.
2) Generated huge number of descriptors are analyzed by "PLS-regression" technique.
Grid-based Descriptor
Probe : sp3 C+ atom is often utilized. (Other probes can be also utilized.)Interaction energy is calculated as sum of van der Waals energy and Coulombic energy.
van der Waals Energy : Tripos Force Field ( 1/r6 1/r12)
Coulombic energy : 1/r, Atomic charge by *Gasteiger-Marsili method (from orbital electronegativity)
Other calculation methods can be used.

* Tetrahedron 1980, 36, 3219.

Grid Point·Enough range to cover all atoms of target compounds·Grid points with too high steric energy is cut-off.
·Grid points with too small standard deviation is eliminated.
·Grid space is usually 1.0-2.0 Å.
PLS-Regression
PLS = Partial Least Square is a regression method
suitable for models with...
·No clear reasonable relation of variables and property.
·Large number of variables (descriptors), often larger thannumber of samples.
Overfitting is a big problem in such a case, if normal multiregression is employed.

4/15

PLS-Regression : m variables (descriptors), n samples, 1 output models
yi : output
xij : variables j of sample itik : latent variables (LC of x)
w : weight coefficients for X to Tc : coefficients
Calculation of w1) For ti1, {w1j} is obtained to maximize
covariance (colinearlity) of {yi} and {ti1}.
2) For ti2, {w2j} is obtained to maximize
covariance (colinearlity) of {yi-c1ti1} and {ti2}.
3) For ti3, {w3j} is obtained to maximize
covariance (colinearlity) of {yi-c1ti1-c2ti2} and {ti3}.
Ref) http://cse.naro.affrc.go.jp/iwatah/index_j.html
PLS4) This procedure is repeated to reach r.

3.2 Application to Asymmetric Catalysis
Taget Reaction and Reaction ConditionsOOcatalyst
NO +ONOCatalysts

O

Experimental results were extrated from repotedpapers.
The most optimized reaction conditions were used.
Differeces in temperature, solvent etc were notconsidered.

Catalyst Structure
Initial Structures : CSD or Built in Spartan
Optimized by PM3tm method (semiempirical MO)
Alignment
Least-square fitting of oxazoline rings.
Validation
1) Internal Cross Validation by LOO.
In some cases, only internal validation is not enough.
Ref) Golbraikh, A; Tropsha, A. J. Mol. Graphics
Modell. 2002, 20, 269.
2) External Validation
At ramdom, 10,15,18,22 were excluded and used asexternal prediction set (LSO).
Best model
Field : both (Steric and Electronic)
Energy cut-off : S/E =30/20 kcal/mol
dielectric functions 1/r2
Probe : C+sp3
Latent variables : 6

Counter anions are not described,
but considered in structure optimizing calculation.

1) All catalysts model : r2cv = 0.833
2) LSO catalysts model : r2cv = 0.785, external r2 = 0.94
Golbraikh Tropsha Criteria : fullfilled
5/15

Visualized Results (STDEV*COEFF contour plot)
All Catalysts

Good Catalyst

BBAA

96% eeContribution of Each Factor to Selectivity
Steric : 60 70%
Steric factor is more important.
Electronic : 30 40%

Bad Catalyst

For Higher Selectivity

BA

A : Steric hindrance should be increased.
B : Steric hindrance should be decreased.

10% ee6/15

4. QM-QSAR : Works of Prof. Marisa C. Kozlowski
4.1 Method
Dixon, S. L.; Merz, K. M. Jr. et. al. "QMQSAR: Utilization of a Semiempirical Probe Potentialin a Field-Based QSAR Method" J. Comp. Chem. 2004, 26, 23.
Compared with CoMFA...
Similar Point : Energy values at grid points are calculated and used as descriptors.
Different Points : Energy is calculated as Probe Interaction Energies (PIE) by quantum mechanical method.
Regression is performed using n-variables regressions by simulated-annealing.
Probe Interaction Energies (PIE)
Probe : a positively charged carbon 2s electron
Potential from nucleus

Potential from electrons

M : Considered atoms: nucleus
si : wave functions of probe: wave functions of basis inc : coeffiecents of LCAO

All calculations are performed PM3 (semiempirical MO) method.

Regression
From several thousands of descriptors, n (2,3,4,5...) descriptors which give good f itting are selected.
How to select optimal descpriptors
Simulated Annealing Ref) Science 1983, 220, 671.
Review) Eur. J. Oper. Re. 1990, 46, 271.
Solutions for optimization problems in NP-hard class.
(impossible to solve in polynomial-time by deterministic algorithmm, if N NP)
Solution time is polinomial, but there is no guarantee to always give the right answer.
Mimic of annealing process.
1) T (Temperature) and initial state is set.
2) "State" is changed to another neighbour state "stochastically" as following.
If next state is better than now, state transition occured.
If next state is worse and temperature is enough high, state transition occured. (A)These operations are repeated.
3) Temperature is decreased.
4) Repeat 2)-3)
By process (A), probability of wrong answer
which is "local minimum" decreased. (Think of start from P.)
EP(A)local minimum(wrong)
Oright answer

S
7/15

4.2 QM-QSAR Approach for Predicting the Selectivity of Asymmetric Alkylation
Kozlowski, M. C. et. al. "Quantum Mechanical Models Correlating Structure with Selectivity: Predicting theEnantioselectivity of -Amino Alcohol Catalysts in Aldehydee Alkylation" J. Am. Chem. Soc. 2003, 125, 6614.
Target Reaction
catalystOPhH

Et2Zn

Et OH
PhH
99% eeTransition StateDAIB
Noyori, R. et. al. J. Am. Chem. Soc. 1986, 108, 6071.
Calculated TS (using Me2Zn) at the RHF/3-21G/Zn level.

Among possible 4 TS (syn/anti, R, S),anti S is the most favored.
Catalysts Set

Noyori, R. et. al. Organometallics 1999, 18, 128.
Catalysts Structure Optimization
Optimized Ground States : planar Zn
not suitableTransition States (PM3)
Aligned fitting Zn-O-Zn atoms

Method·PIEs (descriptors) are correlated with G values which can be converted to ee.
G = RT ln K, where K is enatio metric ratio.
·2 PIEs are selected to give the best fitting model by simulated annealing and normal least-square is used.
G = a + c1(PIE1) + c2(PIE2)
Best Individual Model
Or all acceptable PIEs pairs are weight averaged.
Averaged Model
Effect of Grid Space0.7 Å grid space show good convergence.

CC = Correlation Coefficients
Describe how well the prediction set selectivity
order is calculated.
8/15

Results
AAB
B

averaged model

best model

For good Selectivity.... A: more PIEs
B: less PIEs

PIEs : Electron rich area : decreased
Near Nucleus : inceased

Only minutes of computing gave good models!

4.3 Further Prediction for New Catalysts
M. C. Kozlowski et. al. "A Priori Theoretical Prediction of Selectivity in Asymmetric Catalysis: Design of Chiral
Catalysts by Using Quantum Molecular Interaction Field" Angew. Chem., Int. Ed. 2006, 45, 5502.

MethodImprovement : Ground states of dimeric catalysts gave a good model.

Monomer : Trigonal Zn ... Bad model
Dimer: Tetrahedral Zn ... Good model

Training SetPrediction Set

Trans amino alcohols with tetrasubstituted chiral
centers were included.

9/15

Resultss"Leave-two-out"
All 153 combinations of16 catalyst from T1-T18
gave 153 models.
SD = Standard deviation
CI = 95% confidence interval

Summary of Prediction
Leave-two-models gave better correlatioon.

y=x

4.4 G-QSAR Approach for Asymmetric Lithiation of N-Boc-pyrrolidine
Kozlowski, M. C et. al. "Is the A-Ring of Sparteine Essential for High Enantioselectivity in the Asymmetric LithiationSubstitution of N-Boc-pyrrolidine ?" J. Am. Chem. Soc. 2004, 126, 15473.
Target Reactions
1) chiral aminesBuLi;N
Boc 2) TMSCl+TMS
N
Boc
(S)

TMS
N
Boc
(R)

( )-sparteine : (S) >90% ee(+)-1: (R) >90% ee(+)-1 is easily synthesized andworks as (+)-sparteine surrogate.

chiral aminesABC

DCO'Brien, P. et. al. J. Am. Chem. Soc. 2002, 124, 11870.
BA

Mechanistic Investigation of ( )-Sparteine/sBuLi System
Wiberg. K. B.; Bailey. W. F. J. Am. Chem. Soc. 2001, 123, 8231.
10/15

Is A ring essential ?

What is important for selectivity ?
(S) 95%ee(S) 21% ee(R) 35% eeTraining Set

Results
No good correlation model was obtained by
QM-QSAR with PM3 calculation.
More precise calculation of PIEs was necessary.
G-QSAR : PIEs can be calculated using...
appropriate method (HF, MP2, B3LYP)
appropriate basis set (3-21G, 6-31G*, 6-31+G**...)
with Gaussian program.
Optimized Method
Structure Optimization : HF/3-21G* (ab initio MO)
PIE calculation : B3LYP/6-31G** (DFT)
2-variables model
LOO Cross Validation : r2cv = 0.68
Correlation Coefficient (CC) = 0.82
By optimized model, ( )-4 was predicted to give(R) product in 22-25% ee (exp. (R) 35% ee)

11/15

For better selectivity...
(A) More PIEs
(B) Less PIEs
(B)

These grid points are located above/below A ring.
1) Large group below A ring : good(A)2) Large alkyl group above A ring : Bad(B)3) Ph group above A ring : good(B)(A)

Structure aroud A ring seemed esseintial !

Other application example : Kozlowski, M. C. et. al. Org. Lett. 2006, 8, 1565.

5. Neural Network : Non-Linear Regression Methodology
5.1 Method : Artificial Neural Network Model
Serra, J. M. et. al. "Can artificial neural networks help the experimentation in catalysis?" Catalysis T oday 2003, 81, 393.

General Concept

Input

OutputArtificial Neural Network (ANN) works as "black box",
which gives "Output" from "Input"
even if the correlation is extremely complex and unknown.
"Black box" is programed to mimic a neural network(brain).

Artificial Neural Network
Artificial Neuron Modelinput x

output y

Linear model : y = x
These are called "activation function".
Non-linear model : y = 1/(1+e-x )
x -x
x -x
y = (e e )/(e e )Network Model : Multi-Layer Perceptrons
Weighted sums of previous layers' outputs are used as next inputs.
next input of neuron (i)xi = wi1y1prev + wi2y2prev + wi3y3prev + ... + wijyj prev
This calculations are performed for all next neourons.
input for next neurons : X = t {x1, x2, ... , xi}
previous out put : Yprev = t {y1prev, y2prev, ... , yjprev}

X = WYprev

weight matrix W = {wij}next output Y = g(X)(g : activation function)

Weight matric W should be optimized to give good correlation.
(N-1) Matrixes exist in N layers model.
12/15

Backpropagation - Conceptx1
Training signal
x2
When a teaching signal (X, Y)
is given, weight matrixes W are
modified to minimize errors.

Error
x3

y'1

y1x4

y'2

y2
Teaching signalx5x6
Advantage ·Applicable to many problems where theoretical anaysis or linear-regression is difficult.
Disadvantage ·It it impossible to obtain theoretical or qualitative information from the results.

5.2 QSAR Investigation of Heck Reactions
Farrusseing, D.; Rothenberg G. et. al. "Combinatorial Explosion in Homogeneous Catalysis: Screening 60,000 CrossCoupling Reactions" Adv. Synth. Catal. 2004, 346, 1844.
Target Reactions·412 Reactions were collected to analyze from
reported papers with various conditions.
·Activity : log(TON) and log(TOF)

Descriptors
Initial Set (76 descriptors)
Steric descriptors : MW, Surface, Volume, Tolman's cone angle, Solid angle and related parameters etc
Electronic descriptors : Hammett constant, HOMO, LUMO, GAP, Dipole moment, Chages on ligating atoms etc
Others : Pd loading, Pd precursor, reaction time, Temperature
Selected Descriptors Set (reduced by Relief Algorithm and Principal Component Analysis)
For TON (17 descriptors)
R1(Halide) : HOMO, LOMO, GAP, Socc
R2(Olefin) : LUMO, GAP, dipole, A
Ligand : q, HOMO, LUMO, GAP, Socc
Solvent : qOthers : Temp, Pd loading, Cat. precursor
q : chage on ligating atom

For TOF (20 descriptors)
R1(Halide) : Hammetp(+), Hammetp( ), V
R2(Olefin) : HOMO, LUMO, V, S(ethylene)/S, O, dipole
Ligand : q2, HOMO, LUMO, Socc, A, Rmax
Solvent : Omax
Others : Temp, Pd loading, Cat. precursor, Time
Socc : percentage of sphere occupation

13/15

Tolman's Cone Angle ( , T)

Solid angle(O,)

Tolman, C. A. Chem. Rev. 1977, 77, 313.

Solid angle is reflected by the shape of ligand.
Solid angles for radial profile (Omax, A, Rmax)
Omax

A

A = Area under the curve

White, D. et. al. J. Organomet. Chem. 1994, 478, 205.

Generated ANNs

For TON : 11 nodes and 3 nodes
For TOF : 15 nodes and 10 nodes in the 1st and 2nd hidden layers.

Reluts

A : TON by ANNs
B : TOF by ANNs
C : TON by linear regression model
Inset : residuals

ANNs > Linear regression model ?

14/15

Classification Problem

ModelsTree : Classification Tree Models
LDA : Linear-Discriminant-Analysis
Positive/Negative Threshold
Log(TON) = 2 (TON = 100)
Log(TOF) = 1 (TOF = 10)

Computational Screening of 60,000 Heck Reactions
61 PR3 type ligands

4 olefins

4 aryl-X

5 catalyst precursors

4 solvents

3 Pd loadings = 58,560 conditions

log(TON)PC1 is mainly correlated with Pd loadings andelectronic descriptors of R2.
PC2 is mainlys correlated with ligand electronic
descriptors.

6. Summary
QSAR Approach
Advantages·Short time calculation·Easy to extract what is important·Easy to search a good catalyst

Ab Initio Calculations of Transition States
Advantages·Based on reaction mechanism.
·Many samples are not needed.

Disadvantages·Not based on reaction mechanism·Only statistical estimation·Many samples are necessary for good model.

Disadvantages·Long time calculation·Difficult to predict what is important and good
catalyst without intuition.

15/15