1D pharmacophores in P450 models

During the development of SMARTCyp 2 we realized something profound, site-of-metabolism in CYP2D6 and CYP2C9 can be modelled by reactivity and a simple 1-dimensional pharmacophore. Simply by measuring the distance from the atom of interest to the pharmacophore we get a good variable to use for correcting the energy with.

So how come this simple correction works? It is rather well known that cytochromes P450s are quite flexible. A homology model of CYP2D6 has been shown to be more reliable for site-of-metabolism prediction through docking than the crystal structure, indicating a flexible site. The CYP2C9 structure has also been shown to be able to take on multiple conformations. Hence, a structure based model will require many different conformations to give accurate description of the binding, whereas a non-structure based model (i.e. SMARTCyp) simply implicitly includes this flexibility by assuming that there are very few restraints on the rotation of the substrate. SMARTCyp only uses the relative distance (no. of bonds) to the center of the molecule as a small penalty to the reactivity to prioritize sites that are closer to the end of the molecule than the center.


For the implementation of 1D pharmacophores we devised SMARTS strings that match the pharmacophores of CYP2D6 and CYP2C9 that matches the pharmacophores of the isoforms (strings below from SMARTCyp 2.2).

CYP2D6: positively charged amines
$([N][CX3](=[N])[N]) //guanidine like fragment
$([N^3X3H0]([#6^3])([#6^3])[#6^3]),$([N^3X3H1]([#6^3])[#6^3]),$([N^3X3H2][#6^3]) // primary, secondary, tertiary amines bound to only carbon and hydrogen atoms, not next to sp2 carbon
                                

CYP2C9: carboxylic acids and their bioisosteres
$([O]=[C^2][OH1]) // carboxylic acid oxygen
$([O]=[C^2][C^2]=[C^2][OH1]),$([O]=[C^2][c][c][OH1]) // vinylogous carboxylic acids (e.g. ascorbic acid)
$([n]1:[n]:[n]:[n]:[c]1) // tetrazole 1
$([n]1:[n]:[n]:[c]:[n]1) // tetrazole 2
$([O]=[C^2][N][OH1]) // hydroxamic acid
$([O]=[C^2]([N])[N]) // urea
$([O]=[S][OH1]) // sulfinic and sulfonic acids
$([O]=[PD4][OH1]) // phosphate esters and phosphoric acids
$([O]=[S](=[O])(c)[C][C]=[O]),$([O]=[C][C][S](=[O])(=[O])[c]) // sulfones next to phenyls with carbonyl two bonds away
$([O]=[S](=[O])[NH1][C]=[O]),$([O]=[C][NH1][S](=[O])=[O]) // sulfones bound to nitrogen with carbonyl next to it
$([O]=[C^2][NH1][O]),$([O]=[C^2][NH1][C]#[N]) // peptide with oxygen or cyano group next to nitrogen
$([OH1][c]1[n][o,s][c,n][c]1),$([OH1][n]1[n][c,n][c][c]1),$([OH1][n]1[c][n][c][c]1) // alcohol on aromatic five membered ring
$([O]=[C]1[N][C](=O)[O,S][C,N]1) // carbonyl oxygen on almost conjugated five membered ring
$([O]=[C]1[NH1,O][N]=[N,C][N]1) // carbonyl oxygen on fully conjugated five membered ring
$([nD2]1[nD2][c]([S]=[O])[nD2][c]1),$([nD2]1[c]([S]=[O])[nD2][c][nD2]1),$([nD2]1[c]([S]=[O])[nD2][nD2][c]1) // nitrogens in histidine-like 5-ring with sulfoxide/sulfone next to it
$([O]=[SX4](=[O])[NX3]) // sulfonamides

References:
2D6 model
Patrik Rydberg and Lars Olsen "Ligand-Based Site of Metabolism Prediction for Cytochrome P450 2D6" ACS Med. Chem. Lett., 2012, 3, 69-73
2C9 model
Patrik Rydberg and Lars Olsen "Predicting Drug Metabolism by Cytochrome P450 2C9: Comparison with the 2D6 and 3A4 Isoforms" ChemMedChem, 2012, 7, 1202-1209

SMARTCyp with ChemDoodle web components

Recently we released version 2.01 of SMARTCyp, which is the first version to use ChemDoodle web components.

We have done this for two major reasons, the first one is speed. Creating HTML output by writing text files with javascript is much faster than generating three png files for each molecule. While the old version output is still available in SMARTCyp, the default HTML output is HTML5 using ChemDoodle web components. Not only is it faster, it looks nicer too.

The second reason is platform compatibility, the new HTML output can also respond to input from touch-screens. A touch input on the structures (when using your mobile phone or your tablet device) will display atom numbers in the same way a mouseover effect does it on a standard computer. To get this working we simply made duplicate functions in the javascript. There are mousover and mouseout functions for displaying and not displaying atom numbers. The same functions also exist as touchstart and touchend functions. If you want the details, just check out the resulting HTML you get when you run a molecule on our web server.

When we did this we also changed the sketcher on our web server start page. It used to be a JchemPaint sketcher applet. Since this is java it simply does not work on all mobile devices (and Ipads). Now we instead use the ChemDoodle sketcher, which loads much faster, and is compatible for all mobile devices as well.

Chemdoodle web components tricks #3: Using mouseover functions to do multiple repaints on a canvas

Another step on the way to integrating ChemDoodle web components into SMARTCyp.

Here I show how you can use mouseover, and mouseout functions to redraw the canvas over and over again, and infinite number of times. The trick is to minimize the mouseover and mouseout functions, and use the drawChildExtras function to do all the drawings. Also, the first four lines and the last line in the drawChildExtras function are required to make this work (thanks to Kevin Theisen for helping me with those).

Here is the canvas, and below is the annotated code.



And here is the code:

 var tutorial3_testmol = new ChemDoodle.ViewerCanvas('tutorial3_testmol', 300, 300);
 tutorial3_testmol.specs.atoms_useJMOLColors = true;
 var caffeineMolFile = 'Molecule Name\n  CHEMDOOD08070920033D 0   0.00000     0.00000     0\n[Insert Comment Here]\n 14 15  0  0  0  0  0  0  0  0  1 V2000\n   -0.3318    2.0000    0.0000   O 0  0  0  1  0  0  0  0  0  0  0  0\n   -0.3318    1.0000    0.0000   C 0  0  0  1  0  0  0  0  0  0  0  0\n   -1.1980    0.5000    0.0000   N 0  0  0  1  0  0  0  0  0  0  0  0\n    0.5342    0.5000    0.0000   C 0  0  0  1  0  0  0  0  0  0  0  0\n   -1.1980   -0.5000    0.0000   C 0  0  0  1  0  0  0  0  0  0  0  0\n   -2.0640    1.0000    0.0000   C 0  0  0  4  0  0  0  0  0  0  0  0\n    1.4804    0.8047    0.0000   N 0  0  0  1  0  0  0  0  0  0  0  0\n    0.5342   -0.5000    0.0000   C 0  0  0  1  0  0  0  0  0  0  0  0\n   -2.0640   -1.0000    0.0000   O 0  0  0  1  0  0  0  0  0  0  0  0\n   -0.3318   -1.0000    0.0000   N 0  0  0  1  0  0  0  0  0  0  0  0\n    2.0640   -0.0000    0.0000   C 0  0  0  2  0  0  0  0  0  0  0  0\n    1.7910    1.7553    0.0000   C 0  0  0  4  0  0  0  0  0  0  0  0\n    1.4804   -0.8047    0.0000   N 0  0  0  1  0  0  0  0  0  0  0  0\n   -0.3318   -2.0000    0.0000   C 0  0  0  4  0  0  0  0  0  0  0  0\n  1  2  2  0  0  0  0\n  3  2  1  0  0  0  0\n  4  2  1  0  0  0  0\n  3  5  1  0  0  0  0\n  3  6  1  0  0  0  0\n  7  4  1  0  0  0  0\n  4  8  2  0  0  0  0\n  9  5  2  0  0  0  0\n 10  5  1  0  0  0  0\n 10  8  1  0  0  0  0\n  7 11  1  0  0  0  0\n  7 12  1  0  0  0  0\n 13  8  1  0  0  0  0\n 13 11  2  0  0  0  0\n 10 14  1  0  0  0  0\nM  END\n> \n07-08-2009\n';
 var caffeine = ChemDoodle.readMOL(caffeineMolFile);
 // get the dimension of the molecule
 var size = caffeine.getDimension();
 // find the scale by taking the minimum of the canvas/size ratios
 var scale = Math.min(tutorial3_testmol.width/size.x, tutorial3_testmol.height/size.y);
 // load the molecule first (this function automatically sets scale, so we need to change specs after)
 tutorial3_testmol.loadMolecule(caffeine);
 // change the specs.scale value to the scale calculated, shrinking it slightly so that text is not cut off
 tutorial3_testmol.specs.scale = scale*.8;
 //we are not showing atom numbers right away
 var showatomnumbers = false;
 
 tutorial3_testmol.mouseover = function(){
  //set atom numbers to be displayed
  showatomnumbers = true;
  this.repaint();
 }
 
 tutorial3_testmol.mouseout = function(){
  //set atom numbers to be hidden
  showatomnumbers = false;
  this.repaint();
 } 
 
 tutorial3_testmol.drawChildExtras = function(ctx){
  //four lines to make sure atom numbering, coordinates and scaling work correctly
  ctx.save();
  ctx.translate(this.width/2, this.height/2);
  ctx.rotate(this.specs.rotateAngle);
  ctx.scale(this.specs.scale, this.specs.scale);
  ctx.translate(-this.width/2, -this.height/2);
  //set the font size relative to the scaling of the molecule
  ctx.font = 'bold ' + 2.5*scale + 'px sans-serif';
  //center the atom numbers over the x,y coordinates of the atoms
  ctx.textAlign = 'center';
  ctx.textBaseline = 'middle';
  //iterate through all atoms
  for (var i = 0, ii=caffeine.atoms.length; i<ii; i++) {
   var atom = caffeine.atoms[i];
   if(showatomnumbers){
      //draw a white circle behind the atom number
      ctx.fillStyle = "white";
      ctx.beginPath();
      ctx.arc(atom.x, atom.y, scale*1.8, 0, Math.PI*2, true);
      ctx.fill();
      //draw the atom number in the color
      ctx.fillStyle = ChemDoodle.ELEMENT[atom.label].jmolColor;
      //the "+scale*0.7" is to adjust the text on the y-center of the atoms.
      ctx.fillText(i + 1, atom.x, atom.y);
   }
  } 
  //restore the ctx settings so that application next time will start from scratch
  ctx.restore();
 }
 tutorial3_testmol.repaint();

Why SMARTCyp works - challenges in reactivity models applied to P450s

During the last couple of years I've been developing the SMARTCyp methodology (www.farma.ku.dk/smartcyp) for prediction of site-of-metabolism in drug metabolism mediated by the cytochromes P450 enzyme family. Here I will describe why the SMARTCyp approach works, and why many other models applied to the same problem often fails when they try to include reactivity.

1. SMARTCyp is simple
   The application of the reactivity model is through simple fragment matching towards a library of SMARTS strings. Each SMARTS string describes a number of reference calculation and assigns the average energy of these reference calculations to any atom that matches.
   The advantage of this is that we get rid of the problem of generating many different 3D structures to make sure the conformation is correct, and we get rid of the problem of possible hydrogen bond donors/acceptors messing up reactivity calculations by causing strange structures (this can happen when the oxygen radical attacks a site close to a donor/acceptor).
   
2. SMARTCyp reactivities are very accurate
   For the so far 200+ fragments we have computed the reactivities are as good as can possibly be done today (DFT with B3LYP, large basis sets and transition states verified by frequency calculations). As long as the fragment of interest in a drug molecule is relatively similar, the predicted reactivity will be quite good (>1 kcal/mol error is unusual).
   This is because our calculations use a full heme model, and compute the transition state for the actual reaction mechanism, no simplified model or other assumptions are made. Other reactivity models use intermediates and also do not use a full heme model. Such models are relatively accurate for hydrogen abstraction reactions (that is hydroxylation of aliphatic carbon atoms or dealkylation reactions), but for any other type of reactions the models are usually pretty bad.
   
3. P450 enzymes are flexible
   While it has been shown that CYP3A4 is highly flexible and reactivity is the only major determinant of site-of-metabolism in this enzyme, it's my belief that most other CYPs also are more flexible than most people believe. The crystal structures of several of the human CYPs do not have active sites large enough to explain all the substrates known to be metabolized by them. Hence, the approximation that reactivity is very important is very good for most CYPs (however, it's not sufficient for CYP2D6 and CYP2C9). Unpublished data show that SMARTCyp can be applied to all isoforms except 2D6 and 2C9 with high accuracy. The two exceptions stem from the fact that these are the only two isoforms with charged amino acids in the active site, leading to large contributions from binding. But even these binding contributions can be described by a simple 2D model.
   We recently showed that even for CYP2D6 we can build a good model based on SMARTCyp reactivities using only three descriptors (implemented in SMARTCyp 2.0). While some may claim such a model is too simple, the simplicity is actually an advantage. Simple models are usually robust, the dependence on data sets can be smaller if you have fewer descriptors that you're trying to fit. And hence we have less noise in the data. When compared to ensemble docking supplemented with reactivities from SMARTCyp 1.5, the new SMARTCyp 2D6 model is as accurate,but with a fraction of the computational cost, and much much easier to apply for a medicinal chemist.