Metadata-Version: 2.4
Name: tunnelprofile
Version: 0.1.0
Summary: Tunnel detection and ligand transport energy profiles for proteins, without license-encumbered dependencies
Author-email: Miquel Estevez-Gay <miquel.estevezgay@unibe.ch>
License-Expression: MIT
Project-URL: Homepage, https://github.com/Buller-Lab/tunnelprofile
Project-URL: Documentation, https://github.com/Buller-Lab/tunnelprofile#documentation
Project-URL: Issues, https://github.com/Buller-Lab/tunnelprofile/issues
Keywords: protein tunnels,molecular docking,binding energy,caver,caverdock,openmm,autodock vina,computational enzymology
Classifier: Development Status :: 4 - Beta
Classifier: Intended Audience :: Science/Research
Classifier: Operating System :: POSIX :: Linux
Classifier: Operating System :: MacOS
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Topic :: Scientific/Engineering :: Bio-Informatics
Classifier: Topic :: Scientific/Engineering :: Chemistry
Requires-Python: >=3.10
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=1.24
Requires-Dist: scipy>=1.10
Provides-Extra: chem
Requires-Dist: rdkit>=2023.9; extra == "chem"
Requires-Dist: meeko>=0.6; extra == "chem"
Requires-Dist: vina>=1.2.5; extra == "chem"
Requires-Dist: gemmi>=0.6; extra == "chem"
Provides-Extra: plots
Requires-Dist: matplotlib>=3.7; extra == "plots"
Provides-Extra: dev
Requires-Dist: pytest>=7.4; extra == "dev"
Requires-Dist: ruff>=0.5; extra == "dev"
Requires-Dist: build>=1.0; extra == "dev"
Requires-Dist: twine>=5.0; extra == "dev"
Dynamic: license-file

# 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](#licensing).

```python
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:

| | what it is |
|---|---|
| LB (Vina) | a rigid conformer ensemble docked into every sphere |
| LB (force field) | those poses relaxed with OpenMM, best per sphere |
| UB (force field) | the same, constrained to one continuous trajectory |
| UB (Vina rescored) | the relaxed poses put back on the docking scale |

**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:

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

See [Dynamics](docs/dynamics.md) 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.

**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](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.

```bash
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:

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

<details>
<summary>If charges warn about NAGL on Fedora or RHEL</summary>

NAGL gives AM1-BCC-quality charges and needs torch. Some conda-forge torch
builds ship a `libtorch_cpu.so` that asks for an executable stack, which a
hardened kernel refuses; openff-toolkit then reports only "NAGL is not
available". The run detects this and prints the fix, which is one command:

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

Without it, charges fall back to AM1-BCC (slower, needs ambertools) and then
to Gasteiger (a real drop in quality). The warning says which was used.

</details>

<details>
<summary>What the environment contains, and why each piece is there</summary>

| package | for |
|---|---|
| `numpy`, `scipy` | the tunnel search itself: clearance field, graph, clustering |
| `openmm`, `pdbfixer` | protein repair and protonation, the restrained minimisation, MD |
| `rdkit` | ligand chemistry, conformers, bond-order perception |
| `openff-toolkit`, `openmmforcefields` | ligand and cofactor templates |
| `openff-nagl` | AM1-BCC-quality charges without AmberTools |
| `vina`, `meeko` | docking, and the receptor/ligand PDBQT preparation |
| `gemmi` | mmCIF, and the Chemical Component Dictionary |
| `matplotlib` | the static PNG profile from the batch runner (optional) |

</details>

### 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:

```bash
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:

```bash
pip install tunnelprofile
```

That install runs [example 01](examples/01_tunnels_only.py) 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](docs/dynamics.md):

```bash
mamba install -c conda-forge openmm-plumed
```

## Examples

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

| | | needs |
|---|---|---|
| [01_tunnels_only.py](examples/01_tunnels_only.py) | find and cluster tunnels, write the page and the PyMOL session | pip alone |
| [02_energy_profile.py](examples/02_energy_profile.py) | the full four-curve transport profile | full env |
| [03_trajectory.py](examples/03_trajectory.py) | an MD trajectory, clustered into routes with occupancies | full env |
| [04_reuse_tunnels.py](examples/04_reuse_tunnels.py) | search once, screen many ligands down the same route | full env |
| [05_cofactors_and_metals.py](examples/05_cofactors_and_metals.py) | a metalloenzyme, end to end, with nothing to configure | full env |
| [06_steered_md.py](examples/06_steered_md.py) | pull the ligand out and watch the protein respond | full env |
| [07_variant_screen.py](examples/07_variant_screen.py) | build point mutants, compare bottleneck against barrier | full env |
| [08_ligand_library.py](examples/08_ligand_library.py) | a ligand library across workers | full env |

Or from the command line:

```bash
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:

```bash
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

| | |
|---|---|
| [Examples](examples/README.md) | eight runnable scripts, and which to start from |
| [Inputs](docs/inputs.md) | structures, ligands, starting points, charges, tunnels as objects |
| [Cofactors, metals and disulfides](docs/cofactors.md) | Sage/NAGL parameters, the CCD, non-bonded metal centres |
| [How it works](docs/methods.md) | detection, the two stages, the force-field bounds, units |
| [Reading the results](docs/results.md) | the HTML page and the PyMOL sessions |
| [Dynamics](docs/dynamics.md) | steered MD, well-tempered metadynamics, the MD protocol |
| [Validation](docs/validation.md) | against CAVER, CaverDock, Vina and experiment |
| [Running it](docs/running.md) | timings, GPU, batch screening, the test suite |

## Why not CAVER / CaverDock / MGLTools

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

| role | replaced | used instead | license |
|---|---|---|---|
| tunnel detection | CAVER 3 | own implementation (`tunnel.py`) | this project |
| docking | — | AutoDock Vina 1.2 (Python API) | Apache-2.0 |
| PDBQT preparation | MGLTools / AutoDockTools | Meeko | LGPL-2.1 |
| protein repair, protonation | `pdb4amber` | PDBFixer + OpenMM | MIT |
| ligand chemistry, conformers | Open Babel | RDKit | BSD-3 |
| ligand and cofactor force field | GAFF via `antechamber` (GPL) | OpenFF Sage 2.2.1, typed from SMIRKS | MIT |
| geometry, graphs | — | numpy / scipy | BSD |

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](docs/methods.md)), 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:

| | Sage 2.2.1 | GAFF2 |
|---|---|---|
| drug-like organics, cofactors, porphyrins | yes | yes |
| molecular oxygen (`O=O`) | **no** | yes |

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:

```bash
mamba install -c conda-forge ambertools     # already present as a transitive dep
```
```python
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](docs/cofactors.md).
* 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:

```bash
./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/>:

| field | value |
|---|---|
| PyPI project | `tunnelprofile` |
| Owner | `Buller-Lab` |
| Repository | `tunnelprofile` |
| Workflow | `publish.yml` |
| Environment | `pypi` |

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

```bash
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](LICENSE). Read [Licensing](#licensing) above for what the
dependencies imply.
