Protein-Membrane MD Analysis

What this workflow does

This workflow analyzes a membrane molecular dynamics simulation. The example system is a heteropentameric GABA-gated chloride channel in a DPPC lipid bilayer. The trajectory comes from the MemProtMD project (MDDB accession: A023K).

The workflow fetches the structure, the trajectory, and the topology. It fits the trajectory twice. The first fit puts the membrane in the xy-plane. The second fit centers the protein with an XY translation only.

The workflow then identifies the two membrane leaflets with FATSLiM. Lipyphilic assigns each lipid to a leaflet per frame.

Seven analysis stages read that assignment. They compute the lipid z-positions and the bilayer thickness, both for the whole membrane and around the protein. They compute the deuterium order parameters of the acyl tails, the area per lipid by 2D Voronoi tessellation, the density profile along the membrane normal, the channel pore-radius profile, and the lipid flip-flop events between leaflets.

The compute problem

Membrane analysis is I/O-bound, not compute-bound. This changes what good hardware looks like.

The download stage pulls a full trajectory over the network. It is the slowest single stage on a home connection.

The two fitting stages stream every frame through GROMACS and write a new trajectory. Disk throughput sets their run time. A fast local disk beats a fast CPU here.

The seven analysis stages each read the fitted trajectory again. Seven readers over one large file is the pattern that punishes a network filesystem. mda_hole and gorder_aa are the heaviest of the seven.

Run these seven stages in sequence on a shared NFS mount and you pay the read cost seven times over the slowest possible link.

How Horus solves it

Horus places the analysis stages next to the data. The seven analysis stages are independent of each other. They all depend only on assign_leaflets. Horus dispatches them together on a node with local scratch disk. Each one reads the trajectory from local storage.

The two fitting stages are the ones to co-locate with the download. Keep the download, fit_membrane, and fit_protein on the same executor and the large intermediate trajectories never cross a network boundary. Only the small CSV results come back.

biobb_mem has no osx-arm64 conda package. On Apple Silicon the FATSLiM, Lipyphilic, and MDAnalysis stages run in the Docker executor against a Linux/amd64 biocontainers image. The biobb_io and biobb_analysis stages run natively in the conda executor. Each stage names its executor with one field.

The command strings never change between platforms. You develop the analysis on a Mac against a short trajectory. You then point the same workflow.yaml at a full production trajectory on an HPC node by editing the executor definitions.

Pipeline

download_trajectory        Download structure + trajectory for A023K from MDDB (mddb)
   │
fit_membrane                 Fit trajectory to membrane normal, rot+trans (gmx_image)
   │
fit_protein                   Fit trajectory to protein, XY translation only (gmx_image)
   │
identify_leaflets_fatslim      Identify membrane leaflets (fatslim_membranes)
   │
assign_leaflets                 Per-frame leaflet assignment (lpp_assign_leaflets)
   │
   ├── compute_zpositions           Lipid z-positions, whole membrane (lpp_zpositions)
   ├── compute_zpositions_around     Lipid z-positions around the protein (lpp_zpositions)
   ├── compute_order_parameters       Deuterium order parameters of acyl tails (gorder_aa)
   ├── compute_area_per_lipid          Area per lipid, 2D Voronoi tessellation (fatslim_apl)
   ├── compute_density_profile          Membrane density profile along z (cpptraj_density)
   ├── compute_pore_dimensions           Channel pore-radius profile (mda_hole)
   └── compute_flip_flop                 Lipid flip-flop events across leaflets (lpp_flip_flop)

Inputs and outputs

Inputs

The workflow has no file inputs. It fetches the structure, the trajectory, and the topology for MDDB entry A023K. configs/download_trajectory.yaml sets the node and project ID, the trajectory format, and the frame stride.

Outputs land in results/:

• A023K_10.pdb, A023K_10.xtc, A023K.tpr: the downloaded structure trajectory, and topology.
• A023K_fit.xtc, A023K_fit2.xtc: the membrane-fitted and protein-fitted trajectories.
• leaflets.ndx: the FATSLiM leaflet index groups.
• leaflets_data.csv: the per-frame Lipyphilic leaflet assignment. The flip-flop stage reads this file again.
• zpositions.csv, zpositions_around.csv: the lipid z-position and bilayer thickness data.
• order.csv: the deuterium order parameters.
• apl.csv: the area-per-lipid data.
• density.dat: the density profile along the membrane normal.
• hole.vmd, hole_profile.csv: the channel pore-radius profile.
• flip_flop.csv: the detected lipid flip-flop events.

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 biocontainers image on the first run.

References

• BioExcel Building Blocks (biobb)
• biobb_mem
• FATSLiM
• Lipyphilic
• MDDB
• BioExcel protein-membrane MD analysis tutorial