I wrote about querying research papers with LLMs last February. Since then things have progressed a lot, with significant new LLM results seemingly every week. However, I think it's fair to say there is still no turnkey way to search PubMed, load a corpus of PDFs into an LLM, query, and get good quality results. (If I've missed something great, let me know...)

There are several options:

• The tools I used in my previous article, Paper QA and LlamaIndex, can work well if you already have a corpus of papers. Indexing papers with GPT-4 can get expensive though.

• New RAG (Retrieval Augmented Generation) tools like RAGatouille (which uses the new ColBERT model, now built into LangChain)

• All-in-one local LLM UIs like gpt4all

• Web-based scientific services like elicit.org and scite.ai. I can't tell exactly which papers these tools have access to. scite.ai seems to be focused on "citation statements" (what?) and elicit.org on abstracts. I imagine it is difficult for these services to get access to all full text articles.

As of January 2024, GPT-4 still dominates all open-source models in terms of reasoning performance, but it is expensive for larger projects, so for now I am still exploring and comparing approaches. (For various reasons, I could get neither RAGatouille nor gpt4all to run on my mac!)

Long context

One really useful recent development in LLMs is long context (>100,000 tokens, or >100,000 words, approximately). A few models now have long context: GPT-4 Turbo, Claude, and some open-source models like Anima. By comparision, GPT-3 launched with a context length of only 2,048 tokens — barely enough to hold a few abstracts.

If you are trying to use an LLM to read and analyze documents, you have two main choices: you can fine-tune the model with additional training on a comparatively set of documents, or you can keep the default model and provide the document text directly in context. Perhaps surprisingly, it seems like the latter actually works fine, often better! Long context means you can directly load hundreds of pages of text into the model's memory. Intuitively, it does seem cleaner to keep the reasoning engine (LLM) separate from the domain knowledge (papers).

A recent post on twitter by @DeanCarignan

Notably, smaller contexts can work well with RAG, I assume because RAG pulls out only relevant passages, resulting in a better signal-to-noise ratio.

Text sources

For biomedical research, there are a few sources of text:

• abstracts, available free online via PubMed and Google Scholar
• full text papers (PDFs, HTML or JSON), some of which are available free via PubMed Central (more and more thanks to the NIH Public Access Policy)
• patents, available free via USPTO or Google Patents
• internet forums and discords. I have not actually looked into using this, but there is definitely a lot of high quality discourse on some discords.

BioC

BioC is a very useful format that is new to me. Essentially it's an NCBI-generated structured JSON file with annotated sections for "intro", "methods", etc.

Example BioC JSON

One unfortunate caveat is that I found not all PMIDs are available in BioC, even if the paper is available, so I have to use PMCIDs (PubMed Central publications) instead. No idea why, but I did email them about it.

Script

So here is a simple script that uses the PubMed and OpenAI APIs to rank a set of genes for a given "role" (e.g., ["MAP2K1", "F5"] and "cancer driver"). It searches PubMed, reads the abstracts to figure out if the gene has the specified role, summarizes the results in a JSON object, including a score from 0-10 that can be used for ranking. I found this script worked pretty well for me for ranking a few hundred genes — a task that would have been way too laborious for me to do — and it's a nice short script. The LLM prompt is far from optimized though, so buyer beware!

The long context model (GPT-4 Turbo) means I can simply include all the abstracts (I include the first 30) in the text of my query, and not worry about fine-tuning or RAG. The equivalent Paper QA script would have to additionally include an indexing and retrieval step, but could include the full text where available.

It feels strange to write these LLM scripts because most of the "code" is in English. The prompt engineering part took inspiration from Jeremy Howard from FastAI and this article from finedataproducts, which covers a lot of the same topics as this one, but in more depth.

Note that GPT-4 Turbo is still quite expensive, so although the results are pretty good, if you were attempting to search for all 20,000 genes in the human proteome, it would get expensive. Also, since you pay per token, the longer the context you provide, the more expensive it gets — proper RAG is probably a better idea for bigger projects!

Example results

Here are some results from running the script for the role "cancer driver". Note, the Summary field is an opportunity for the LLM to think "step-by-step". Surprisingly, the phrase "cancer driver" is not that commonly found with specific genes so some genes return no result.

Here are some other example use-cases for this script; in my experience, pretty common types of tasks in biotech:

• Analyze a set of genes and rank them by their importance in cancer metabolism, immune response, lupus, etc.
• Analyze a set of compounds and rank them by their toxicity, potency, etc.
• Analyze a set of bacterial strains and rank them by their presence in human microbiome, potential use as probiotics, etc.

With these kinds of queries, you'll likely get reasonable results, with the obvious candidates at the top (and perhaps a few surprises), comparable to what you might get from an inexperienced intern with unlimited time.

import json
import re
import requests

from pathlib import Path
from textwrap import dedent
from time import sleep
from urllib.parse import quote

from openai import OpenAI
client = OpenAI()

# Requires OpenAI API Key: https://openai.com/blog/openai-api
# os.environ["OPENAI_API_KEY"] = "sk-xxx"

SEARCH_PUBMED_URL = "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=pubmed&term={params}&retmax={max_abstracts}{api_key_param}"
GET_ABSTRACT_URL = "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=pubmed&id={pmid}&retmode=text&rettype=abstract{api_key_param}"
# https://www.ncbi.nlm.nih.gov/research/bionlp/APIs/BioC-PubMed/
BIOC_PMCID_URL = "https://www.ncbi.nlm.nih.gov/research/bionlp/RESTful/pmcoa.cgi/BioC_json/{pmcid}/unicode"

DEFAULT_MAX_ABSTRACTS = 30

# NCBI recommends that users post no more than three URL requests per second and limit large jobs
# Failure to comply with this policy may result in an IP address being blocked from accessing NCBI.
# Add an API key if you want to download faster
NCBI_API_KEY = None

def remove_ref_blank_entries_and_offsets(obj):
    """
    Recursively traverse through the JSON object (dict or list) and:
    1. Remove any node if it or any of its nested structures contains a dict with 'section_type': 'REF'.
    2. Remove any key-value pairs where the value is an empty list or empty dictionary.
    3. Remove any key-value pairs where the key is 'offset'.
    """
    if isinstance(obj, dict):
        # Check if any value of this dict is a nested dict with 'section_type': 'REF'
        if any(isinstance(v, dict) and v.get('section_type') == 'REF' for v in obj.values()):
            return None
        else:
            # No nested dict with 'section_type': 'REF', recursively process each key-value pair
            return {k: remove_ref_blank_entries_and_offsets(v) for k, v in obj.items() if k != 'offset' and
                       remove_ref_blank_entries_and_offsets(v) is not None and v != [] and v != {}}
    elif isinstance(obj, list):
        # Recursively process each item in the list
        return [remove_ref_blank_entries_and_offsets(item) for item in obj
                if remove_ref_blank_entries_and_offsets(item) is not None and item != [] and item != {}]
    else:
        # Return the item as is if it's not a dict or list
        return obj


def main(genes, role, max_abstracts=DEFAULT_MAX_ABSTRACTS, ncbi_api_key=NCBI_API_KEY):
    """Download abstracts and fulltext from pubmed, generate a prompt, query OpenAI.
    """
    Path("downloads").mkdir(exist_ok=True)
    Path("gpt_output").mkdir(exist_ok=True)

    api_key_param = "&api_key={ncbi_api_key}" if ncbi_api_key is not None else ""

    for gene in genes:
        abstracts = []
        fulltexts = []

        params = quote(f'({gene}) AND ("{role}")')

        pmtxt = requests.get(SEARCH_PUBMED_URL.format(params=params, max_abstracts=max_abstracts, api_key_param=api_key_param)).text
        pmids = re.findall("<Id>(\d+)</Id>", pmtxt)
        sleep(0.3)

        for pmid in pmids:
            if Path(f"downloads/abstract.pmid_{pmid}.txt").exists():
                abstracts.append(open(f"downloads/abstract.pmid_{pmid}.txt").read())
                if Path(f"downloads/fulltext.pmid_{pmid}.json").exists():
                    fulltexts.append(json.load(open(f"downloads/fulltext.pmid_{pmid}.json")))
                continue

            abstract = requests.get(GET_ABSTRACT_URL.format(pmid=pmid, api_key_param=api_key_param)).text
            open(f"downloads/abstract.pmid_{pmid}.txt", 'w').write(abstract)
            abstracts.append(abstract)
            sleep(0.3)

            pmcid = re.findall("^PMCID:\s+(PMC\S+)$", abstract, re.MULTILINE)
            assert len(pmcid) <= 1, "there should be max one PMCID per paper"
            pmcid = pmcid[0] if len(pmcid) == 1 else None

            if pmcid is not None:
                fulltext_request = requests.get(BIOC_PMCID_URL.format(pmcid=pmcid))
                if fulltext_request.status_code == 200:
                    fulltext_json = remove_ref_blank_entries_and_offsets(fulltext_request.json())
                    json.dump(fulltext_json, open(f"downloads/fulltext.pmid_{pmid}.json", 'w'), indent=2)
                    fulltexts.append(fulltext_json)
                sleep(0.3)

        # If there are only a couple of papers then we can include the fulltext
        # 250k characters should be well below 100k tokens; 500k characters is too many
        abstract_len = len("\n".join(abstracts))
        fulltext_len = len("\n".join(json.dumps(j) for j in fulltexts))
        if abstract_len + fulltext_len < 250_000:
            abstracts_txt = "\n".join(abstracts) + "\n" + "\n".join(json.dumps(j) for j in fulltexts)
        else:
            abstracts_txt = "\n".join(abstracts)

        gpt_out = f"gpt_output/{gene}_{role.replace(' ','_')}.txt"
        if len(abstracts_txt) > 1_000 and not Path(gpt_out).exists():
            completion = client.chat.completions.create(
                model="gpt-4-1106-preview",
                messages=[
                {"role": "system",
                "content": dedent("""
                You are an autoregressive language model that has been fine-tuned with instruction-tuning and RLHF.
                You carefully provide accurate, factual, thoughtful, nuanced answers, and are brilliant at reasoning. 
                If you think there might not be a correct answer, you say so. 
                Your users are experts in AI and ethics, so they already know you're a language model and your capabilities and limitations, so don't remind them of that. 
                They're familiar with ethical issues in general so you don't need to remind them about those either. 
                Don't be verbose in your answers, but do provide details and examples where it might help the explanation.

                Your users are also experts in science, and especially biology, medicine, statistics. 
                Do NOT add any details about how science or research works, tell me to ask my doctor or consult with a health professional.
                Do NOT add any details that such an expert would already know.
                """)},
                {"role": "user", 
                "content": dedent(f"""
                Help me research whether this gene has the role of {role}: {gene}. 
                Here are the top abstracts from pubmed, and any full text that will fit in context:

                {abstracts_txt}

                Read and synthesize all the evidence from these papers.
                Prefer evidence from good journals and highly cited papers, and papers that are recent.
                If there is one standout result in a top journal like Nature, Science or Cell, focus on that.
                There should usually be one primary result and most of the evidence should depend on that.

                Include a score out of 10, where 1 would be a single paper with weak evidence,
                5 would be a mid-tier journal with a single paper with believable evidence,
                and 10 would be multiple Nature, Science or Cell papers with very strong evidence.
                A score above 7 is exceptional, implying at least the result has been replicated and is trustworthy.

                Make sure the results are valid json with fields exactly as in the examples below.
                Do NOT change the fields, change the order, or add any other fields.
                The "Gene" json entry MUST match the query gene. Here, {gene}.                
                Include AT LEAST ONE reference for each entry on evidence (unless there is no evidence), e.g., (Smith et al., Nature, 2021). ALWAYS include the journal name.
                To help you think step-by-step about your response, use the FIRST "Summary" entry to summarize the question
                and the relevant available evidence in your own words, in one hundred words or fewer.
                Do NOT output markdown, just raw json. The first character should be {{ and the last character should be }}.

                Here are some examples of json you might output:

                {{
                "Summary": "ABC1 has no evidence of {role}. The experiments are not convincing. There is some evidence for ABC2, which may indicate something."; // string
                "Gene": "ABC1", // string
                "Cancer type": "", // string
                "Evidence": "No evidence as {role}", // string
                "Score": 0 // integer from 0-10
                }}

                {{
                "Summary": "There is some evidence for ABC2, which may indicate something. The evidence is medium, e.g. a paper in a mid-tier journal" // string
                "Gene": "DEF2", // string
                "Cancer type": "Bladder cancer", // string
                "Evidence": "A CRISPR screen identified DEF2 as {role} in bladder cancer (Jones et al., Nature Cancer, 2001)", // string
                "Score": 5 // integer from 0-10
                }}

                {{
                "Summary": "GHI3 has been cited as {role} many times. The experiments are convincing. The evidence is strong: a recent Nature paper and a Science paper" // string
                "Gene": "GHI3" // string
                "Cancer type": "Colorectal cancer" // string
                "Evidence": "A cell-based assay using a small molecule inhibitor identified GHI3 as {role} in colorectal cancer (Smith et al., Science, 2022, Thomson et al., Nature, 2019)", // string
                "Score": 8 // integer from 0-10
                }}
                """)
                }
            ]
            )
            open(gpt_out, 'w').write(completion.choices[0].message.content)
        else:
            # just log null results
            open(gpt_out, 'w').write("")

if __name__ == "__main__":
    main(genes = ["MAP2K1", "F5"], role = "cancer driver")

For a long time, the dream has been to be able to test code on your laptop and transparently scale it to infinite compute on the cloud. There are many, many tools that can help you do this, but modal, a new startup, comes closer to a seamless experience than anything I've used before.

There is a really nice quickstart, and they even include a generous $30/month to help get your feet wet.

pip install modal
python3 -m modal setup
git clone https://github.com/hgbrian/biomodals
cd biomodals
modal run modal_omegafold.py --input-fasta modal_in/omegafold/insulin.fasta

OmegaFold

Here's an example of how to use modal to run OmegaFold. OmegaFold is an AlphaFold/ESMFold/ColabFold-like algorithm. It is much easier to run than AlphaFold (which needs 2TB+ of reference data!) and it performs well according to a recent benchmark by 310.ai.

import glob
from subprocess import run
from pathlib import Path
from modal import Image, Mount, Stub

FORCE_BUILD = False
MODAL_IN = "./modal_in"
MODAL_OUT = "./modal_out"

stub = Stub()

image = (Image
         .debian_slim()
         .apt_install("git")
         .pip_install("git+https://github.com/HeliXonProtein/OmegaFold.git", force_build=FORCE_BUILD)
        )

@stub.function(image=image, gpu="T4", timeout=600,
               mounts=[Mount.from_local_dir(MODAL_IN, remote_path="/in")])
def omegafold(input_fasta:str) -> list[tuple[str, str]]:
    input_fasta = Path(input_fasta)
    assert input_fasta.parent.resolve() == Path(MODAL_IN).resolve(), f"wrong input_fasta dir {input_fasta.parent}"
    assert input_fasta.suffix in (".faa", ".fasta"), f"not fasta file {input_fasta}"

    run(["mkdir", "-p", MODAL_OUT], check=True)
    run(["omegafold", "--model", "2", f"/in/{input_fasta.name}", MODAL_OUT], check=True)

    return [(pdb_file, open(pdb_file, "rb").read())
            for pdb_file in glob.glob(f"{MODAL_OUT}/**/*.pdb", recursive=True)]

@stub.local_entrypoint()
def main(input_fasta):
    outputs = omegafold.remote(input_fasta)

    for (out_file, out_content) in outputs:
        Path(out_file).parent.mkdir(parents=True, exist_ok=True)
        if out_content:
            with open(out_file, 'wb') as out:
                out.write(out_content)

Hopefully the code is relatively self-explanatory. It's just Python code with code to set up the docker image, and some decorators to tell modal how to run the code on the cloud.

There are only really two important lines: installing OmegaFold with pip:

.pip_install("git+https://github.com/HeliXonProtein/OmegaFold.git", force_build=FORCE_BUILD)

and running OmegaFold:

run(["omegafold", "--model", "2", f"/in/{input_fasta.name}", MODAL_OUT], check=True)

The rest of the code could be left unchanged and reused for many bioinformatics tools. For example, to run minimap2 I would just add:

.run_commands("git clone https://github.com/lh3/minimap2 && cd minimap2 && make")

Finally, to run the code:

modal run modal_omegafold.py --input-fasta modal_in/omegafold/insulin.fasta

Outputs

Modal has extremely nice, almost real-time logging and billing.

CPU and GPU usage for an OmegaFold run. Note how it tracks fractional CPU/GPU use.

Billing information for the same run, split into CPU, GPU, RAM.

This OmegaFold run cost me 3c and took about 3 minutes. So if I wanted to run OmegaFold on a thousand proteins I could probably do the whole thing for ~$100 in minutes. (I'm not sure since my test protein is short, and I don't know how many would run in parallel.) In theory, someone reading this article could go from never having heard of modal or OmegaFold to folding thousands proteins in well under an hour. That is much faster than any alternative I can think of.

Modal vs Docker

Why not just use Docker? Docker has been around forever and you can just make a container and run it the same way! Modal even uses Docker under the hood!

There are significant differences:

• Dockerfiles are weird, have their own awful syntax, and are hard to debug;
• you have to create and manage your own images;
• you have to manage running your containers and collecting the output.

I have a direct equivalent to the above OmegaFold modal script that uses Docker, and it includes:

• a Dockerfile (30 lines);
• a Python script (80 lines, including copying files to and from buckets);
• a simple bash script to build and push the image (5 lines);
• and finally a script to run the Dockerfile on GCP (40 lines, specifying machine types, GPUs, RAM, etc).

Also, Docker can be slow to initialize, at least when I run my Docker containers on GCP. Modal has some impressive optimizations and cacheing (which I do not understand). I find I am rarely waiting more than a minute or two for scripts to start.

Modal vs snakemake, nextflow, etc

There is some overlap here, but in general the audience is different. Nextflow Tower may be a sort-of competitor, I have not tried to use it.

The advantages of these workflow systems over modal / Python scripts are mainly:

• you can separate your process into many atomic steps;
• you can parameterize different runs and keep track of parameters;
• you can cache steps and restart after a failure (a major advantage of e.g., redun).

However, for many common tasks like protein folding (OmegaFold, ColabFold), genome assembly (SPAdes), read mapping (minimap2, bwa), there is one major step — executing a binary — and it's unclear if you need to package the process in a workflow.

Pricing

Modal is priced per second so you don't pay for more than you use.

Modal runs (transparently) on AWS, GCP, and Oracle. Another blogpost claims that Modal adds around a 12% margin. However, since Modal charges per cycle, I think it's possible you could end up saving quite a bit of money? If you run your own Docker image on a VM, you may end up paying a lot for idle time; for example, if your GPU is utilized for a fraction of the time (as with AlphaFold). It's very unclear to me how modal makes this work (or if I am misunderstanding something), but it's really nice to not have to worry about maximizing utilization.

One challenge with modal — and all cloud compute — for bioinformatics is having to push large files around (e.g., TB of sequencing reads). If you want to do that efficiently you may have to look into the modal's Volumes feature so larger files only get uploaded once.

Conclusion

Not too long ago, I felt I had to do almost everything on a linux server, since my laptop did not have enough RAM, or any kind of GPU. Now with the M-series MacBook Pros, you get a GPU and as much RAM as high end server (128GB!) I still need access to hundreds of cores and powerful GPUs, but only for defined jobs.

I am pretty excited about modal's vision of "serverless" compute. I think it's a big step forward in how to manage compute. Their many examples are illustrative. Never having to think about VMs or even having to choose a cloud is a big deal, especially these days when GPUs are so hard to find (rumor has it Oracle has spare A100's, etc!) Although it's an early startup, there is almost no lock-in with modal since everything is just decorated Python.

I made a basic biomodals repo and added some examples.

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!