Classical Molecular Interaction Potentials with CMIP

What this workflow does

This workflow computes Classical Molecular Interaction Potentials from protein structures. It uses the CMIP code through the BioExcel Building Blocks. It runs three independent case studies in one graph.

Case 1 computes Molecular Interaction Potential grids on Lysozyme (PDB: 1AKI). It produces a positive grid, a negative grid, and a neutral grid. It also titrates structural waters and ions onto the surface.

Case 2 computes protein-ligand interaction energies for the EGFR kinase domain bound to Erlotinib (PDB: 4HJO). This case generates an AMBER topology with tleap, runs an energy minimization with sander, and parameterizes the ligand with ACPype and GAFF.

Case 3 computes protein-protein interaction energies between the SARS-CoV-2 receptor-binding domain and human ACE2. It reads a pre-run AMBER structure and topology.

The compute problem

The three cases have almost nothing in common at the hardware level.

Case 1 is a set of grid calculations. Each MIP grid is a separate CMIP call. The three grids are independent. They need CPU and memory proportional to the grid volume.

Case 2 is a preparation chain with a real MD stage inside it. The sander_mdrun minimization is the slowest single stage. AMBER tools want a Linux build. ACPype adds a chemistry toolchain on top.

Case 3 is a set of grid boxes and masked interaction calculations. It skips the setup because the AMBER files already exist.

Run all three cases in sequence on one machine and you wait for the AMBER stage before the Lysozyme grids finish. The cases do not depend on each other, so that wait buys you nothing.

How Horus solves it

Horus runs the three cases at the same time. There is no edge between them in the graph. Case 1 and Case 3 finish while Case 2 is still in tleap.

The executor routing matters here more than in most biobb workflows. Both biobb_cmip and biobb_amber lack an osx-arm64 conda package. On Apple Silicon, the cmip_run, cmip_titration, leap_gen_top, and sander_mdrun stages run in the Docker executor against a quay.io/biocontainers Linux/amd64 image. The biobb_io, biobb_structure_utils, and biobb_chemistry stages run natively in the conda executor.

Two executor definitions sit at the top of workflow.yaml. Each stage points at one of them. The runtime.command strings never change. The same file runs on a Linux workstation with a single conda executor and on a Mac with the mixed strategy.

To move sander_mdrun to an HPC node, change one field. The grid stages stay where they are. You do not restructure the workflow to relocate one expensive stage.

Pipeline

Case 1, Lysozyme MIP & titration
fetch_pdb_lysozyme        Download 1AKI from the PDB (pdb)
   │
cmip_prepare_pdb_lysozyme Add CMIP charges/elements (cmip_prepare_pdb)
   ├── cmip_titration_lysozyme ──► cat_pdb_titration   Structural waters/ions (cmip_titration / cat_pdb)
   ├── cmip_run_mip_pos       Positive MIP grid (cmip_run)
   ├── cmip_run_mip_neg       Negative MIP grid (cmip_run)
   └── cmip_run_mip_neutral   Neutral MIP grid (cmip_run)

Case 2, EGFR–Erlotinib protein-ligand interaction
fetch_pdb_complex ──► remove_pdb_water ──► extract_heteroatoms ──► reduce_add_hydrogens ──► acpype_params_ac
   (download 4HJO, strip water, extract AQ4 ligand, add H, GAFF params via ACPype)
   │
leap_gen_top ──► sander_mdrun ──► amber_to_pdb   AMBER topology, minimization, PDB conversion
   │
cmip_prepare_structure_complex ──► remove_ligand
   ├── cmip_ignore_residues_protein   Mask protein for interaction calc (cmip_ignore_residues)
   └── cmip_run_egfr_interaction      EGFR–Erlotinib interaction energies (cmip_run)

Case 3, RBD–hACE2 protein-protein interaction (from a pre-run MD structure)
cmip_prepare_structure_md (Files/RBD-hACE2-ZN.pdb + .top)
   ├── extract_chain_hace2 / extract_chain_rbd     Split complex by chain (extract_chain)
   ├── cmip_run_rbd_box / cmip_run_hace2_box / cmip_run_complex_box   Grid boxes (cmip_run)
   ├── cmip_ignore_residues_rbd ──► cmip_run_rbd_interaction     RBD vs hACE2 energies
   └── cmip_ignore_residues_hace2 ──► cmip_run_hace2_interaction hACE2 vs RBD energies

Inputs and outputs

Inputs

The workflow fetches 1AKI and 4HJO from the PDB. Case 3 reads two bundled files: Files/RBD-hACE2-ZN.pdb and Files/RBD-hACE2-ZN.top. The workflow does not generate those two files.

Outputs land in results/:

• 1aki.pdb, 1aki.cmip.pdb, 1aki.wat_ions.pdb, 1aki.tit.pdb: the Lysozyme preparation and titration structures.
• 1aki.mip_pos.cube, 1aki.mip_neg.cube, 1aki.mip_neutral.cube: the MIP grids in Gaussian cube format.
• 4hjo.pdb, AQ4.pdb, AQ4.reduce.H.pdb, AQ4params.*: the EGFR complex and the Erlotinib parameters.
• structure.leap.prmtop, structure.leap.crd, structure.amber-min.pdb: the AMBER topology and the minimized complex.
• 4hjo.EGFR.energies.log, 4hjo.EGFR.ener.byat.out: the EGFR-Erlotinib interaction energies.
• RBD-hACE2-ZN.cmip.pdb and the per-chain CMIP structures: the split complex.
• RBD.box.json, hACE2.box.json, COMPLEX.box.json: the grid box definitions.
• RBD.energies.log, hACE2.energies.log and the matching .byat.out files: the protein-protein interaction energies.

Run the workflow

Install the horus-runtime one time:

uv sync

If you do not have uv, install it first:

curl -LsSf https://astral.sh/uv/install.sh | sh

You can also install the packages with pip:

pip install horus-runtime horus-environments

Then run the workflow:

uv run horus run workflow.yaml

On Apple Silicon you also need Docker. Horus pulls the CMIP and AMBER images on the first run.

References

• BioExcel Building Blocks (biobb)
• biobb_cmip
• CMIP
• AMBER / AmberTools
• ACPype
• BioExcel CMIP tutorial