MDL aromaticity model doesn't recognize aromatic atoms in HIS and TRP

[j-wags]

This is a "bookmark" issue for something that we won't necessarily fix (because it isn't quite broken), but which can be used as an aggregator for likely future user reports.

Summary

The MDL aromaticity model doesn't recognize TRP and HIS sidechains as being aromatic. We are limited to using the MDL aromaticity model because it's the only one that both OpenEye and RDKit implement. This may be a problem because some of our small molecule and protein force fields depend on aromaticity information to assign physical parameters, and because in general users will be surprised by the assignment of aromaticity to these sidechain atoms.

Reproducing case

Histidine isn't aromatic under the MDL aromaticity model:

from rdkit import Chem
rdmol = Chem.MolFromSmiles('C1=C(NC=N1)C[C@@H](C(=O)O)N')
Chem.Kekulize(rdmol, clearAromaticFlags=True)
Chem.SetAromaticity(rdmol, Chem.AromaticityModel.AROMATICITY_MDL)
print(sum([1 for bond in rdmol.GetBonds() if bond.GetIsAromatic()]))

0

But if we switch the model to the RDKit default, histidine is aromatic

from rdkit import Chem
rdmol = Chem.MolFromSmiles('C1=C(NC=N1)C[C@@H](C(=O)O)N')
Chem.Kekulize(rdmol, clearAromaticFlags=True)
Chem.SetAromaticity(rdmol, Chem.AromaticityModel.AROMATICITY_DEFAULT)
print(sum([1 for bond in rdmol.GetBonds() if bond.GetIsAromatic()]))

5

For tryptophan, the MDL aromaticity model identifies the 6-membered ring as aromatic, but not the 5-membered ring.

Consequences

• SMIRNOFF protein force fields using the MDL model will need to either define HIS and TRP sidechain parameters without using aromaticity, or explicitly match to aliphatic interpretations of the sidechain. Due to issues with resonance, we already do the former for HIS, but we match to a kekulized representation of the five-membered ring of TRP [1].
• This may complicate handling of protein-specific torsions (and potentially post-translational modifications) in ways I can't really articulate well.
• Toolkit users may be surprised that their TRPs and HISes don't have their aromaticity labeled "correctly".

[1] An example TRP bond parameter from the ff14sb port version 0.0.2, where part of the sidechain is aromatic but part is aliphatic

<Bond smirks="[H][c:1]1c(c(c2c([c:2]1[H])C(=C(N2[H])[H])C([H])([H])[C@@]([H])([C]=O)[N][H])[H])[H]" length="1.4 * angstrom" k="938.0 * angstrom**-2 * mole**-1 * kilocalorie" id="A14SB-MainChain_TRP-CA_CA"></Bond>

Current plan

I'm going to plan on removing aromaticity/formal charge/bond order from TRP in our ff14sb port (and instead just to ~ bonds): openforcefield/amber-ff-porting#39

I'll update the OpenFF biopolymer tests to expect exactly the output from using the MDL aromaticity model.