tunnelprofile

Energy profiles for ligands moving through protein tunnels — tunnel detection, per-sphere docking, force-field refinement and steered MD, built entirely from permissively licensed components.

It is a replacement for the CAVER / CaverDock toolchain, which is free for academic use but not commercially. Everything here can be used commercially, including the parts under GPL — see Licensing.

from tunnelprofile import analyse

analyse("protein.pdb", ligand="BrCCBr", start=["A:108", "A:132", "A:272"],
        methods=["tunnels", "lb", "ub"], out="results")

One call gives you the routes, an energy profile along each, an interactive HTML page and PyMOL sessions, zipped up to carry to a workstation.

What it can do

Find tunnels. A clearance field, a cost-weighted shortest path, medial-axis centring and spline smoothing — no CAVER installation. Validated sphere by sphere against a CAVER reference on LinB, and against 10 reviewed tunnels across three structures.

Follow them through a trajectory. A multi-model PDB is a trajectory: frames are superposed on a common core, tunnels found in each, and the results clustered into routes with an occupancy and one representative each. Clustering is machinery, not a result — you are shown the route, not the twenty curves behind it.

Profile the cost of getting through. Four curves, all in kcal/mol:

Check the barrier dynamically. Steered MD (methods=["smd"]) pulled along the tunnel, with an adaptive force that re-aims at the next sphere as the ligand moves.

Well-tempered metadynamics is also implemented, on a path collective variable that turns with the tunnel so the substrate stays free to rotate — but it is a library-level API, not a methods= option, because a run is tens of minutes and its output is a free-energy surface rather than a profile the HTML page can show:

from tunnelprofile.metad import MetaDOptions, run_metadynamics
result = run_metadynamics(build, tunnel, start_pose_A, MetaDOptions())

See Dynamics for how to get build and start_pose_A, and for what the numbers do and do not mean.

Handle a real deposited structure. Cofactors, metals, halides and disulfides are kept and parameterised rather than deleted — see below.

Take a predicted structure as it comes. Boltz, AlphaFold3 and Chai name chains for what they hold (FE, AKG, SUB), leave the residues inside as LIG1/LIG2, and name atoms however they like. All of that is handled without an argument, and a site the predictor packed shut is reported with the width of the tightest point rather than a shrug — see Inputs.

Screen. Search once and reuse the routes for any number of ligands; build point mutants and compare bottleneck against barrier; run a ligand library across workers with the protein work done once.

Show you the result. One interactive HTML page per run — pick tunnel, frame, ligand and method, normalise the energy axis, swap spheres for ångströms — plus PyMOL sessions for the structures, the ligand poses and the steered-MD trajectories, and one zip containing all of it.

Cofactors, metals and disulfides

Deleting a cofactor opens the pocket it fills and drops the charges the substrate feels going past, so nothing is deleted for want of a template. Anything the protein force field cannot type is parameterised with OpenFF Sage and NAGL charges, its chemistry taken from the PDB Chemical Component Dictionary and matched to your structure by atom name. Ionisation is applied, not assumed — 2-oxoglutarate is deposited as the neutral diacid and used as the dianion that actually chelates the iron.

A metal inside a cofactor is split out and treated non-bonded: the organic part typed as though the metal were not there, the metal given the ion's own Lennard-Jones parameters, and the coordination shell held by distance restraints. On myoglobin that reproduces the six-coordinate iron — four pyrrole nitrogens, the proximal histidine, the bound dioxygen.

Where the CCD leaves a charge blank, which it does for every metal, you have to say what it is: cofactor_charge={"HEM": -2}. Reading a blank as neutral is how a ferric heme quietly becomes something else.

Disulfides are detected and bonded; free cysteines keep their thiol hydrogen.

Details, and the two limitations that remain, in docs/cofactors.md.

Install

The chemistry stack has no reliable wheels, so conda-forge is the supported way in. Either solver works; mamba is much faster.

git clone https://github.com/Buller-Lab/tunnelprofile.git
cd tunnelprofile

mamba env create -f environment.yml     # or: conda env create -f environment.yml
mamba activate tunnelprofile            # or: conda activate tunnelprofile

pip install -e .                        # or: pip install tunnelprofile

environment.yml pins nothing tighter than a lower bound and solves to about 270 packages. Check it worked:

python -c "import tunnelprofile; print(tunnelprofile.__version__)"
tunnelprofile --help
python examples/01_tunnels_only.py

pip install patchelf
patchelf --clear-execstack "$CONDA_PREFIX/lib/libtorch_cpu.so"

No GPU?

The energy profiles run on the CPU, roughly 30x slower — the LinB example takes 9 minutes instead of 17 seconds. CPU is never chosen automatically, because a silent drop to it is orders of magnitude slower with no visible cause. Ask for it explicitly, once:

export TUNNELPROFILE_PLATFORM=CPU

or per run: --platform CPU, or analyse(..., platform="CPU"). Tunnel detection does not use OpenMM and is unaffected.

Just the tunnels, without the chemistry

Tunnel detection, clustering, the HTML page and the PyMOL sessions need only numpy and scipy, and nothing chemical is imported until a method asks for it:

pip install tunnelprofile

That install runs example 01 as-is. The energy profiles need the conda environment above.

PLUMED, for OPES only

Not needed for anything above, and note that the conda-forge package does not currently ship the opes module — see Dynamics:

mamba install -c conda-forge openmm-plumed

Examples

Runnable, from the repository root. Full index and guidance in examples/README.md.

Or from the command line:

tunnelprofile analyse protein.pdb --ligand "BrCCBr" \
    --start "A:108,A:132,A:272" --methods tunnels,lb,ub --out results

analyse and batch are the whole command set. --methods chooses what runs: any of tunnels, lb, ub, smd.

Reusing a previous run

Every run writes tunnels.json and keeps the prepared frames beside it, so a different ligand, a different route or different settings cost only the profiling:

tunnelprofile analyse protein.pdb --start "A:108,A:132" --methods tunnels --out run1
tunnelprofile analyse --tunnels run1/tunnels.json --ligand "CCBr" \
    --methods lb,ub --tunnel 2 --out run2

The frames are copied into run1/frames/ and the saved ensemble points there. That is not tidiness: a tunnel is a set of coordinates, so it only means anything alongside the structure it was found in. Saving the ensemble while its frames sat in a scratch directory that then gets deleted produces a file that loads fine and fails the moment anything touches a frame -- so the load checks, and says which files are missing.

Outputs

One directory, and one zip of it to carry to a workstation:

results/
  profiles.html         every profile, interactive: pick tunnel, frame, ligand, method
  tunnels.json          the routes, reusable with tunnels= (see below)
  frames/               the prepared structures the tunnels belong to
  structures/           always: protein, tunnel spheres, bottlenecks, view.pml
  energies/             if lb or ub ran: the above plus the ligand pose per sphere
  smd/                  if smd ran: pulling trajectories, protein included
results.zip

Each PyMOL directory has a view.pml; pymol view.pml loads that session with the tunnels grouped and the representative shown. profiles.html is one page for the whole run, not one per folder.

Documentation

Why not CAVER / CaverDock / MGLTools

Those are free for academic use but not for commercial use. Everything here is replaced with a permissive equivalent:

LGPL-2.1 (Meeko) and MPL-2.0 (gemmi, pulled in by Meeko) are library licenses that do not impose obligations on your own code.

Licensing

This project is MIT, nothing in the stack restricts commercial use, and the defaults call no GPL code.

The defaults are GPL-free, and that is tested

Ligands and cofactors are typed with OpenFF Sage 2.2.1 from SMIRKS and charged with NAGL, both MIT, neither needing an external program. The protein is ff14SBonlysc (see How it works), which ships with openmmforcefields (MIT). The suite builds WelO5 — metal, halide, cosubstrate, the lot — with antechamber and sqm made unreachable, so this is checked rather than claimed.

ambertools is not a dependency of this package. It will still be installed in your conda environment, because openmmforcefields declares it, and there is nothing we can do about that short of dropping openmmforcefields. Installed is not invoked: the default path never runs it.

And using it would be fine anyway

GPL restricts distribution, not use. The GPL-3 text shipped with AmberTools:

This License explicitly affirms your unlimited permission to run the unmodified Program. The output from running a covered work is covered by this License only if the output, given its content, constitutes a covered work.

You may make, run and propagate covered works that you do not convey, without conditions.

So running antechamber commercially carries no obligation, and charges it computes are output, not a covered work. This is categorically different from CAVER, CaverDock and MGLTools, which restrict use itself to academic contexts — that restriction is what this project exists to escape, and GPL does not impose it. Within AmberTools, note also that the GAFF2 parameters are public domain (dat/leap/parm/gaff2.dat); only the programs are GPL.

The one case with an obligation is redistribution: put AmberTools inside a container or installer you hand to someone else and you are conveying GPL-3 code, and must pass on its source and licence — which is public, so in practice a pointer and a notice. It does not reach your own code.

When you may still want GAFF2

Sage covers organic chemistry. It does not cover everything, and the gap is easy to hit:

On myoglobin that costs the bound dioxygen: the heme is typed and the iron restrained to its four pyrrole nitrogens and the proximal histidine, but OXY is dropped and the sixth coordination position with it. The run names the component and the remedy rather than failing obscurely. To get it back:

mamba install -c conda-forge ambertools     # already present as a transitive dep

analyse(..., ligand_ff="gaff-2.11", cofactor_ff="gaff-2.11",
             charge_method="am1bcc", cofactor_charge_method="am1bcc")

That is a deliberate opt-in to GPL tooling, for the reasons above.

None of this is legal advice; it is what the licence files say and how these arrangements are generally understood. For a commercial deliverable, have your own legal people read it.

Known limitations

• GBn2 implicit solvent, not explicit water. Adequate for ranking; not a free energy. No water-mediated hydrogen bonds, and no real viscosity, so egress kinetics from the dynamics modules are not physical.
• The force-field barrier is pose-sampling limited: repeating the profile with different docking seeds moves the upper bound by +/- 6 kcal/mol, and raising the pose budget from seeds=16, keep_top=1 to seeds=32, keep_top=3 does not fix it. Treat a single force-field barrier as an estimate, not a measurement. This dominates every other uncertainty in stage 2, including the GB cutoff, which was checked and cleared.
• Freezing protein atoms does not speed up MD — OpenMM still evaluates their forces. --mobile-radius and --rigid-protein control accuracy and stability, not cost.
• The metadynamics barrier on LinB (25.35 +/- 3.89 kcal/mol) exceeds the +16.4 kcal/mol ceiling that kcat allows, and the PMF does not bind the substrate at all. The restrained pocket cannot open its 1.33 A bottleneck for a ligand needing 2.2-2.5 A. A larger mobile shell and softer backbone restraints are the untested next step; until then read these barriers as a property of the model.
• Metadynamics transitions are one-way here (0.6 +/- 0.5 recrossings per 10 ns walker), so the reweighted free energy has no statistics at the barrier and only the bias-inversion estimate is usable.
• SMD forward and reverse work distributions overlap poorly, so BAR carries a large uncertainty that is not yet computed.
• Only one system, one ligand and one tunnel have been run this way.
• OPES is implemented but unusable: the installed PLUMED has no opes module.
• Tunnel detection is validated quantitatively against one CAVER reference (LinB WT) plus a functional check on L177W. More reference tunnels would strengthen this.
• Multiple-tunnel enumeration is now checked against 10 reviewed tunnels on three LinB structures, but all three are the same fold; other folds are untested.
• Auxiliary tunnels narrower than the probe are invisible by construction — see the probe discussion above before trusting a tunnel count.
• In the whole-box search test the ensemble often finds a different pose at similar energy (RMSD 1.3–3.7 Å) until the seed budget is large. Inside a sphere this matters much less, but raise --seeds for unusually large or flexible ligands.
• NAGL charges fall back to AM1-BCC via AmberTools when torch is unavailable; --charge-method gasteiger is the fast, less accurate option.
• Stage 2 energies are interaction energies, not potentials of mean force. Use them to compare positions along one tunnel, not to compare ligands.
• Meeko cannot build a receptor template for some cofactors -- a heme is one -- so those structures get no docking. The run falls back to force-field profiles from geometric starting poses and says so; see Cofactors.
• A metal inside a cofactor is treated non-bonded: no charge transfer, no covalency, and the coordination geometry is held where it was deposited rather than derived.

Continuous integration

Every push runs lint and minimal-install (3.10, 3.11, 3.12) — about a minute, and the second one is the guard that matters most: it proves the package still works from pip alone, which one stray top-level import would break invisibly.

The chemistry suite runs weekly and on demand, not per push. GitHub's runners have no GPU, the energy tests are roughly 30x slower on the CPU platform, and even the subset marked "not slow" was still going after an hour. Run it by hand from Actions → tests → Run workflow when a change touches the chemistry, and locally with pytest -q before a release — which scripts/release.sh does for you.

Releasing

The version lives in one place, tunnelprofile/__init__.py; pyproject.toml reads it. One command does the rest:

./scripts/release.sh 0.2.0 --dry     # checks and builds, changes nothing
./scripts/release.sh 0.2.0           # bumps, tests, tags, pushes
gh release create v0.2.0 --generate-notes

Publishing the GitHub release triggers .github/workflows/publish.yml, which uploads to PyPI over Trusted Publishing — no API token stored in the repository or in Actions secrets. Configure that once at https://pypi.org/manage/account/publishing/:

To rehearse against TestPyPI first — worth doing, because a version number on PyPI can never be reused:

export TWINE_PASSWORD='pypi-...'     # a TestPyPI token
./scripts/release.sh 0.2.0 --test

The script refuses to run on a dirty tree, off main, with a version that is already tagged, or with one already on PyPI.

Citing and contributing

Issues and pull requests welcome at https://github.com/Buller-Lab/tunnelprofile. If you use this in published work, cite the repository and the tools it stands on: AutoDock Vina, OpenMM, RDKit, Meeko, PDBFixer, OpenFF and AmberTools.

License

MIT — see LICENSE. Read Licensing above for what the dependencies imply.