For this week’s journal club I presented this paper by Jian et al from the Zhang Lab. The paper tackles the problem of refining protein structure models using molecular dynamics (MD).

The most successful protocols for building protein structure models have been template-based methods. These involve selecting whole or parts of known protein structures and assembling them to form an initial model of the target sequence of interest. This initial model can then be altered or refined to (hopefully) achieve higher accuracy. One method to make these adjustments is to use molecular dynamics simulations to sample different conformations of the structure. A “refined” model can then be taken as a low-energy state that the simulation converges to. However, whilst physics-based potentials are effective for certain aspects of refinement (e.g. relieving clashes between side chain atoms), the task of actually improving overall model quality has so far proved to be too ambitious. 

The Method

In this paper entitled “Atomic-Level Protein Structure Refinement Using Fragment-Guided Molecular Dynamics Conformation Sampling,” Jian et al demonstrate that current MD refinement methods have little guarantee of making any improvement to the accuracy of the model. They therefore introduce a technique of supplementing physics-based potentials with knowledge about fragments of structures that are similar to the protein of interest.

The method works by using the initial model to search for similar structures within the PDB.  These are found using two regimes. The first is to search for global templates by assessing the TMscore of structures to the whole initial model. The second is to search for fragments of structures by dividing the initial model into continuous 3 secondary structure elements. From these sets of templates and the initial model, the authors can generate a bespoke potential for the model based on the distances between Cα atoms. By doing this, additional information about the likely global topology of the protein can be incorporated into a molecular dynamics simulation. The authors claim that this enables the MD energy landscape is therefore reshaped from being “golf-course-like” being “funnel-like”.  Essentially, the MD simulations are guided to sample conformations which are likely (as informed by the fragments) to be close to the target protein structure. 

A schematic of the FG-MD refinement procedure

 Does it work?

As a full solution to the problem of protein structure model refinement, the results are far from convincing. Quality measures show improvement in only the second or third decimal place from the initial model to the refined model. Also, as might be expected, the degree to which the model quality is improved is dependent on the accuracy of the initial of the model.

However, what is important about this paper is that, although small, the improvements made do exist in a systematic fashion. Previously, attempts to refine a model using MD not only failed to improve its accuracy but would be likely to reduce its quality. Fragment-guided MD (FG-MD) and the explicit inclusion of a hydrogen bonding potential, is not only able to improve the conformations of side chains but also improve (or at least not destroy) the global backbone topology of a model.

Dependence of the energy of a model on its TM-Score to the native structure. In black is the energy as measure using the AMBER99 energy function. In grey is the corresponding funnel-like shape of the FG-MD energy function.

This paper therefore lays the groundwork for the development of further refinement methods that incorporate the knowledge from protein structure fragments with atomically detailed energy functions. Given that the success of the method is related to the accuracy of the initial model, there may be scope for developing similar techniques to refine models of specific proteins where modelling quality is already good. e.g. antibodies.

Last wednesday I talked about this paper which is about antimicrobial peptides. First I should say even though many of the authors are French this particular group was unknown to me and it was only the subject tackled which caught my eye – no patriotism, as it should become obvious in the following paragraphs. More precisely, being a bit of an outlier in OPIG (as I spend all my time in the biochemistry building doing Molecular Dynamics (MD) simulations) I was curious to get feedback from OPIG members on the sequence processing method developed in this article.

So what are antimicrobial peptides (AMP) and what do we want to do with them?

These are peptides able to kill bacteria in a rather unspecific manner: instead of targeting a particular membrane receptor AMP appear to lyse bacterial membrane via poration mechanisms which are not totally understood but relies on electrostatic interactions and lead to cell leakage and death. This mode of action means it is difficult for bacteria to develop resistance (they can’t change their entire membrane) and AMP are extensively studied as promising antibiotics. However, the therapeutic potential of AMP is severely limited by the fact that they are also often toxic to eukaryotic cells – this is the downside of their rather unspecific membrane lytic activity.

If AMP are to become useful antibiotic it is therefore important to maximise their activity against bacterial cells while limiting it against eukaryotic cells. This is precisely the problem tackled by this paper.

The ratio between the desired and not desired activity of AMP can be measured by the therapeutic index (TI) which is defined in this paper as:

TI = HC50 / MIC

where HC50 is the AMP concentration needed to kill 50% of red blood cells and MIC is the minimal inhibitory concentration against bacteria. The higher the TI the better the AMP for therapeutic use and the authors presents a method to identify mutations (single or double) in order to increase the TI of a given sequence.

What did the authors do?

Their method relies on a new way to describe AMP properties, which the authors called the sequence moment and which was built to reflect the asymmetry along the peptide chain regarding the property of interest. The rationale underlying this is that there is a functional asymmetry in the peptide chain, with the N terminal region being “more important for activity and selectivity of AMPs” – this claim is based on a reference from 2000 in which it is shown that in the peptides studied mutations in the N-terminus region are more likely than that in the C-terminus region to result in a decrease in AMP activity.

The derivation of the sequence moment is interesting and illustrated in the figure below. The AMP sequence is projected on a 90 degrees arc, with the N terminus on the y axis. The metric of interest is then plotted for each residue as a vector (tiny arrows in the figure) with the same origin as the arc and the length and direction of which depends on the metric value and the orientation. The mean value of the metric for the AMP sequence is then obtained by summing the positional vectors together (big arrows).

With this method the authors are able to get a mean value for the metric of interest which also contains information on the distribution along the AMP sequence of the metric of interest.

Interestingly, the descriptor used by the authors is actually the angle difference between the arrows obtained when calculating the sequence moment using two different hydrophobicity scales (Janin’s and Guy’s) – see red and blue arrows in the figure above.

The authors claim that the cosine of this angle, termed descriptor D, correlates well with the TI of AMP according to the following linear relation:

TI = 50.1 – 44.8*D

Long story short the authors then describe how they implemented a software, called Mutator and available here, which takes a sequence as input and based on a set of 26 best AMP (defined as having a TI > 20 and less than 70% sequence homology between them) suggest single or double mutations to improve its TI based on the method and relationship above. They then test their predictions experimentally for 2 peptides (Ascaphin-8 and XT-7) and it appears that the analogues suggested by the Mutator do have a much higher TI than their parent sequence.

What do we think about it?

Although this paper presents a method which seemingly gives impressive results I have 2 major problems with it.

The first one would be that the D descriptor has no rationale nor clear signification which makes interpreting the results difficult. Which sequence property does this descriptor capture? At best it can be said that when D is small the two hydrophobicity scales differ widely for the sequence studied (arrows at a 90 degrees angle), whereas they agree when D is close to 1. It would then be necessary to go back to how the scales were derived to understand what is being picked-up here, if anything.

Second, there is no proper statistical analysis of the significance of the results obtained. As noted by OPIG members the peptides studied are all fairly similar,and the less than 70% pairwise identity rule used for the training does not guarantee much diversity. Essentially the algorithm is thus trying to make sequences which are not very different from the training set into even more so similar sequences and the type of residues mutated to achieve is limited due to the fact AMPs are rather short sequences and enriched in particular residues. Therefore one might be able to argue that any such mutation is rather likely to lead to an improvement of the TI – especially if the input sequence is chosen specifically for not being optimised for the desired activity (broad spectrum AMP). It would be important to design a proper null hypothesis control and measure experimentally whether the TI index of the analogues obtained are statistically significantly lower than those obtained with the Mutator software.

In summary it might sound a bit harsh but my personal opinion is that this paper is the kind of paper people who run Virtual Screening like a black box (see JP’s excellent previous post) will enjoy. Copy and paste a sequence, hit run and it seemingly gives you result. If you look hard enough for “a” descriptor that “correlate” with your data you will find one – especially if A) you don’t define what “to correlate” means, and B) your descriptor doesn’t have to mean anything. So the trouble is no one has the beginning of a clue about what is going on and , what is worse, before being able to think about it it would be first necessary to run proper controls in order to assess whether anything is actually going on.

So, my turn to present at Group Meeting Wednesdays – and I decided to go for the work by Scior et al. about pitfalls in virtual screening.  As a general comment, I think this paper is well written and tackles most of the practical problems when running a virtual screening (VS) exercise.  Anyone, who intends to either develop a method in this field or else is planning to run a virtual screening  exercise should read it.  I’ve often heard the phrase “virtual screening doesn’t work”, and that comes almost exclusively from people who run computational experiments as a black box, without understanding what is going on and by accepting all defaults for a specific protocol.  This paper highlights what to watch out for.  Of the author list, I’ve only met with Andreas Bender once at an MGMS meeting a few years back – his main PhD work was on molecular similarity.

The article describes pitfalls associated with four areas; expectations and assumptions; data design and content; choice of software; and conformational sampling as well as ligand and target flexibility.  The authors start off by arguing that the expectations are too high; people just run a VS experiment and expect to find a potent drug.  But this a rare occurrence indeed.  Below is a set of notes for their main points.

Erroneous assumptions and expectations

1. High expectations: main goal is to identify novel bioactive chemical matter for the particular target of interest.  Highly potent compounds desirable but not required.  Expectations too high.  Lead; single digit µM and Hit; < 25µM
   
2. Stringency of queries: strict vs loose search criteria.  Strict; no diversity, few good results returned.  Loose; returns many false positives.   Removal of one feature at a time – a kind of pharmacophoric feature bootstrapping which highlights which features are important.
3. Difficulty in binding pose prediction.  Taken from their reference [46], “For big (> 350 Daltons) ligands, however, there is currently very little evidence. We hope investigators will come forward with crystallographic confirmation of docking predictions for higher molecular weight compounds to shed more light on this important problem.”  This point is really interesting and tools, such as Gold, even have an FAQ entry which addresses this question.
4. Water: hydrogen bonds are mediated by water which are often visible in the crystal structure.  Hard to predict exact number, position and orientation.  Realism of model at a cost of computational resources.
5. Single vs. multiple/allosteric binding pockets: Sometimes binding site is not known yet we always assume that ligand binds to one specific place.
6. Subjectivity of Post VS compound:  The result of a VS experiment is a ranking list of the whole screening database.  Taking top N results in very similar compounds – so some sort of post processing is usually carried out (e.g. clustering, but even a subjective manual filtering).  Difficult to reproduce across studies.
7. Prospective validation: the benchmarking of VS algorithms is done retrospectively.  Test on an active/decoy set for a particular target.  Only putative inactives.  Rarely validated in an external prospective context.
8. Drug-likeness: most VS experiments based on Lipinski Ro5 – not more than 5 hydrogen bond donors (nitrogen or oxygen atoms with one or more hydrogen atoms), not more than 10 hydrogen bond acceptors (nitrogen or oxygen atoms), a molecular mass less than 500 daltons, An octanol-water partition coefficient log P not greater than 5.  But these apply to oral bioavailability.  Lots of drugs fall out of this scope; intravenous drugs, antibiotic, peptidic drugs.  VS validated on Lipinski space molecules.
9. Diversity of benchmark library vs diversity of future , prospective vs runs: Library must fit the purpose of the experiment.  Most VS validation experiments are on commercially available libraries ~ small fraction of chemical space.  Type of screening library must be closely related to objective of VS campaign.   Results have to be transferable between runs. Validation on specific target family? If the goal is lead optimization combinatorial libraries are attractive.  Natural versus synthesized compounds; different chemical space.

Data design and content

1. Incomparability of benchmark sets: some datasets for docking studies, others for ligand based VS – incomparable methods.  In general 2D methods better than 3D (surprising!).  In 2D methods fingerprints outperform 3D methods.  Same datasets for validation of different methods – hard to reproduce any study otherwise. 
2. Limited comparability of performance metrics: Tag along on previous point; performance measurement used for different measurements should be the same. Mean EF risky because of ratio between actives to inactive molecules. ROC curves a problem because of early and late performance – use of BedROC (different importance to early and late stages of retrieved list of compounds). EF = (number of actives / number of expected) for a given % of the database
3. Hit rate in benchmark data sets; small libraries not good enough. Typical VS hit rates ~0.01% – 0.14%. Analogue bias; actives all look very similar to each other. Artificial enrichment; easy to tell between actives and decoys. Recent study found that for ligand based VS using no. of atoms gives half the VS performance.
4. Assay Comparability and Technology: Properly designed datasets such as MUV, use of similarity distributions to remove anything very similar to each other.  Remove problematic molecules like autofluoroscence .  MUV uses data from pubchem; different bioassays from different groups hence different quality.  Choices of targets; cutoffs; parameters; etc.  “Ideal VS benchmark deck will never happen.”
5. Bad molecules as actives: No real activity but either reactive or aggregating molecules in the assay which gives up a false positive; PAINs Pan assay interfering substances or frequent hitters. Small number of actives compared to inactives, false positives worse than false negatives.
6. Putative inactive compounds as decoys. The decoys are actually actives. 
7. Feature weights: LBVS based on a single query fails to identify important parts of the the molecule, e.g. benzamidine warhead in factor Xa inhibitors

Choice of Software

1. Interconverting chemical formats; errors or format incompatibilities.  Information lost or altered; or when using same format across different software (e.g. chirality, hybridization, and protonation states).
2. Molecule preparation; query molecules must be preprocessed exactly the same way as the structures in the database being screened to ensure consistency (e.g. partial charge calculation)
3. Feature definition: Specific rules which are sometimes left out of pharmacophoric definition. e.g. O, N in oxazole do not both behave as a HBA. Watch out for tautomers, protonation state, and chirality
4. Fingerprint selection and algorithmic implementation: different implementations of same fingerprint MACCS result in different fingerprints. Choice of descriptors; which ones to pick? Neighbourhood? Substructure?
5. Partial charges: Mesomeric effects; formal +1 charge spread over guanidine structure.
6. Single predictors versus ensembles: no single method works best in all cases. Consensus study; apply multiple methods and combine results.

Conformational sampling as well as ligand and target flexibility

1. Conformational coverage: four main parameters: (i) sampling algorithms and their specific parameters; (ii) strain energy cutoffs (iii) maximum number of conformations per molecule (iv) clustering to remove duplicates
2. Defining bioactive conformations: most ligands have never been co-crystallized with their primary targets and even fewer have been cocrystallized with counter targets. Same ligand might bind to different proteins in vastly different conformations. How easy is it to reproduce the cognate conformation? Also ligand changes shape upon binding. Minimum energy conformations are a common surrogate.
3. Comparing conformations: definitions of identity thresholds. 0.1 < rmsd < 0.5 excellent; 0.5 < rmsd < 1.0 good fit; 1.0 < rmsd < 1.5 acceptable; 1.5 < rmsd < 2.0 less acceptable; >2.0 not a fit in terms of biological terms. All atoms vs fragments RMSD makes direct comparison hard.
4. Size of conformational ensemble; trade off between computational cost and sampling breadth. Conformer generator may not generate bioactive conformation. How many conformations required to have bioactive one. Many bioactive conformations might exist.
5. Ligand flexibility – hard upper limit for no. of conformations. Conformer sizes depend mostly on number of rotatable bonds. Conformer generation tools don’t work well on some classes of molecules e.g. macrocycles
6. High energy conformations – high energy conformers (or physically unrealistic molecules; e.g. a cis secondary amide) detrimental to VS experiments. 3D pharmacophore searches sometimes result in matching strained structure; but 70% of ligands bind at strain energies below 3kcal/mol (stringent). 
7. Target flexibility – target flexibility – can do simple things like sidechain rotation, but nothing major like backbone flexibility. Sometimes docking to multiple structures snapshots resulting from molecular dynamics
8. Assumption of ligand overlap – lots of 3D shape based VS attempt to maximize the overlap between ligands – but based on X-ray structures this is not always the case (different ligands may occupy slight different regions of the binding pocket).
9. Missing positive controls – Strict cutoff stops you from retrieving positive controls in your Virtual Screening experiment. Selectivity (lower number of false postives)/ sensitivity (larger percentage of true positives) cutoff needs to be determined appropriately.

In conclusion, VS can be run by a monkey – but if that is the case expect bad results. Careful database preparation, judicious parameter choices, use of positive controls, and sensible compromises between the different goals one attempts to obtain are required. VS probabilistic game – careful planning and attention to detail increases probability of success.

In my Part II project, I have been comparing pairs of homologous alpha-helical transmembrane proteins from the PDB, looking specifically at how often kinks are conserved. Kinks are of particular interest in membrane protein structure as they can introduce flexibility for movement between conformations and are therefore important for function.

Kinks typically involve some disruption of the usual pattern of hydrogen bonding in the helix, so proline is frequently found at the location of helix kinks due to its lack of an N-H group. However, there are many kinks where proline is not present, and mutation of a proline residue to alanine at the site of a kink does not necessarily lead to removal of the kink. This led to the suggestion that all kinks are originally caused by prolines, but consequent changes to the surrounding structure mean that the kink can be retained even when a mutation occurs at a later point. Therefore I have particularly been looking at the incidence of prolines around kinks in the pairs of homologous helices.

In order to find kinks in helices I have been using two different programs: Kink Finder, written by Henry, another member of the group, and MC-HELAN. Kink Finder fits cylinders to a helix to measure angles at each residue, while MC-HELAN finds sections of ideal helix and labels any overlapping sections as kinks. These contrasting methods both indicate the extent of the kink by giving an angle between the axes fitted to sections of helix on either side. TM-align and MP-T have both been used to align the pairs of proteins and find homologous helices, and I am also looking at the multiple sequence alignments provided by MP-T to question whether proline can always be found in a homolog at the position of the kink, even if it is not present in the protein itself.

This work will hopefully give a better indication of when kinks should be predicted, particularly when using a program such as MEDELLER to predict membrane protein structure by homology modelling. If it were known in which circumstances a kink could be confidently predicted as conserved, the accuracy of the prediction of membrane regions could be further improved.