Metal coordination in dock config file.

[electronicsbyjulie]

A metal ion may now be specified in the config file to be added to the PDB before docking. It will act like an additional ligand, but with its coordinating residues frozen in place around it. The best-binding algorithm has been updated to treat a metal and its coordination residues as a single monolithic alignment target with a property of metallicity that favors sulfur atoms in ligands.

Also included are: increased flexibility of the ligand; reduced availability of backbone atoms for hydrogen bonding; an improvement to the viewer to avoid hiding side chains behind backbone cartoons.

[electronicsbyjulie]

🟩🟩🟩

[primaryodors]

Couple things with this PR...

The "sulfur kludge" is kind of antithetical to the design of PrimaryDock. Code shouldn't be tied to any specific element, just like we try not to tie code to specific amino acids (glycine being the exception since no other glycine-like α-amino acid is possible). Instead, we tie code to properties, such as the histidine flip which should work for any imidazole-containing AA etc. No element, other than perhaps hydrogen, is unique enough in all its properties to warrant coding for directly like this. How would a selenium-containing ligand behave in the real world under the conditions where this sulfur kludge would apply? It would be far better to specify "chalcogen within a specific electronegativity range" or even just "chalcogen" with electronegativity-specific feature multiplication factors rather than tying the code to a single element.

Similarly, commenting out the Cu...O coordination strength from bindings.dat is a poor way to deal with oxygen atoms coordinating with metal ions. Why wouldn't a ligand's hydroxy group coordinate to a Cu+ or a Cu2+? Could be the simulation is correct. If the journal article ran their simulation and found it does not coordinate, then there must be a documented reason why that can be used to inform PrimaryDock just as it was used to inform their docker. In that case, finding the answer and applying it to the code is more proper than willy nilly removing an interatomic binding potential simply because it was inconvenient.

Whatever reason would cause a hydroxy group to coordinate to other hydroxyls instead of to a metal, should also be the reason for not best-binding it to a metal "glom", cf. the sulfur kludge above.

As for the hydrogen bonding kludge, it can stay for now given the programming complexity involved in making H-bonding competitive. Going to start a new issue soon for adding that capability.

[primaryodors]

This may be of help:

https://chem.libretexts.org/Courses/Saint_Marys_College_Notre_Dame_IN/CHEM_342%3A_Bio-inorganic_Chemistry/Readings/Metals_in_Biological_Systems_(Saint_Mary's_College)/Oxygen_transport_and_Storage/Blue_blood_(Invertebrates_molluscs_and_arthropods)/Haemocyanin

Hard/Soft Acid/Base Theory (HSAB)

Hard/Soft Acid/Base (HSAB) Theory allows for predictions of which acids and bases (metals and ligands) prefer to interact and which are most stable. In HSAB Theory, acids are electron acceptors and are referred to as acids/metal ions, and bases are electron donors and are referred to as ligands. These acids and bases can be categorized as hard or soft based on their charge density (charge to size ratio), polarizability (ability for an induced dipole to occur), and the covalent vs ionic nature of their interactions. Hard acids and bases have a higher charge density, are less polarizable, and/or have bonds that are more ionic; these tend to have a high charge and a small radius. On the other hand, soft acids and bases have a lower charge density, are more polarizable, and/or they have bonds that are more covalent; these tend to have a low charge and a large radius. HSAB can be used to predict which different metal-ligand pairs prefer to interact; hard acids prefer to coordinate with hard bases, soft acids prefer to coordinate with soft bases, and borderline acids prefer borderline bases. Essentially, the more “like” the characteristics of a metal and ligand, the stronger the preference of the ligand for that metal. This means that if there is a soft acid present, Au+ for example, in the active site of a molecule, it will prefer to bind to a soft base, such as S2− .10 In general, when you move down the periodic table, softness increases, and as you move right softness increases. Additionally, all of the bases are nonmetal electron donors (typically anions), and the acids are mostly metal ions.

As seen in figure 2, the active site of hemocyanin contains two copper ions, each bound to three histidine (His) residues. When the Cu(I) ion binds to oxygen it is oxidized to Cu(II). Cu(I) is a soft acid, Cu(II) is a borderline acid, and the imidazole ring of histidine is a borderline base. Atoms or molecules that are classified as having intermediate characteristics are called borderline, as they do not lie at either extreme of the hard/soft spectrum; they are harder than metals that are categorized as “soft” and softer than metals categorized as “hard”. Nitrogens are typically classified as hard bases, however if a nitrogen is contained within an aromatic ring, the entire ring is more polarizable. For this reason, imidazole nitrogen atoms are classified as borderline bases. Because the copper ion alternates between a soft acid (Cu(I)) and a borderline acid (Cu(II)), it makes sense that it is bound to a borderline base. Borderline histidine is a better match for borderline Cu(II) than soft Cu(I), however His is still a better match for Cu(I) than a set of ligands that is harder. It is important to recognize, however, that there are other factors that affect metal ion selectivity, such as ionic size and ligand field stabilization energy, that must also be considered.

[electronicsbyjulie]

https://doi.org/10.1002/anie.201007100

[electronicsbyjulie]

This article in favor of HSAB gives numbers that would seem to predict that Cu2+ would preferentially coordinate to O rather than S: https://doi.org/10.1021%2Fja00364a005

[primaryodors]

The question of HSAB vs. other methods warrants further study.

For this PR, it is sufficient to code the following: Using the electronegativities of the ligand atom (El) and of the metal (Em), compute 0.5 + 0.5 * cos((El + Em) / 2 - 2.25) and use that as the multiplier for the ligand-metal coordination strength.

[electronicsbyjulie]

The latest commit assumes that the metal ion is something with a higher electronegativity such as Cu. For olfaction, this might be a safe assumption, but for versatility, the best-binding algorithm should be rewritten to use potential binding energies between atoms instead of guessing with constant values for things like polarity, pi capability, charge, etc.

[primaryodors]

Okay to go ahead with the potential binding energies method in this PR if it isn't too big a programming task.

[electronicsbyjulie]

This is working currently:

PROT pdbs/OR2/OR2M3.upright.pdb
TEMPLATE pdbs/OR51/OR51E2.8f76.pdb
BRIDGE 6.57 45.51
LIG sdf/3-mercapto-2-methylpentan-1-ol.sdf

# https://doi.org/10.1007/s00018-019-03279-y
CEN RES 3.33 3.37 3.40 5.41 5.42 6.55
MCOORD Cu +1 5.42 5.43 5.46

SIZE 7.0 7.5 7.0
SOFT TMR2 TMR4 TMR5 TMR6 TMR7

EXCL 1 57       # Head, TMR1, and CYT1.
WET

POSE 10
ELIM 99
PROGRESS

FLEX 1
FLXR 3.33

ITERS 50

OUT output/%p_%l.active.dock
APPENDPROT

[electronicsbyjulie]

This is the nonsense that it's pulling for glom pairing currently (using a thiolated ligand and after removing Met206(5.46) from metal coordination):

Final pair assignments:
atom_glom[ C1 C2 C3 C4 C5 C6 H9 H10 H11 H12 H13 H14 H15 H16 H17 H18 H19 H20 ]-residue_glom[ Cys202 Cys203 ]
atom_glom[ O7 H21 ]-residue_glom[ Cys202 Cys203 ]
atom_glom[ S8 ]-residue_glom[ Arg261 ]

[electronicsbyjulie]

This PR is turning into another monster.

[electronicsbyjulie]

TODO:

Histidine flips and tyrosine flips (HH atom).

All the places where multimol_conform() is currently being called, should be replaced with calls to conform_molecules() instead. In cases where background molecules or molecules allowed to clash are specified, dedicated flags can be defined in MovabilityType.

Once all calls have been switched over, the monstrous multimol_conform() function and its myriad compiler constants can be removed.

All metal-binding code not pertaining to this PR should also be removed.

[electronicsbyjulie]

Electronegativity changes with ionization charge. There's a formula in this paper: https://doi.org/10.1021/jp062886k that can be used to calculate an ion's electronegativity on the fly, so long as we add another dat file for information about ions, ionization energies, and crystalline ionic radii.

[primaryodors]

Ion data would make a great topic for a future PR. Would like to get this PR wrapped up and ready to merge.

[electronicsbyjulie]

🟩🟩🟩

[electronicsbyjulie]

Let's also make sure the MCOORD command in Pepteditor works the same as the MCOORD parameter in PrimaryDock, i.e. they should both use the new metal coordination code.

Later, we'll be using a command line utility to identify and remove obsolete functions, and the above step will be helpful for that.

They're different formats with different capabilities. The amount of development involved would deserve its own PR. Recommend adding a deprecated notice to the Pepteditor MCOORD and ultimately replacing it with some kind of all-purpose atom/molecule coordination/prosthetic group feature, i.e. something that could add not only metal ions but also retinal, heme, etc.

[primaryodors]

Shouldn't deprecate a feature without having a replacement ready.

[primaryodors]

🟩🟩🟩

[electronicsbyjulie]

Drat. The side chain flexion stuff from #217 seems to have come undone with this PR. Even metal coordinating residues are flexing willy nilly.

[primaryodors]

This is the kind of thing that should be tested before marking the PR ready for review!