Self-avoiding walk

Usage examples and full code reference for SAW. For the underlying theory, see Self-avoiding walk.

Quick start

from afrc.polymer_models.saw import SAW

model = SAW('MASNDYTQQATQSYGAYPTQPGQGYSQQSSQPYG')

model.get_mean_end_to_end_distance()
model.get_root_mean_squared_end_to_end_distance()
model.get_mean_radius_of_gyration()

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

The dimensions are tuned with the prefactor argument (default 5.5 A), which can be passed to any of the methods:

for pref in (4.5, 5.5, 6.5):
    print(pref, model.get_mean_end_to_end_distance(prefactor=pref))

See also the demo/demo_SAW.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 = SAW('MASNDYTQQATQSYGAYPTQPGQGYSQQSSQPYG')

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

A Gaussian approximation: the mean-squared distances are exact, but the distance distributions are Gaussian and there is no excluded volume (see Self-avoiding walk).

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

Code reference

class afrc.polymer_models.saw.SAW(seq: str, p_of_r_resolution: float = 0.05)

Self-avoiding walk, using the des Cloizeaux scaling form as implemented by O’Brien et al. [1]. This model was developed by Jhullian ‘J’ Alston.

The scaling exponent is fixed at the good-solvent value, held on the object as self.nu = 0.598. That single value sets both the chain-length dependence of the size scale (\(R_{ee} = \texttt{prefactor}\,N^{\nu}\)) and the universal \(R_g/R_e\) ratio, so the two always describe the same chain. To vary the exponent, use NuDepSAW.

This is a composition-independent reference model: the sequence is only used to set the number of residues. The overall size is set by the prefactor argument accepted by every method (default 5.5 A).

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] Le Guillou, J. C., & Zinn-Justin, J. (1977). Critical Exponents for the n-Vector Model in Three Dimensions from Field Theory. Physical Review Letters, 39(2), 95-98.

__init__(seq: str, p_of_r_resolution: float = 0.05) → None

Create a SAW 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.

check_ensemble(conformations: NDArray[float64], prefactor: float = 5.5) → EnsembleReport

Check how well an ensemble reproduces the self-avoiding walk.

The exact checks are the mean-squared distance of every residue pair (in bands of separation) and the root-mean-square radius of gyration and first-to-last distance that follow from them. The SAW’s analytical (des Cloizeaux) end-to-end distribution and its radius of gyration from the universal SAW ratio are reported as context, together with a note that ensembles from sample_conformations() are a Gaussian approximation.

Parameters:

conformations (np.ndarray) – Array of shape [n_conformations x N x 3] (in Angstroms).
prefactor (float) – Size scale in Angstroms, so that the root-mean-square distance between residues k apart is prefactor * k**0.598. Default is 5.5 A.

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.
SAWException – If a parameter is out of range.

get_end_to_end_distribution(prefactor=5.5)

Return the end-to-end distance distribution.

Because prefactor can change between calls, the distribution is recomputed every time rather than cached.

Parameters:: prefactor (float) – Size scale in Angstroms, such that the root-mean-square end-to-end distance is prefactor * N**0.598. Default is 5.5 A; this should be tuned to match excluded-volume simulations for quantitative work.
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:: SAWException – If prefactor is not positive.

get_mean_end_to_end_distance(prefactor=5.5)

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)\).

Parameters:: prefactor (float) – Size scale in Angstroms (see get_end_to_end_distribution). Default is 5.5 A.
Returns:: The mean end-to-end distance (in Angstroms).
Return type:: float
Raises:: SAWException – If prefactor is not positive.

get_mean_radius_of_gyration(prefactor=5.5)

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

This is obtained from the root-mean-square end-to-end distance via the universal ratio

\[\frac{\langle R_g^2 \rangle}{\langle R_e^2 \rangle} = \frac{\gamma(\gamma + 1)}{2(\gamma + 2\nu)(\gamma + 2\nu + 1)}\]

using the object’s gamma and nu attributes, so the exponent here is by construction the one that sets the size scale. 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.

Parameters:: prefactor (float) – Size scale in Angstroms (see get_end_to_end_distribution). Default is 5.5 A.
Returns:: The root-mean-square radius of gyration (in Angstroms).
Return type:: float
Raises:: SAWException – If prefactor is not positive.

get_mean_squared_distance_map(prefactor: float = 5.5) → NDArray[float64]

Return the model’s mean-squared distance between every pair of residues.

Applying the model’s size scaling to every pair of residues, \(\langle r_{ij}^2 \rangle = (\texttt{prefactor}\, |i - j|^{\nu})^2\).

Parameters:: prefactor (float) – Size scale in Angstroms, so that the root-mean-square distance between residues k apart is prefactor * k**0.598. Default is 5.5 A.
Returns:: Symmetric [N x N] matrix of mean-squared distances (in Angstroms squared), with zeros on the diagonal.
Return type:: np.ndarray
Raises:: SAWException – If a parameter is out of range.

get_root_mean_squared_end_to_end_distance(prefactor=5.5)

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, and equals prefactor * N**0.598.

Parameters:: prefactor (float) – Size scale in Angstroms (see get_end_to_end_distribution). Default is 5.5 A.
Returns:: The root-mean-square end-to-end distance (in Angstroms).
Return type:: float
Raises:: SAWException – If prefactor is not positive.

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

Generate 3D conformations (one bead per residue) approximating the self-avoiding walk.

This is a Gaussian approximation: the bead coordinates are drawn exactly from the Gaussian chain whose mean-squared inter-residue distances are the model’s, \((\texttt{prefactor}\, |i - j|^{\nu})^2\), so the size and scaling of every inter-residue distance are right. The distance distributions are Gaussian rather than the SAW’s des Cloizeaux form, and there is no excluded volume, so beads can overlap. (Genuine self-avoiding conformations would need an explicit simulation.)

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.
prefactor (float) – Size scale in Angstroms, so that the root-mean-square distance between residues k apart is prefactor * k**0.598. Default is 5.5 A.

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.
SAWException – If a parameter is out of range.

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

Generate a self-avoiding walk ensemble (Gaussian approximation) 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.
prefactor (float) – Size scale in Angstroms, so that the root-mean-square distance between residues k apart is prefactor * k**0.598. Default is 5.5 A.

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.
SAWException – If a parameter is out of range.