Well-Prepared Proteins and How to Use Them
better protein preparation through ML; Gnina & covalent docking; MM/GBSA refinement; fast MD and more trajectory analysis; synthetic RBFE intermediates; resubmitting FEP graphs without rerunning legs
We’re excited to share a variety of updates to Rowan’s structure-based drug discovery toolkit today. As we work to support more and more biotech/pharma teams, we’ve been finding lots of ways to increase our platform’s capabilities and depth.
(To help keep things manageable, this newsletter is organized in pipeline order, starting with protein preparation and ending with FEP.)
Protein Preparation and Protein pKa
Getting a crystal structure ready for simulation (whether docking, molecular dynamics, or FEP) is no easy task. Protein structures deposited in the PDB are often missing residues from flexible sections of a polymer chain as well as atoms from flexible side chains (not to mention hydrogens). The limitations of crystallographic methods mean there’s always work that needs to be done before a protein is simulation ready.
This is the problem of protein preparation; we previously offered a solution to this problem with PDBFixer, a set of tools built to get proteins ready for molecular dynamics (MD) simulations.
Today, we’re releasing a new workflow aimed at improving this process. At a high level, the workflow does a three-step process:
First, we add missing residues and any other missing heavy atoms. To add missing atoms, we support PDBFixer (as previously) as well as Boltz-2. When Boltz-2 is selected, we refold all the polymer chains using the input structure as a template. We then run an MD relaxation restrained to the input structure to enforce stronger template adherence. Co-folding models like Boltz are great at filling in missing flexible regions with plausible guesses, enabling preparation of difficult systems that PDBFixer struggles to prepare.
Next, we add hydrogens. To predict the per-residue protonation states, we support OpenMM’s Modeller (as previously) as well as the protonate_utils package from Pat Walters, which builds on the Hydride package. We’re interested in adding more protein-pKa methods in the near future.
Finally, we optimize hydrogen positions, iteratively flipping polymer side chains to avoid any donor–donor and acceptor–acceptor hydrogen-bond clashes. The flipping is followed by a quick hydrogen-only MM relaxation to ensure that the workflow returns a fully MD-ready system.
Gnina and Covalent Docking
We’ve also added a new docking engine, Gnina. Unlike the empirical physics-based scoring functions used in most docking programs (like Vina), Gnina scores poses with a convolutional neural network trained on experimental protein-ligand structures. In practice, this makes Gnina somewhat slower and more accurate than other docking programs: a variety of benchmarks have shown that Gnina outperforms Vina and other docking methods at pose retrieval and binder enrichment (see e.g. 1, 2, 3, 4).
Gnina also supports covalent docking, which can be used to model the docking of compounds that form a covalent bond with a protein. In practice, this means a user must additionally specify the atom in the ligand (the warhead atom) and the atom on the targeted residue to run covalent docking—Gnina will link these atoms together to generate the initial docked pose before running additional conformation-exploration steps. (You can read more about how this works in the Gnina 1.3 publication.)
Gnina can be run through Rowan’s docking workflows by selecting “Gnina” in the docking-settings dropdown. For covalent docking, input the atom indices of the reactive protein and ligand atoms; for non-covalent docking, leave these fields blank.
The “covalent docking” process modeled by Gnina and other non-covalent docking methods is not fully physical and produces a single hybrid score, rather than the multistep process of binding and bond formation undergone by real covalent inhibitors. We’re working on a QM-based method to simulate the individual mechanistic steps of covalent drug binding, which we expect to be more physical and to better capture warhead SAR than covalent docking; if this is of interest to you, reach out to discuss more with our team.
Co-Folding Non-Canonical Amino Acids and Nucleotides
It’s now possible to run co-folding with non-canonical residues and nucleotides. These can be added directly into an input protein sequence and are available with all 3 co-folding models in Rowan (Boltz, OpenFold, and Chai).
To run co-folding with a non-canonical residue/nucleotide, simply include the desired PDB Chemical Component Dictionary (CCD) code in parentheses in the protein sequence: for instance, the sequence “AGA(SEP)G” requests the alanine–glycine–alanine–(O-phosphoserine)–glycine pentapeptide.
Use of CCD codes allows co-folding methods to model D amino acids, methylated residues, different terminal caps, phosphorylation, and many more biochemical processes of interest to scientists and drug hunters. For more on the CCD and how to find CCD codes for your non-canonical residue/nucleotide, see this article.
MM/GBSA Single Points
We’ve added the ability to refine poses with molecular mechanics and generalized-Born-surface-area implicit-solvent calculations (MM/GBSA for short). For docking, analogue docking, and protein–ligand co-folding, Rowan now optimizes all ligand poses using MM/GBSA implicit-solvent calculations and reports the final single-point MM/GBSA interaction energy. We’re doing this for a few reasons:
Poses generated by docking and co-folding are often unreliable or unphysical, and a quick MM/GBSA optimization can lead to dramatically more physical poses. (Previous Rowan pose refinement was ligand-only; this approach takes protein–ligand interactions into account.)
The MM/GBSA single-point energy also serves as an external sanity check on poses. While docking and co-folding will often find a way to put a pose somewhere when requested, poor poses often have negligible or even positive interaction energy.
Running these MM/GBSA optimizations also makes the systems much better behaved for downstream simulations, including MD and SAPT interaction analysis.
(Note that MM/GBSA interaction energies are not an effective way to predict binding affinity! However, they can be useful in ranking poses or triaging bad binders.)
MM/GBSA optimizations are run when users select “Local opt. and MM/GBSA” on the docking, analogue docking, and protein-ligand co-folding workflows.
Molecular Dynamics (MD) Speed Improvements
As more customers run MD workloads through Rowan, we’ve been working to make sure that our MD is as performant as possible. Increasing the timestep used for MD is one of the easiest ways to increase performance. To prevent numerical error, the timestep is typically required to be about an order-of-magnitude less than the fastest motion in the system. An O–H bond has a vibrational period of about 10 fs, implying that the MD timestep must be 1 fs or less.
This timestep can be increased through a technique called hydrogen-mass repartitioning (HMR), where the hydrogen mass is artificially increased to slow the bond vibration and allow for longer MD timesteps. Rowan previously used a hydrogen mass of 1.5 Daltons, allowing for 2 fs timestep. We’ve now changed the default hydrogen mass to 3 Daltons and the default timestep to 4 fs, resulting in a 2x increase in MD speed. Users can now tune these settings themselves if they prefer to disable HMR and use shorter timesteps.
In real life, this translates to dramatically faster MD runtimes. On a dasatinib•ABL kinase system (PDB 2GQG) with 55,517 atoms, Rowan was able to solvate, parameterize, minimize, equilibrate, and run a 10-ns MD simulation in 18 minutes total, reaching speeds of >1 µs/day on our standard L40S GPU hardware. We aim to make it routine to run short MD simulations on structures from docking, co-folding, or experiment, and will continue to look for ways to make these simulations faster and cheaper.
More MD Analysis
We’re also adding more outputs to Rowan’s protein MD and pose-analysis MD workflows, with the goal of making it easier to quickly gain useful structural insights. Here’s what we’ve added:
MM/GBSA. Rowan now reports MM/GBSA interaction energies for protein–ligand trajectories, making it easier to compare poses while accounting for solvent effects. Average MM/GBSA scores can be used as an endpoint-based binding-affinity prediction method; these scores are almost always quantitatively incorrect, but can sometimes be very useful for quickly ranking binders before running FEP calculations. (See our post comparing binding-affinity methods for more details.)
Radius of gyration. Rowan now reports the isotropic radius of gyration along trajectories, giving a simple readout of protein compactness and large-scale conformational change.
SASA and polar SASA. Rowan also now computes total and polar solvent-accessible surface area (SASA), helping track changes in solvent exposure and buried area during MD simulations.
Trajectory clustering. Geometrical clustering of trajectory frames is a convenient way to analyze protein conformational motion. Users can specify a fixed number of clusters to use k-means clustering or provide an RMSD threshold for Butina-like frame clustering—once computed, cluster centroids can be downloaded, and the change on cluster can be viewed as a function of time to understand transition timescales.
RMSF. Per-residue root-mean-squared fluctuation values are now calculated for proteins in trajectory post-processing, which can be used to highlight flexible loops and stable regions. Highly mobile regions (like loops or disordered regions) typically have high RMSF values, while regions that are rigid or relatively frozen have low RMSF values.
Hydration-site analysis. We’re also releasing a preview of Rowan’s hydration-site analysis feature, which detects what regions of the binding interface are most occupied by waters, how many unique waters pass through these sites, and whether waters at these mediate hydrogen bonds between protein atoms and ligand atoms (water bridges). This is a big feature and we’re planning to release more on this—including benchmarks and case studies—but subscribing Rowan users can run this starting today.
These properties can be viewed through new tabs and graphs in Rowan’s existing MD trajectory workflow.
Synthetic RBFE Intermediates
While RBFE calculations are a powerful way to accurately predict relative binding affinities, they require high overlap between a series. Modifications that look simple to a chemist (switching alkyl groups, adding or removing substituents, and so on) can sometimes lead to very difficult alchemical transformations, which in turn makes FEP slow and inaccurate.
Rowan FEP can now automatically identify “missing” compounds from an RBFE map based on maximum-common-substructure analysis and substituent enumeration, and automatically add these compounds to the graph:
This technique follows recent work from Linde Schoenmaker and co-workers, and can be enabled by selecting “Generate intermediate ligands” on the workflow submission screen.
While these sorts of synthetic intermediates are rarely useful on carefully prepared FEP benchmark sets, where there’s already good overlap between the ligands under study, we expect that this will be helpful in real-world use cases with more diverse chemotypes.
FEP Resubmission
In real-world lead optimization, scientists often go through rapid iterative design cycles. A medicinal chemist might design 20 analogues, simulate them with FEP, and then want to simulate another 20 analogues based on the results from the first set. Unfortunately, as discussed above, RBFE simulations work better on large sets of closely related compounds—so getting accurate results for the second set of analogues often requires re-running the first set of analogues, doubling the cost and runtime of the simulation.
Rowan FEP now allows users to resubmit finished FEP runs and add or remove ligands without re-running previously completed legs. This allows a few different useful features:
Iteratively adding ligands. If you’ve already run a big RBFE graph and want to add a few more ligands, you previously had to re-run the entire graph. Now, you can simply resubmit, add the new ligands, and run only the handful of legs needed to connect these ligands to the rest of the graph.
Removing bad ligands or poses. If a few of the ligands proved problematic or there were poorly converged legs, you can resubmit and remove these ligands to quickly recalculate the resulting well-behaved graph.
This feature was requested by a lot of our early FEP users, and we think it enables faster and more natural integration of FEP into real-world design cycles.
To enable this feature, simply resubmit an existing RBFE run to Rowan’s RBFE graph workflow. The previous FEP run will appear as a “seed graph” which serves as a template for the new run: the graph can be modified or deleted as desired. (The original FEP run will be saved and unchanged, allowing for clean record-keeping and audit trails.)
We’re excited to explore further ways to integrate FEP more naturally into existing drug-discovery workflows and make it possible to run FEP routinely at scale—if you’re interested in exploring this with our team, please reach out.
 | A guest post by
|
 |
|