Mutation Free Energy Calculations

What this workflow does

This workflow computes a fast-growth non-equilibrium mutation free energy. It uses GROMACS and pmx. The system is Staphylococcal nuclease (PDB: 1STN). The mutation is Ile10 to Ala, and the reverse.

The workflow starts from pre-equilibrated wild-type (stateA) and mutant (stateB) trajectories. It extracts snapshots from each one. It then models the alchemical hybrid structure for each state with pmxmutate, and builds a hybrid dual topology with pmxgentop.

The mutation introduces dummy atoms. The workflow freezes them with an index group, minimizes them, and runs an NPT equilibration for each state. It then runs fast thermodynamic integration in both directions: WT to Mut and Mut to WT.

The last stage, pmxanalyse, reads the accumulated dHdl work values. It estimates the free-energy difference with the Crooks Gaussian Intersection, BAR, and Jarzynski estimators.

The compute problem

Free-energy calculation is the most compute-heavy workflow in this collection, and for a specific reason: the estimators need statistics.

One TI transition is a short GROMACS run. It costs minutes. But CGI, BAR, and Jarzynski need dozens of independent transitions per direction before the estimate converges. A real ΔΔG study runs 80 or more transitions. That is 80 short GROMACS jobs that share nothing and can all run at once.

The rest of the graph is cheap. The pmxmutate, pdb2gmx, pmxgentop, and make_ndx stages take seconds. The final pmxanalyse stage reads text files and writes a number and a plot.

So the shape is one wide, embarrassingly parallel band of GPU work between two thin CPU layers. A single workstation runs that band in series and wastes days.

How Horus solves it

Horus assigns an executor per stage. Put the structure preparation and the analysis on your local CPU. Put the minimization, the equilibration, and the TI transitions on a GPU node or on an HPC scheduler. The scheduler runs the replicates concurrently.

The runtime.command string does not change when you move a stage. You change one executor: field. The same workflow.yaml runs on a laptop for a smoke test and on a cluster for the production sample.

The biobb_pmx package has no osx-arm64 conda package. On Apple Silicon, the pmx stages run in a quay.io/biocontainers/biobb_pmx linux/amd64 image through the Docker executor. The GROMACS stages keep running natively in the conda environment, because biobb_gromacs has an arm64 build. Horus mixes the two in one workflow.

Horus also caches finished stages. The equilibration output does not rerun when you only change the estimator options in configs/pmxanalyse.yaml.

Pipeline

gmx_trjconv_str_ens_stateA    Extract snapshots from WT (stateA) equilibrium trajectory
gmx_trjconv_str_ens_stateB    Extract snapshots from Mutant (stateB) equilibrium trajectory
   │
pmxmutate_stateA               Model stateA mutated structure (Ile10 to Ala, WT to Mut)
pmxmutate_stateB               Model stateB mutated structure (Ala10 to Ile, Mut to WT)
   │
pdb2gmx_stateA / pdb2gmx_stateB      Build GROMACS topology per state (pdb2gmx)
   │
pmxgentop_stateA / pmxgentop_stateB   Generate hybrid dual-topology parameters (pmxgentop)
   │
make_ndx_stateB                 Create FREEZE index for stateB dummy atoms (make_ndx)
   │
grompp_min_stateB ──► mdrun_min_stateB ──► gmx_energy_min_stateB   Minimize stateB dummy atoms
   │
grompp_eq_stateA ──► mdrun_eq_stateA ──► gmx_energy_eq_stateA      NPT equilibration, stateA
grompp_eq_stateB ──► mdrun_eq_stateB ──► gmx_energy_eq_stateB      NPT equilibration, stateB
   │
grompp_ti_stateA ──► mdrun_ti_stateA    Fast thermodynamic integration, WT to Mut
grompp_ti_stateB ──► mdrun_ti_stateB    Fast thermodynamic integration, Mut to WT
   │
pmxanalyse                      Compute ΔΔG (CGI / BAR / Jarzynski) from dhdl work values

Inputs and outputs

Inputs ship in pmx_tutorial/:

• stateA_1ns.xtc, stateA.tpr: the wild-type equilibrium trajectory.
• stateB_1ns.xtc, stateB.tpr: the mutant equilibrium trajectory.
• dhdlA.zip, dhdlB.zip: bundles of pre-computed replicate dHdl curves.

The pmxanalyse stage reads the two zip bundles. The estimators need dozens of independent transitions to converge. The single tiA_dhdl.xvg and tiB_dhdl.xvg files that this workflow produces show the method. They are not the statistical sample.

Outputs land in results/:

• mutA.pdb, mutB.pdb: the modelled mutant structures.
• pmxA_top.zip, pmxB_top.zip: the hybrid dual-topology files.
• eqoutA.gro, eqoutB.gro: the equilibrated structures.
• tiA.gro, tiA_dhdl.xvg, tiB.gro, tiB_dhdl.xvg: the TI output per direction.
• pmx.txt: the free-energy estimate and the statistics.
• pmx.plots.png: the work-distribution and convergence plots.

Run the workflow

Install the horus-runtime and the plugins 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, start Docker before you run the workflow. Horus pulls the biobb_pmx image on the first run.

References

• BioExcel Building Blocks (biobb)
• biobb_pmx
• pmx
• Official pmx Sardinia 2018 tutorial
• BioExcel mutation free energy tutorial