Worm-like chain (Zhou)

Usage examples and full code reference for WormLikeChain. For the underlying theory, see Worm-like chain (Zhou).

Quick start

from afrc.polymer_models.wlc import WormLikeChain

# defaults: persistence length lp = 3.0 A, segment size aa_size = 3.8 A
model = WormLikeChain('MASNDYTQQATQSYGAYPTQPGQGYSQQSSQPYG')

model.get_mean_end_to_end_distance()
model.get_root_mean_squared_end_to_end_distance()

# full end-to-end distribution (distances, probabilities)
r, p = model.get_end_to_end_distribution()

Explore the effect of chain stiffness by varying the persistence length:

for lp in (2.0, 3.0, 4.0):
    wlc = WormLikeChain('MASNDYTQQATQSYG', lp=lp)
    print(lp, wlc.get_mean_end_to_end_distance())

See also the demo/demo_WormLikeChain.ipynb notebook for a worked, plotted example.

Generating 3D ensembles

Generate explicit conformations (one bead per residue), write them to disk, and check how well they reproduce the model:

model = WormLikeChain('MASNDYTQQATQSYGAYPTQPGQGYSQQSSQPYG', lp=3.0)

xyz = model.sample_conformations(n=5000, seed=1)                 # [n x N x 3] array, Angstroms
model.save_ensemble('ens', n=5000, seed=1)                       # ens.pdb + ens.xtc (needs mdtraj)
print(model.check_ensemble(xyz).format())                  # the model-check report

Or from the command line: afrc-ensemble -s <sequence> -m wlc (see Command-line tools).

Code reference

class afrc.polymer_models.wlc.WormLikeChain(seq: str, p_of_r_resolution: float = 0.05, lp: float = 3.0, aa_size: float = 3.8)

Worm-like chain model, as implemented by Zhou (2004).

This is a composition-independent reference model: the sequence is only used to set the number of residues, and hence the contour length \(L_c = N b\). It should agree closely with the O’Brien model (WormLikeChain2), but unlike that model it does not provide a radius of gyration.

The underlying expression is a series expansion in \(L_p/L_c\) and \(r/L_c\), so it is only accurate when the contour length comfortably exceeds the persistence length (for the default parameters, chains of more than ~10-20 residues). No probability is assigned beyond the contour length, and a chain too short for the expansion to have any valid region raises a WLCException when the distribution is requested.

References

[1] Zhou, H.-X. (2004). Polymer models of protein stability, folding, and interactions. Biochemistry, 43(8), 2141-2154.

__init__(seq: str, p_of_r_resolution: float = 0.05, lp: float = 3.0, aa_size: float = 3.8) → None

Create a WormLikeChain object.

Parameters:

seq (str) – Amino acid sequence. Only its length is used, and it is not validated.
p_of_r_resolution (float) – Grid spacing (in Angstroms) used for the distribution. Default is 0.05 A.
lp (float) – Persistence length, in Angstroms. We use a default of 3.0 A, although 4 A is also common in the literature. Must be > 0.
aa_size (float) – Contour length per residue (called \(b\) in the literature), in Angstroms. The default of 3.8 A is the Cα-Cα distance. Must be > 0.

Raises:

WLCException – If lp or aa_size is not positive.

check_ensemble(conformations: NDArray[float64]) → EnsembleReport

Check how well an ensemble reproduces the worm-like chain.

The exact checks are the mean-squared distance of every residue pair (in bands of separation), the root-mean-square radius of gyration and first-to-last distance that follow from them, and finite extensibility (no pair further apart than its contour length). The expected values are those of the discretized chain sample_conformations() draws from, which differs from the continuous worm-like chain by at most 0.02% (the report notes the exact figure). The model’s analytical end-to-end distribution (Zhou 2004) is reported as context.

Parameters:: conformations (np.ndarray) – Array of shape [n_conformations x N x 3] (in Angstroms).
Returns:: The checks; report.passed says whether the ensemble is model-like and report.format() gives a printable report.
Return type:: EnsembleReport
Raises:: AFRCException – If the conformations do not match this chain.

get_end_to_end_distribution()

Return the end-to-end distance distribution.

The distribution is computed on first use and then cached.

Returns:: (distances, probabilities), where distances are in Angstroms and the probabilities are a normalized probability mass function (they sum to 1).
Return type:: tuple of np.ndarray
Raises:: WLCException – If the contour length is too short relative to the persistence length for the Zhou expansion to be valid anywhere.

get_mean_end_to_end_distance()

Return the mean end-to-end distance, \(\langle R_e \rangle\).

This is the expectation over the end-to-end distribution, \(\sum r P(r)\).

Returns:: The mean end-to-end distance (in Angstroms).
Return type:: float

get_mean_squared_distance_map() → NDArray[float64]

Return the exact mean-squared distance between every pair of residues.

For beads k residues apart (contour length \(L = k b\)) this is the exact worm-like chain result \(\langle r^2 \rangle = 2 L_p L - 2 L_p^2 (1 - e^{-L/L_p})\).

Returns:: Symmetric [N x N] matrix of mean-squared distances (in Angstroms squared), with zeros on the diagonal.
Return type:: np.ndarray

get_root_mean_squared_end_to_end_distance()

Return the root-mean-square end-to-end distance, \(\sqrt{\langle R_e^2 \rangle}\).

This is the square root of \(\sum r^2 P(r)\) over the end-to-end distribution.

Returns:: The root-mean-square end-to-end distance (in Angstroms).
Return type:: float

sample_conformations(n: int = 1000, seed: int | Generator | None = None) → NDArray[float64]

Generate 3D conformations (one bead per residue) of the worm-like chain.

Each residue is split into short straight sub-segments, each bending away from the previous one with a fixed mean cosine, and only the bead at the end of each residue is kept. That correlation is either the worm-like chain’s own, \(e^{-s/L_p}\) for sub-segments of length \(s\), or - for chains that are flexible on the scale of a residue - one that reproduces its long-range size exactly with far fewer sub-segments (see afrc.ensemble.worm_like_chain_discretization()). Either way the number of sub-segments is chosen so the mean-squared distances match the continuous worm-like chain to within 0.02% at every separation. This samples the worm-like chain itself; the model’s analytical end-to-end distribution (Zhou 2004) is an approximation to it.

Parameters:

n (int) – Number of conformations. Default is 1000.
seed (int, np.random.Generator or None) – Seed (or generator) for reproducible ensembles. If None (default) a fresh, unpredictable seed is used.

Returns:

Array of shape [n x N x 3] with the bead coordinates, in Angstroms, each conformation centred on the origin.

Return type:

np.ndarray

Raises:

AFRCException – If n is not a positive integer.

save_ensemble(filename: str, n: int = 1000, seed: int | Generator | None = None, pdb_only: bool = False) → NDArray[float64]

Generate a worm-like chain ensemble and write it to disk.

Writes <filename>.pdb (topology and first conformation) and <filename>.xtc (every conformation; needs mdtraj), or with pdb_only=True a single multi-model <filename>.pdb. See sample_conformations() for how the conformations are generated.

Parameters:

filename (str) – Output path without an extension (a trailing .pdb or .xtc is dropped).
n (int) – Number of conformations. Default is 1000.
seed (int, np.random.Generator or None) – Seed (or generator) for reproducible ensembles.
pdb_only (bool) – Write a single multi-model PDB instead of a PDB/XTC pair. Default is False.

Returns:

The conformations that were written, shape [n x N x 3] (Angstroms).

Return type:

np.ndarray

Raises:

AFRCException – If n is invalid, the sequence has non-standard amino acids, mdtraj is missing (for a PDB/XTC pair), or the PDB format limits are exceeded.