Adaptyv is a newish startup that sells high-throughput protein assays. The major innovations are (a) they tell you the price (a big innovation for biotech services!) (b) you only have to upload protein sequences, and you get results in a couple of weeks.

A typical Adaptyv workflow might look like the following:

• Design N protein binders for a target of interest (Adaptyv has 50-100 pre-specified targets)
• Submit your binder sequences to Adaptyv
• Adaptyv synthesizes DNA, then protein, using your sequences
• In ~3 weeks you get affinity measurements for each design at a cost of $149 per data-point

This is an exciting development since it decouples "design" and "evaluation" in a way that enables computation-only startups to get one more step towards a drug (or sensor, or tool). There are plenty of steps after this one, but it's still great progress!

The Adaptyv binder design competition

A couple of months ago, Adaptyv launched a binder design competition, where the goal was to design an EGFR binder. There was quite a lot of excitement about the competition on Twitter, and about 100 people ended up entering. At around the same time, Leash Bio launched a small molecule competition on Kaggle, so there was something in the air.

PAE and iPAE

For this competition, Adaptyv ranked designs based on the "PAE interaction" (iPAE) of the binder with EGFR.

PAE (Predicted Aligned Error) "indicates the expected positional error at residue x if the predicted and actual structures are aligned on residue y". iPAE is the average PAE for residues in the binder vs target. In other words, how accurate do we expect the relative positioning of binder and target to be? This is a metric that David Baker's lab seems to use quite a bit, at least for thresholding binders worth screening. It is straightforward to calculate using the PAE outputs from AlphaFold.

Unusually, compared to, say, a Kaggle competition, in this competition there are no held-out data that your model is evaluated on. Instead, if you can calculate iPAE, you know your expected position on the leaderboard before submitting.

The original paper Adaptyv reference is Improving de novo protein binder design with deep learning and the associated github repo has an implementation of iPAE that I use (and I assume the code Adaptyv use.)

Confusingly, there is also a metric called "iPAE" mentioned in the paper Systematic discovery of protein interaction interfaces using AlphaFold and experimental validation. It is different, but could actually be a more appropriate metric for binders?

At the end of last month (August 2024), there was a new Baker lab paper on Ras binders that also used iPAE, in combination with a few other metrics like pLDDT.

Experiments

A week or so after the competition ended, I found some time to try a few experiments.

Throughout these experiments, I include modal commands to run the relevant software. If you clone the biomodals repo it should just work(?)

iPAE vs Kd

The consensus seems to be that <10 represents a decent iPAE, but in order for iPAE to be useful, it should correlate with some physical measurement. As a small experiment, I took 55 PDB entries from PDBbind (out of ~100 binders that were <100 aas long, had an associated Kd, and only two chains), ran AlphaFold, calculated iPAE, and correlated this to the known Kd. I don't know that I really expected iPAE to correlate strongly with Kd, but it's pretty weak.

PDBbind Kd vs iPAE correlation

# download the PDBbind protein-protein dataset in a more convenient format and run AlphaFold on one example
wget https://gist.githubusercontent.com/hgbrian/413dbb33bd98d75cc5ee6054a9561c54/raw -O pdbbind_pp.tsv
tail -1 pdbbind_pp.tsv
wget https://www.rcsb.org/fasta/entry/6har/display -O 6har.fasta
echo ">6HAR\nYVDYKDDDDKEFEVCSEQAETGPCRACFSRWYFDVTEGKCAPFCYGGCGGNRNNFDTEEYCMAVCGSAIPRHHHHHHAAA:IVGGYTCEENSLPYQVSLNSGSHFCGGSLISEQWVVSAAHCYKTRIQVRLGEHNIKVLEGNEQFINAAKIIRHPKYNRDTLDNDIMLIKLSSPAVINARVSTISLPTAPPAAGTECLISGWGNTLSFGADYPDELKCLDAPVLTQAECKASYPGKITNSMFCVGFLEGGKDSCQRDAGGPVVCNGQLQGVVSWGHGCAWKNRPGVYTKVYNYVDWIKDTIAANS" > 6har_m.fasta
modal run modal_alphafold.py --input-fasta 6har_m.fasta --binder-len 80

Greedy search

This is about the simplest approach possible.

• Start with EGF (53 amino acids)
• Mask every amino acid, and have ESM propose the most likely amino acid
• Fold and calculate iPAE for the top 30 options
• Take the best scoring iPAE and iterate

Each round takes around 5-10 minutes and costs around $4 on an A10G on modal.

# predict one masked position in EGF using esm2
echo ">EGF\nNSDSECPLSHDGYCL<mask>DGVCMYIEALDKYACNCVVGYIGERCQYRDLKWWELR" > esm_masked.fasta
modal run modal_esm2_predict_masked.py --input-fasta esm_masked.fasta

# run AlphaFold on the EGF/EGFR complex and calculate iPAE
echo ">EGF\nNSDSECPLSHDGYCLHDGVCMYIEALDKYACNCVVGYIGERCQYRDLKWWELR:LEEKKVCQGTSNKLTQLGTFEDHFLSLQRMFNNCEVVLGNLEITYVQRNYDLSFLKTIQEVAGYVLIALNTVERIPLENLQIIRGNMYYENSYALAVLSNYDANKTGLKELPMRNLQEILHGAVRFSNNPALCNVESIQWRDIVSSDFLSNMSMDFQNHLGSCQKCDPSCPNGSCWGAGEENCQKLTKIICAQQCSGRCRGKSPSDCCHNQCAAGCTGPRESDCLVCRKFRDEATCKDTCPPLMLYNPTTYQMDVNPEGKYSFGATCVKKCPRNYVVTDHGSCVRACGADSYEMEEDGVRKCKKCEGPCRKVCNGIGIGEFKDSLSINATNIKHFKNCTSISGDLHILPVAFRGDSFTHTPPLDPQELDILKTVKEITGFLLIQAWPENRTDLHAFENLEIIRGRTKQHGQFSLAVVSLNITSLGLRSLKEISDGDVIISGNKNLCYANTINWKKLFGTSGQKTKIISNRGENSCKATGQVCHALCSPEGCWGPEPRDCVSCRNVSRGRECVDKCNLLEGEPREFVENSECIQCHPECLPQAMNITCTGRGPDNCIQCAHYIDGPHCVKTCPAGVMGENNTLVWKYADAGHVCHLCHPNCTYGCTGPGLEGCPTNGPKIPSI" > egf_01.fasta
modal run modal_alphafold.py --input-fasta egf_01.fasta --binder-len 53

One of the stipulations of the competition is that your design must be at least 10 amino acids different to any known binder, so you must run the loop above 10 or more times. Of course, there is no guarantee that there is a single amino acid change that will improve the score, so you can easily get stuck.

After 12 iteratations (at a cost of around $50 in Alphafold compute), the best score I got was 7.89, which would have been good enough to make the top 5. (I can't be sure, but I think my iPAE calculation is identical!) Still, this is really just brute-forcing EGF tweaks. I think the score was asymptoting, but there were also jumps in iPAE with certain substitutions, so who knows?

Unfortunately, though the spirit of the competition was to find novel binders, the way iPAE works means that the best scores are very likely to come from EGF-like sequences (or other sequences in AlphaFold's training set).

Adaptyv are attempting to mitigate this issue by (a) testing the top 200 and (b) taking the design process into account. It is a bit of an impossible situation, since the true wet lab evaluation happens only after the ranking step.

Bayesian optimization

Given an expensive black box like AlphaFold + iPAE, some samples, and a desire to find better samples, one appropriate method is Bayesian optimization.

Basically, this method allows you, in a principled way, to control how much "exploration" of new space is appropriate (looking for global minima) vs "exploitation" of variations on the current best solutions (optimizing local minima).

Bayesian optimization of a 1D function

The input to a Bayesian optimization is of course not amino acids, but numbers, so I thought reusing the ESM embeddings would be a decent, or at least convenient, idea here.

I tried both the Bayesian Optimization package and a numpyro Thompson sampling implementation. I saw some decent results at first (i.e., the first suggestions seemed reasonable and scored well), but I got stuck either proposing the same sequences over and over, or proposing sequences so diverged that testing them would be a waste of time. The total search space is gigantic, so testing random sequences will not help. I think probably the ESM embeddings were not helping me here, since there were a lot of near-zeros in there.

This is an interesting approach, and not too difficult to get started with, but I think it would work better with much deeper sampling of a smaller number of amino acids, or perhaps a cruder, less expensive, evaluation function.

ProteinMPNN

ProteinMPNN (now part of the LigandMPNN package), maps structure to sequence (i.e., the inverse of AlphaFold). For example, you can input an EGF PDB file, and it will return a sequence that should produce the same fold.

I found that for this task ProteinMPNN generally produced sequences with low confidence (as reported by ProteinMPNN), and as you'd expect, these resulted in low iPAEs. Some folds are difficult for ProteinMPNN, and I think EGF falls into this category. To run ProteinMPNN, I would recommend Simon Duerr's huggingface space, since it has a friendly interface and includes an AlphaFold validation step.

ProteinMPNN running on huggingface

# download a EGF/EGFR crytal structure and try to infer a new sequence that folds to chain C (EGF)
wget https://files.rcsb.org/download/1IVO.pdb
modal run modal_ligandmpnn.py --input-pdb 1IVO.pdb --extract-chains AC --params-str '--seed 1 --checkpoint_protein_mpnn "/LigandMPNN/model_params/proteinmpnn_v_48_020.pt" --chains_to_design "C" --save_stats 1 --batch_size 5 --number_of_batches 100'

RFdiffusion

RFdiffusion was the first protein diffusion method that showed really compelling results in generating de novo binders. I would recommend ColabDesign as a convenient interface to this and other protein design tools.

The input to RFdiffusion can be a protein fold to copy, or a target protein to bind to, and the output is a PDB file with the correct backbone co-ordinates, but with every amino acid labeled as Glycine. To turn this output into a sequence, this PDB file must then be fed into ProteinMPNN or similar. Finally, that ProteinMPNN output is typically folded with AlphaFold to see if the fold matches.

Although RFdiffusion massively enriches for binders over random peptides, you still have to screen many samples to find the really strong binders. So, it's probably optimistic to think that a few RFdiffusion-derived binders will show strong binding, even if you can somehow get a decent iPAE.

In my brief tests with RFdiffusion here, I could not generate anything that looked reasonable. I think in practice, the process of using RFdiffusion successfully is quite a bit more elaborate and heuristic-driven than anything I was going to attempt.

Figure 1 from De novo design of Ras isoform selective binders, showing multiple methods for running RFdiffusion

# Run RFdiffusion on the EGF/EGFR crystal structure, and diffuse a 50-mer binder against chain A (EGFR)
modal run modal_rfdiffusion.py --contigs='A:50' --pdb="1IVO"

Other things

A few other strategies I thought might be interesting:

• Search FoldSeek for folds similar to EGF. The idea here is that you might find a protein in another organism that wants to bind EGFR. I do find some interesting human-parasitic nematode proteins in here, but decided these were unlikely to be EGFR binders.
• Search NCBI for EGF-like sequences with blastp. You can find mouse, rat, chimp, etc. but nothing too interesting. The iPAEs are lower than human EGF, as expected.
• Search the patent literature for EGFR binders. I did find some antibody-based binders, but as expected for folds that AlphaFold cannot solve, the iPAE was low.
• Delete regions of the protein with low iPAE contributions to increase the average score. I really thought this would work for at least one or two amino acids, but it did not seem to. I did not do this comprehensively, but perhaps there are no truly redundant parts of this small binder?

Conclusion

All the top spots on the leaderboard went to Alex Naka, who helpfully detailed his methods in this thread. (A lot of this is similar to what I did above, including using modal!) Anthony Gitter also published an interesting thread on his attempts. I find these kinds of threads are very useful since they give a sense of the tools people are using in practice, including some I had never heard of, like pepmlm and Protrek.

Finally, I made a tree of the 200 designs that Adaptyv is screening (with iPAE <10 in green, <20 in orange, and >20 in red). All the top scoring sequences are EGF-like and cluster together. (Thanks to Andrew White for pointing me at the sequence data). We can look forward to seeing the wet lab results published in a couple of weeks.

Tree of Adaptyv binder designs

Molecular dynamics (MD) means simulating the forces acting on atoms. In drug discovery, MD usually means simulating protein–ligand interactions. This is clearly a crucial step in modern drug discovery, yet MD remains a pretty arcane corner of computational science.

This is a different problem to docking, where molecules are for the most part treated as rigid, and the problem is finding the right ligand orientation and the right pocket. Since in MD simulations the atoms can move, there are many more degrees of freedom, and so a lot more computation is required. For a great primer on this topic, see Molecular dynamics simulation for all (Hollingsworth, 2018).

What about deep learning?

Quantum chemical calculations, though accurate, are too computationally expensive to use for MD simulation. Instead, "force fields" are used, which enable computationally efficient calculation of the major forces. As universal function approximators, deep nets are potentially a good way to get closer to ground truth.

Analogously, fluid mechanics calculations are very computationally expensive, but deep nets appear to do a good job of approximating these expensive functions.

A deep net approximating Navier-Stokes

Recently, the SPICE dataset (Eastman, 2023) was published, which is a reference dataset that can be used to train deep nets for MD.

We describe the SPICE dataset, a new quantum chemistry dataset for training potentials relevant to simulating drug-like small molecules interacting with proteins. It contains over 1.1 million conformations for a diverse set of small molecules, dimers, dipeptides, and solvated amino acids.

This dataset has enabled new ML force fields like Espaloma (Takaba, 2023) and the recent Allegro paper (Musaelian, 2023), where they simulated a 44 million atom system of a HIV capsid. Interestingly, they scaled their system as high as 5120 A100's (which would cost $10k an hour to run!)

There are also hybrid ML/MM approaches (Rufa, 2020) based on the ANI2x ML force field (Devereux, 2020).

All of this work is very recent, and as I understand it, runs too slowly to replace regular force fields any time soon. Despite MD being a key step in drug development, only a small number of labs (e.g., Chodera lab) appear to work on OpenMM, OpenFF, and the other core technologies here.

Doing an MD simulation

I have only a couple of use-cases in mind:

• does this ligand bind this protein in a human cell?
• does this mutation affect ligand binding in a human cell?

Doing these MD simulations is tricky since a lot of background is expected of the user. There are many parameter choices to be made, and sensible options are not obvious. For example, you may need to choose force fields, ion concentrations, temperature, timesteps, and more.

By comparison, with AlphaFold you don't need to know how many recycles to run, or specify how the relaxation step works. You can just paste in a sequence and get a structure. As far as I can tell, there is no equivalent "default" for MD simulations.

A lot of MD tutorials I have found are geared towards simulating the trajectory of a system for inspection. However, with no specific numerical output, I don't know what to do with these results.

Choosing an approach

There are several MD tools out there for doing protein–ligand simulations, and calculating binding affinities:

• Schrodinger is the main player in computational drug discovery, and a mature suite of tools. It's not really suitable for me, since it's expensive, geared toward chemists, designed for interactive use over scripting, and not even necessarily cutting-edge.
• OpenEye also appears to be used a lot, and has close ties to open-source. Like Schrodinger, the tools are high quality, mostly interactive and designed for chemists.
• HTMD from Acellera is not open-source, but it has a nice quickstart and tutorials.
• GROMACS is open-source, actively maintained, and has tutorials, but is still a bit overwhelming with a lot of boilerplate.
• Amber, like GROMACS, has been around for decades. It gets plenty of use (e.g., AlphaFold uses it as a final "relaxation" step), but is not especially user-friendly.
• OpenMM seems to be where most of the open-source effort has been over the past five years or so, and is the de facto interface for a lot of the recent ML work in MD (e.g., Espaloma). A lot of tools are built on top or OpenMM:
  • yank is a tool for free energy binding calculations. Simulations are parameterized by a yaml file.
  • perses is also used for free energy calculation. It is pre-alpha software but under active development — e.g., this recent paper on protein–protein interaction. (Note, I will not claim to understand the differences between yank and perses!)
  • SEEKR2 is a tool that enables free energy calculation, among other things.
  • Making it rain is a suite of tools and colabs. It is a very well organized repo that guides you through running simulations on the cloud. For example, they include a friendly colab to run protein–ligand simulations. The authors did a great job and I'd recommend this repo broadly.
  • BAT, the Binding Affinity Tool, calculates binding affinity using MD (also see the related GHOAT).

OpenMM quickstart for protein simulation

Since I am not a chemist, I am really looking for a system with reasonable defaults for typical drug development scenarios. I found a nice repo by tdudgeon that appears to have the same goal. It uses OpenMM, and importantly has input from experts on parameters and settings. For example, I'm not sure I would have guessed you can multiply the mass of Hydrogen by 4.

This keeps their total mass constant while slowing down the fast motions of hydrogens. When combined with constraints (typically constraints=AllBonds), this often allows a further increase in integration step size.

I forked the repo, with the idea that I could keep the simulation parameters intact but change the interface a bit to make it focused on the problems I am interested in.

Calculating affinity

I am interested in calculating ligand–protein affinity (or binding free energy) — in other words, how well does the ligand bind the protein. There's a lot here I do not understand, but here is my basic understanding of how to calculate affinity:

• Using MD: This is the most accurate way to measure affinity, but the techniques are challenging. There are "end-point" approaches (e.g., MM/PBSA) and Free Energy Perturbation (FEP) / alchemical approaches. Alchemical free energy approaches are more accurate, and have been widely used for years. (I believe Schrodinger were the first to show accurate results (Wang, 2015).) Still, I found it difficult to figure out a straightforward way to do these calculations.
• Using a scoring function: This is how most docking programs like vina or gnina work. Docking requires a very fast, but precise, score to optimize.
• Using a deep net: Recently, several deep nets trained to predict affinity have been published. For example, HAC-Net is a CNN trained on PDBbind. This is a very direct way to estimate binding affinity, and should be accurate if there is enough training data.
  
  
  

The SQM/COSMO docking scoring function (Ajani, 2017)

Unfortunately, I do know know of a benchmark comparing all the above approaches, so I just tried out a few things.

Predicting cancer drug resistance

One interesting but tractable problem is figuring out if a mutation in a protein will affect ligand binding. For example, let's say we sequence a cancer genome, and see a mutation in a drug target, do we expect that drug will still bind?

There are many examples of resistance mutations evolving in cancer.

Common cancer resistance mutations (Hamid, 2020)

Experiments

BRAF V600E is a famous cancer target. Vemurafenib is a drug that targets V600E, and L505H is known to be a resistance mutation. There is a crystal structure of BRAF V600E bound to Vemurafenib (PDB:3OG7). Can I see any evidence of reduced binding of Vemurafenib if I introduce an L505H mutation?

PDB:3OG7, showing the distance between vemurafenib (cyan) and L505 (yellow)

I ran a simple simulation: starting with the crystal structure, introduce each possible mutation at position 505, allow the protein–ligand system to relax, and check to see if the new protein–ligand interactions are less favorable according to some measure of affinity.

I first used gnina's scoring function, which is fast and should be relatively precise (in order for gnina to work!) The rationale here was that the "obstruction" due to the resistance mutation would be detectable as the new atom positions of the amino acid and ligand would lead to a lower affinity.

Estimated affinity given mutations at position 505 in 3OG7

Nope. The resistance mutation has higher affinity (realistically, there are no distinguishable differences for any mutation).

We also know that MEK1 V215E acts as a resistance mutation to PD0325901, and the PDB has a crystal structure of MEK1 bound to PD0325901 (PDB:70MX).

Estimated affinity given mutations at position 215 in 70MX

Again, I can't detect any difference in affinity due to the resistance mutation.

HAC-Net

I also tried a deep-learning based affinity calculator, HAC-Net. HAC-Net has a nice colab and is relatively easy to run Dockerized.

The HAC-Net colab gives me a pKd of 8.873 for 3OG7 (wild-type)

Estimated pKd given mutations at position 505 in 3OG7 using HAC-Net

I still see no difference in affinity with HAC-Net.

Each of these simulations (relaxing a protein–ligand system with solvent present) took a few minutes on a single CPU. If I wanted to simulate a full trajectory, which could be 50 nanoseconds or longer, it would take hundreds or thousands of times as long.

Conclusions

On the one hand, I can run state-of-the-art MD simulations pretty easily with this system. On the other hand, I could not discriminate cancer resistance mutations from neutral mutations.

There are several possible reasons. Most likely, the short relaxation I am doing is insufficient and I need longer simulations. The simulations may also be insufficiently accurate, either intrinsically or due to poor parameterization. It's difficult to feel confident in these kinds of simulations, since there's no simple way to verify anything.

If anyone knows of any obvious fix for the approach here, let me know! Probably the next thing I would try would be adapting the Making It Rain tools, which include (non-alchemical) free energy calculation. For some reason the Making It Rain colab specifies "This notebook is NOT a standard protocol for MD simulations! It is just simple MD pipeline illustrating each step of a simulation protocol." which begs the questions, why not and where is such a notebook?

I do think that enabling anyone to run such simulations — specifically, with default parameters blessed by experts for common use-cases — would be a very good thing.

There are already several cancer drug selection companies like Oncobox, so maybe there should be a company doing this kind of MD for predicting cancer resistance. Maybe there is and I just have not heard of it?

Addendum: modal labs

I have been experimenting with modal labs for running code like this, where there are very specific environment requirements (i.e., painful to install libraries) and heavy CPU/GPU load. Previously, I would have used Docker, which is fundamentally awkward, and still requires manually provisioning compute. Modal can be a joy to use and I'll probably write up a separate blogpost on it.

To do your own simulation (bearing in mind all the failed experiments above!), you can either use my MD_protein_ligand colab or if you have a modal account, clone the MD_protein_ligand repo and run

mkdir -p input && modal run run_simulation_modal.py --pdb-id 3OG7 --ligand-id 032 --ligand-chain A

This basic simulation (including solvent) should cost around 10c on modal. That means we could relax all 5000 protein–ligand complexes in the PDB for around $500, perhaps in just a day or two (depending on how many machines modal allows in parallel). I'm not sure if there's any point to that, but pretty cool that things can scale up so easily these days!

Last year I wrote a post about computational tools for drug development. Since then quite a lot has happened, especially the appearance of several generative models based on equivariant neural networks for drug design. This article is a sort of update to that post, and also a collection of colabs I have found or developed over the past year or so that can be stitched together to design drugs.

The tools listed here are focused on the earliest stages of drug development. Specifically:

1. Generate a new molecule or use a virtual screen to find one
2. Evaluate the molecule's potential as a drug
3. Synthesize or purchase the molecule

1a. Generate a new molecule

Pocket2Mol: Protein structure → small molecules

Pocket2Mol is one of a new crop of generative models that start with a binding pocket and generate molecules that fit in the pocket. VD-Gen (with animation below) is similar to Pocket2Mol but since it has no code available so I cannot tell if it works.

I created a Pocket2Mol colab that enables easy molecule generation. The input is a PDB file and 3D coordinates to search around. The centroid of a bound ligand in the PDB file can serve as the coordinates. The output is a list of generated molecules. I rank the generated molecules with gnina, a fast, easy to install, and relatively accurate way to measure binding affinity.

ColabDesign: Protein structure → peptide binders

ColabDesign is an extremely impressive project that democratizes a lot of the recent work in generative protein models. As you may guess from the name, ColabDesign — and sister project ColabFold — allow anyone to run them via colab.

For example, even the recent state-of-the-art RFDiffusion algorithm from the Baker lab has been incorporated and made available as an RFDiffusion colab.

ColabDesign can generate proteins that conform to a specific shape or reference protein backbone (structure → sequence, i.e., AlphaFold in reverse). It can also generate peptides that can bind to a specific protein.

The same group also developed AfCycDesign, which uses a nice trick to get Alphafold to fold cyclic peptides — an increasingly important drug type, and often an alternative to antibodies. There is an AfCycDesign colab too!

In this animation it is attempting to generate a cyclic peptide covid spike protein binder, as published recently.

1b. Run a virtual screen

gnina: PDB structure + ligand → posed ligand + binding affinity

gnina is the deep-learning–based successor to the extremely popular smina, itself an AutoDock Vina fork. There is a minimal implementation available as a gnina colab. gnina is a successor to smina and appears to be strictly better. gnina can be used for virtual screening and performs very well at that.

PointVS is a new, deep-learning approach. Like many recent methods, it uses an EGNN (equivariant graph neural network). PointVS's performance is impressive, comparable to gnina, but it's unclear which method runs faster. Uni-Dock is a new GPU-accelerated Autodock Vina that claims impressive speed but unfortunately there is no open code.

PointVS and gnina perform comparably

DiffDock: protein structure + ligand → posed ligand

DiffDock is a generative diffusion model with impressive performance. Given a protein and a ligand, DiffDock can try to find the best pose for that ligand. I created a DiffDock colab that runs DiffDock and again ranks the results using gnina.

DiffDock is a pose prediction method and is not designed to do virtual screening per se. It returns a "confidence" score that correlates with smina/gnina affinity.

DiffDock confidence and smina/gnina binding affinities correlate fairly well

As for which algorithm produces superior results, it's unclear, but according to a recent review, gnina comes out ahead.

2. Evaluate molecule properties

SMILES property predictor: molecule → predicted properties

A molecule the binds to a protein has one property necessary to become a drug, but that is far from the only thing that matters. We can also predict a molecule's intrinsic properties using machine learning.

I created a Smiles to properties colab based on chemprop, soltrannet, and SMILES2Caption. Most of the training data is from MoleculeNet. For some reason there was no pre-trained network available so I had to retrain from scratch.

The graphs generated show the predicted properties of your SMILES as compared to FDA-approved drugs (grey bounds). In the example below, acetaminophen appears as the most toxic of the group, as expected.

SwissADME also provides an excellent property prediction service, though anecdotally analyzing proprietary molecules on someone else's server is not commonly done in drug development.

Espaloma Charge: molecule → charge

Espaloma Charge is a standalone part of Chodera lab's Espaloma system. I include it here as they make available a very simple Espaloma Charge colab that will assign charges for your molecule.

HAC-Net: Protein structure + posed ligand → binding affinity

HAC-Net is a new deep learning method specifically designed to estimate protein–ligand binding affinity.

The HAC-Net colab takes as input a protein and posed ligand and returns a dissociation constant (pKd) and a nice pymol image. This method may or may not be more accurate than gnina — I have not attempted to benchmark.

3. Synthesize or purchase molecule

Postera Manifold: SMILES → molecule or building blocks

Postera Manifold is a machine learning tool that helps figure out how to synthesize small molecules. It can supply the building blocks and reactions necessary for synthesis, or there is a "Have PostEra make it for you" button (at least according to their documentation, but I could not find it!) However, usually synthesis of a new small molecule is a custom, CRO-driven process, that could be expensive and take a long time depending on the molecule's complexity.

Biomatik: SMILES → purchasable molecule

For peptides — which are linear, modular molecules — synthesis is much easier. A service like Biomatik (or GenScript, or many others) can supply peptides for as little as $80 per. They can use proteinogenic or non-proteinogenic amino acids, and can cyclize the peptide too.

Small World: SMILES → purchasable molecule

Small world allows you to search for molecules in the largest compound databases: ZINC, Mcule, Enamine REAL... There is also a nice, unofficial Small World Python API.

If you find a close match to your designed molecule, you could evaluaate it to see if it works as a binder, or use it as a starting point for modification.

Orforglipron

As an example of using these tools, I'll use a new GLP-1 agonist, orforglipron, that showed promising results in a phase 2 study published in NEJM just this week.

First, I can generate an image of the protein GLP-1 with a different bound ligand using my pdb2png colab. (PDB:7S15, "Pfizer small molecule bound to GLP-1").

I can look at the molecule and its charges using the Espaloma Charge colab. (Nice image, but to be honest this doesn't mean much to me!)

I can get the molecule's SMILES from pubchem, and search for it using Small World.

orforglipron in Small World: dark green is a match; light green is a mismatch

A little surprisingly, I find a matching molecule that is only two edits from orforglipron in the ZINC20 "for sale" set. It is labeled as orforglipron in MCE and interestingly the text says it was "extracted from a patent". It's unclear to me why it's not an exact match. Most sources charge $2k for 1mg, but there is one source that charges $1k for 5mg.

I can evaluate its properties with SwissADME, and it's interesting to note how far outside the ideal orforglipron is, in terms of size (883 Da!) and predicted solubility.

I can use the DiffDock colab to see how well it might bind to GLP-1

The best binding score I could get after three attempts was a DiffDock confidence of -2 (very unconfident) and a gnina affinity of -6 (poor-to-mediocre affinity). Because this molecule is so large, there are many possible conformers and it may be impossible to adequately sample them adequately. It does at least bind to the same place as the Pfizer small molecule.

When I feed this top pose from DiffDock into the HAC-Net colab, I get a pKd of 11, which is very high affinity.

The tools for drug design are getting very powerful, and colab is a fantastic way to make them widely available and easy to run. Even when I have a command-line equivalent set up, I often find running the colab to be quicker and easier. In theory, Docker provides the same kind of advantages, but it's never quite so easy, as you still have to provision compute, disk, etc. Having easy access to A100's (on the extremely affordable Colab Pro) is not bad either.