Worm-like chain (O’Brien)

Usage examples and full code reference for WormLikeChain2. For the underlying theory, see Worm-like chain (O’Brien).

Quick start

from afrc.polymer_models.wlc2 import WormLikeChain2

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

model.get_mean_end_to_end_distance()
model.get_root_mean_squared_end_to_end_distance()
model.get_mean_radius_of_gyration()       # closed-form root-mean-square Rg

r, p = model.get_end_to_end_distribution()

Vary the persistence length:

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

Note

WormLikeChain2 requires the chain’s contour length (len(seq) * aa_size) to be at least one persistence length, otherwise a WLC2Exception is raised.

See also the demo/demo_WormLikeChain2.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 = WormLikeChain2('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 wlc2 (see Command-line tools).

Code reference

class afrc.polymer_models.wlc2.WormLikeChain2(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 O’Brien et al. (2009).

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\). Unlike the Zhou model (WormLikeChain) the O’Brien expression enforces finite extensibility exactly and stays well behaved for long chains, and this model also provides a closed-form radius of gyration.

References

[1] O’Brien, E. P., Morrison, G., Brooks, B. R., & Thirumalai, D. (2009). How accurate are polymer models in the analysis of Forster resonance energy transfer experiments on proteins? The Journal of Chemical Physics, 130(12), 124903.

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

Create a WormLikeChain2 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:

WLC2Exception – If lp or aa_size is not positive, or if the contour length (len(seq) * aa_size) is shorter than the persistence length.

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 (O’Brien 2009) and the Benoit-Doty radius of gyration are 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). No probability lies at or beyond the contour length.
Return type:: tuple of np.ndarray

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_radius_of_gyration()

Return the root-mean-square radius of gyration, \(\sqrt{\langle R_g^2 \rangle}\).

O’Brien et al. give \(\langle R_g^2 \rangle\) in closed form (they do not say so explicitly, but the expression is a mean-square). With \(C_2 = 1/(2L_p)\),

\[\langle R_g^2 \rangle = \frac{L_c}{6 C_2} - \frac{1}{4 C_2^2} + \frac{1}{4 C_2^3 L_c} - \frac{1 - e^{-L_c/L_p}}{8 C_2^4 L_c^2}\]

which is the Benoit-Doty worm-like chain result \(L_c L_p/3 - L_p^2 + 2L_p^3/L_c - 2L_p^4/L_c^2 (1 - e^{-L_c/L_p})\). It reduces to \(L_c L_p/3\) for a flexible chain and to the rigid-rod value \(L_c^2/12\) when \(L_c \ll L_p\).

Note that despite the method name this is the root-mean-square radius of gyration, not \(\langle R_g \rangle\). The name is kept for consistency with the other models.

Returns:: The root-mean-square radius of gyration (in Angstroms).
Return type:: float

References

[1] O’Brien, E. P., Morrison, G., Brooks, B. R., & Thirumalai, D. (2009). How accurate are polymer models in the analysis of Forster resonance energy transfer experiments on proteins? The Journal of Chemical Physics, 130(12), 124903.

[2] Benoit, H., & Doty, P. (1953). Light scattering from non-Gaussian chains. The Journal of Physical Chemistry, 57(9), 958-963.

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 (O’Brien 2009) 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.