Analytical Flory Random Coil

Usage examples and full code reference for AnalyticalFRC. For the underlying theory, see Analytical Flory Random Coil.

Quick start

Build an object from a sequence and read off ensemble-average dimensions:

from afrc import AnalyticalFRC

P = AnalyticalFRC('MASNDYTQQATQSYGAYPTQPGQGYSQQSSQPYGQQSYSGYSQSTDTSGYG')

P.get_mean_end_to_end_distance()                 # mean Re (A)
P.get_mean_radius_of_gyration()                  # mean Rg (A)
P.get_mean_hydrodynamic_radius()                 # mean Rh (A), Kirkwood-Riseman
P.get_mean_hydrodynamic_radius('nygaard')        # mean Rh (A), Nygaard et al. (2017)

Pull out full probability distributions (each returns (distances, probabilities)):

re_r, re_p = P.get_end_to_end_distribution()
rg_r, rg_p = P.get_radius_of_gyration_distribution()

# distribution between two specific residues
d_r, d_p = P.get_interresidue_distance_distribution(4, 40)

Inter-residue quantities and whole-chain maps (residue indices start at 0):

P.get_mean_interresidue_distance(4, 40)  # mean distance between two residues (A)
P.get_contact_fraction(4, 40, 10.0)      # fraction of time they are within 10 A
P.get_internal_scaling()                 # [|i-j|, mean distance] profile
P.get_distance_map()                     # n x n mean inter-residue distances
P.get_contact_map(15.0)                  # contact fractions at a 15 A threshold
P.get_pre_profile(0)                     # expected PRE profile for a spin label at residue 0

Draw a size-matched sample (e.g. to compare against a simulation trajectory):

samples = P.sample_end_to_end_distribution(n=5000)

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

Generating 3D ensembles

The AFRC can also generate explicit 3D conformations, one bead per residue, drawn exactly from the model (see Analytical Flory Random Coil for how). Save them as a PDB/XTC pair for analysis or visualization, or work with the coordinates directly:

# writes p53.pdb (topology + first conformation) and p53.xtc (all 5000 conformations)
xyz = P.save_ensemble('p53', n=5000, seed=1)

# or just the coordinates: an [n x N x 3] array in Angstroms
xyz = P.sample_conformations(n=5000, seed=1)

Writing the XTC file needs mdtraj (pip install "afrc[ensemble]"). The pair loads directly into mdtraj or SOURSOP:

import mdtraj as md
traj = md.load('p53.xtc', top='p53.pdb')

from soursop.sstrajectory import SSTrajectory
protein = SSTrajectory('p53.xtc', 'p53.pdb').proteinTrajectoryList[0]

Note that adjacent beads are not at a fixed 3.8 Å spacing and can overlap: the ensemble follows the AFRC’s own inter-residue distance distributions, including for neighbouring residues. Pass pdb_only=True to save_ensemble() to write every conformation to a single multi-model PDB file instead (no mdtraj needed, at most 9999 conformations).

Checking an ensemble

check_ensemble() reports how well an ensemble reproduces the AFRC - every quantity the model fixes exactly, each with a standard error estimated from the ensemble:

report = P.check_ensemble(xyz)
report.passed          # True for a model-like ensemble
print(report.format())

The same thing is available from the command line, for this and every other model, with afrc-ensemble - see Command-line tools.

Code reference

class afrc.AnalyticalFRC(seq: str, adaptable_P_res: bool = False)

The Analytical Flory Random Coil for a single amino acid sequence.

This is the main object the package provides. You build it from one sequence and then ask it polymer questions - mean dimensions, full distributions, inter-residue distances, contact fractions, PRE profiles and so on.

from afrc import AnalyticalFRC

protein = AnalyticalFRC('KFGGPRDQGSRHDSEQDNSDNNTIFVQGLG')
protein.get_mean_radius_of_gyration()

seq

The sequence (upper-case).

Type:

str

p_of_r_resolution

Grid spacing (in Angstroms) used for every distribution.

Type:

float

worm_like_chain

A Zhou worm-like chain with the same number of residues (default lp and aa_size), provided for convenience.

Type:

WormLikeChain

Notes

Nothing is computed until it is requested, so building the object is cheap. Methods that visit every pair of residues (distance maps, contact maps, internal scaling and the Kirkwood-Riseman \(R_h\)) build an [n x n] matrix of sub-chains on first use, which is the most expensive step, and then keep it. Methods that need one pair, or one row of pairs (the PRE profile), build only the segments they need.

Inter-residue quantities between residues \(i\) and \(j\) treat the segment between them as a chain of \(|i - j|\) residues, taking its composition from the sequence between the two. Whole-chain quantities use all \(N\) residues, so get_mean_interresidue_distance(0, N-1) is slightly smaller than get_mean_end_to_end_distance().

__init__(seq: str, adaptable_P_res: bool = False) → None

Create an AnalyticalFRC object from an amino acid sequence.

Parameters:

• seq (str) – Amino acid sequence (case insensitive). Only the 20 standard amino acids are allowed.

• adaptable_P_res (bool) – If False (default) every distribution uses a fixed grid spacing of 0.05 A. If True the spacing is set to \(d_{max}/500\), where \(d_{max} = 3.7N\) approximates the contour length. This changes the discretization only, not the model.

Raises:

AFRCException – If seq is not a string, is empty, or contains a non-standard amino acid.

check_ensemble(conformations: NDArray[float64]) → EnsembleReport

Check how well an ensemble reproduces the AFRC’s statistics.

Every quantity the AFRC fixes exactly is measured in the ensemble, with a standard error estimated from the ensemble, and compared with the model: the root-mean-square radius of gyration, the first-to-last bead distance, the Kirkwood-Riseman hydrodynamic radius, the mean-squared and mean inter-residue distances (in bands of sequence separation), and the full distribution of three representative distances. The AFRC’s separately calibrated \(\langle R_g \rangle\) and its whole-chain \(\langle R_e \rangle\) are reported as context. This is the report that afrc-ensemble prints.

Parameters:

conformations (np.ndarray) – Array of shape [n_conformations x N x 3] (in Angstroms), e.g. from sample_conformations().

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 sequence.

get_contact_fraction(R1, R2, threshold)

Return the fraction of the time two residues are closer than a threshold.

This is the cumulative probability of the inter-residue distance distribution below threshold. With a threshold of around 5 A this gives the expected contact fraction for a pair of residues under the AFRC, which is a useful normalization factor.

Parameters:

• R1 (int) – Index of the first residue (zero-based).

• R2 (int) – Index of the second residue (zero-based).

• threshold (float) – Distance threshold (in Angstroms).

Returns:

The contact fraction, between 0 and 1. A residue is always in contact with itself, so R1 == R2 returns 1.0.

Return type:

float

Raises:

AFRCException – If either residue index is invalid or threshold is not a number.

get_contact_map(threshold, symmetric_map=False)

Return the contact fraction for every pair of residues.

Parameters:

• threshold (float) – Distance threshold (in Angstroms) - see get_contact_fraction().

• symmetric_map (bool) – If True, return the full symmetric matrix. If False (default), only the upper triangle (including the diagonal) is filled and the lower triangle is zero.

Returns:

An [n x n] matrix of contact fractions, each between 0 and 1. The diagonal is 1.

Return type:

np.ndarray

Raises:

AFRCException – If threshold is not a number.

get_distance_map(calculation_mode='scaling law', symmetric_map=False)

Return the mean inter-residue distance for every pair of residues.

Parameters:

• calculation_mode (str) – Either 'scaling law' (default), which uses \(\langle r_{ij} \rangle = R_0 |i - j|^{0.5}\), or 'distribution', which takes the expectation over each pair’s distance distribution. The two agree to within about 0.3%.

• symmetric_map (bool) – If True, return the full symmetric matrix. If False (default), only the upper triangle is filled and the lower triangle is zero.

Returns:

An [n x n] matrix of mean inter-residue distances (in Angstroms).

Return type:

np.ndarray

Raises:

AFRCException – If calculation_mode is not recognized.

get_end_to_end_distribution()

Return the end-to-end distance (\(R_e\)) distribution.

This is the Gaussian chain distribution (see [Rubinstein2003])

\(P(r) = 4\pi r^2 \left( \frac{3}{2\pi \langle r^2 \rangle} \right)^{3/2} e^{-\frac{3 r^2}{2 \langle r^2 \rangle}}\)

with \(\sqrt{\langle r^2 \rangle} = R_0^{rms} N^{0.5}\).

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

get_internal_scaling(calculation_mode='scaling law')

Return the internal scaling profile.

This is the mean distance between residues that are \(|i - j|\) apart in sequence, averaged over every such pair. A log-log fit of the profile gives a slope of 0.5.

Parameters:

calculation_mode (str) – Either 'scaling law' (default) or 'distribution' - see get_distance_map().

Returns:

An [n-1 x 2] matrix with one row per sequence separation. The first column is \(|i - j|\) (1 to n-1) and the second is the mean distance (in Angstroms) at that separation.

Return type:

np.ndarray

Raises:

AFRCException – If calculation_mode is not recognized.

get_interresidue_distance_distribution(R1, R2)

Return the distance distribution between two residues.

Parameters:

• R1 (int) – Index of the first residue (zero-based).

• R2 (int) – Index of the second residue (zero-based). The order of R1 and R2 does not matter.

Returns:

(distances, probabilities), where distances are in Angstroms and the probabilities are a normalized probability mass function (they sum to 1). If R1 == R2 all the weight sits at zero.

Return type:

tuple of np.ndarray

Raises:

AFRCException – If either residue index is invalid.

get_mean_end_to_end_distance(calculation_mode: str = 'scaling law') → float

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

Parameters:

calculation_mode (str) – Either 'scaling law' (default), which uses \(\langle R_e \rangle = R_0 N^{0.5}\), or 'distribution', which takes the expectation over the end-to-end distribution. The two agree to within about 0.3%.

Returns:

The mean end-to-end distance (in Angstroms).

Return type:

float

Raises:

AFRCException – If calculation_mode is not recognized.

get_mean_hydrodynamic_radius(calculation_mode='kirkwood-riseman')

Return the mean hydrodynamic radius, \(R_h\).

In 'kirkwood-riseman' mode (default) this is

\[R_h = \left\langle \frac{1}{r_{ij}} \right\rangle_{i \neq j}^{-1}\]

where the average runs over every pair of residues and, for each pair, over its Gaussian distance distribution. The inner average has the closed form \(\langle 1/r_{ij} \rangle = \sqrt{6 / (\pi \langle r_{ij}^2 \rangle)}\), so the result is exact for the model. Note that this is the mean of the inverse distance, which is what the Kirkwood-Riseman equation needs - not the inverse of the mean distance, which for a Gaussian chain is larger by a factor of \(4/\pi\). This is the same form used by Nygaard et al. [1], Pesce et al. [3] and SOURSOP’s mode='kr'.

In 'nygaard' mode we instead apply the empirical \(R_g\)-to- \(R_h\) conversion of Nygaard et al. [1] to the mean radius of gyration.

Parameters:

calculation_mode (str) – Either 'kirkwood-riseman' (default) or 'nygaard'.

Returns:

The mean hydrodynamic radius (in Angstroms).

Return type:

float

Raises:

AFRCException – If calculation_mode is not recognized, or the chain has fewer than two residues (Kirkwood-Riseman then has no pairs to average over, and the Nygaard denominator \(N^{0.60} - N^{0.33}\) is zero).

References

[1] Nygaard, M., Kragelund, B. B., Papaleo, E., & Lindorff-Larsen, K. (2017). An efficient method for estimating the hydrodynamic radius of disordered protein conformations. Biophysical Journal, 113(3), 550-557.

[2] Kirkwood, J. G., & Riseman, J. (1948). The intrinsic viscosities and diffusion constants of flexible macromolecules in solution. The Journal of Chemical Physics, 16(6), 565-573.

[3] Pesce, F., Newcombe, E. A., Seiffert, P., Tranchant, E. E., Olsen, J. G., Grace, C. R., Kragelund, B. B., & Lindorff-Larsen, K. (2023). Assessment of models for calculating the hydrodynamic radius of intrinsically disordered proteins. Biophysical Journal, 122(2), 310-321.

get_mean_interresidue_distance(R1, R2, calculation_mode='scaling law')

Return the mean distance between two residues.

Parameters:

• R1 (int) – Index of the first residue (zero-based).

• R2 (int) – Index of the second residue (zero-based). The order of R1 and R2 does not matter.

• calculation_mode (str) – Either 'scaling law' (default), which uses \(\langle r_{ij} \rangle = R_0 |i - j|^{0.5}\), or 'distribution', which takes the expectation over the pair’s distance distribution.

Returns:

The mean distance between R1 and R2 (in Angstroms), or 0.0 if R1 == R2.

Return type:

float

Raises:

AFRCException – If either residue index or calculation_mode is invalid.

get_mean_interresidue_radius_of_gyration(R1, R2, calculation_mode='scaling law')

Return the mean radius of gyration of the segment between two residues.

Parameters:

• R1 (int) – Index of the first residue (zero-based).

• R2 (int) – Index of the second residue (zero-based). The order of R1 and R2 does not matter.

• calculation_mode (str) – Either 'scaling law' (default), which uses \(\langle R_g \rangle = R_0^{g} |i - j|^{0.5}\), or 'distribution', which takes the expectation over the segment’s \(R_g\) distribution.

Returns:

The mean radius of gyration of the segment (in Angstroms), or 0.0 if R1 == R2.

Return type:

float

Raises:

AFRCException – If either residue index or calculation_mode is invalid.

get_mean_radius_of_gyration(calculation_mode: str = 'distribution') → float

Return the mean radius of gyration, \(\langle R_g \rangle\).

Parameters:

calculation_mode (str) – Either 'distribution' (default), which takes the expectation over the \(R_g\) distribution, or 'scaling law', which uses \(\langle R_g \rangle = R_0^{g} N^{0.5}\) with the composition-weighted prefactor. The two agree to well within 0.1%.

Returns:

The mean radius of gyration (in Angstroms).

Return type:

float

Raises:

AFRCException – If calculation_mode is not recognized.

get_mean_squared_distance_map() → NDArray[float64]

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

For residues \(i\) and \(j\) this is \(\langle r_{ij}^2 \rangle = (R_0^{rms})^2 |i - j|\), with \(R_0^{rms}\) the average per-residue prefactor over the segment between them - exactly the value used by the inter-residue distance distributions. Because every AFRC inter-residue distance is Gaussian, this matrix fully specifies the model’s pairwise statistics (and is what sample_conformations() draws from). Unlike get_distance_map() it does not build the inter-residue matrix; every pair is computed at once from a running sum of the per-residue prefactors, so it is fast even for long sequences.

Returns:

Symmetric [n x n] matrix of mean-squared distances (in Angstroms squared), with zeros on the diagonal.

Return type:

np.ndarray

get_pre_profile(label_position, tau_c=4, t_delay=12, R_2D=14, W_H=3769911184.3077517, sample_size=10000)

Return the expected paramagnetic relaxation enhancement (PRE) profile for a spin label at label_position.

For every residue we draw sample_size distances to the label from the AFRC, convert each to a transverse relaxation rate

\[\Gamma_2 = \frac{K}{r^6} \left( 4\tau_c + \frac{3\tau_c}{1 + \omega_H^2 \tau_c^2} \right)\]

with \(K = 1.23 \times 10^{-32}\) cm⁶ s^-2, and then average the intensity ratio \(R_{2D} e^{-\Gamma_2 t} / (R_{2D} + \Gamma_2)\) over those samples. Averaging per conformer (rather than converting a mean distance) matters, because relaxation depends very non-linearly on distance.

Note that the model does not account for the spin-label linker. It gives the profile expected if the chain behaved as an AFRC chain.

Parameters:

• label_position (int) – Index of the labelled residue (zero-based).

• tau_c (float) – Effective correlation time, in nanoseconds. Typically between 1 and 30. Default is 4.

• t_delay (float) – Total duration of the INEPT delays in the PRE experiment, in milliseconds. This depends on the pulse sequence, but is typically 1-30 ms for an HSQC. Default is 12.

• R_2D (float) – Transverse relaxation rate of the backbone amide protons in the diamagnetic protein, in Hz (per second). Around 10 is typical. Default is 14.

• W_H (float) – Proton Larmor frequency as an angular frequency, in rad/s (\(\omega_H = 2\pi\nu_H\)). For a 600 MHz magnet pass 2*np.pi*600e6 (about 3.77e9), not 600000000; passing the linear frequency makes the dispersive term \((2\pi)^2 \approx 39.5\) times too small and over-estimates \(\Gamma_2\) by roughly 10% at the default tau_c. This matches the convention used by SOURSOP’s SSPRE class. Default is 2*np.pi*600e6 (a 600 MHz magnet).

• sample_size (int) – Number of distances drawn per residue. Larger values give a smoother profile at the cost of more compute. Default is 10000.

Returns:

A 3-element list:

[0] - residue indices (np.ndarray, starting at 0)

[1] - the PRE intensity-ratio profile (one value between 0 and 1 per residue). The labelled residue itself is 0.

[2] - the per-conformer relaxation rates \(\Gamma_2\) (one np.ndarray of length sample_size per residue, in s^-1). These are infinite for the labelled residue.

Return type:

list

Raises:

AFRCException – If label_position is not a valid residue index.

References

[1] Meng, W., Lyle, N., Luan, B., Raleigh, D. P., & Pappu, R. V. (2013). Experiments and simulations show how long-range contacts can form in expanded unfolded proteins with negligible secondary structure. Proceedings of the National Academy of Sciences, 110(6), 2123-2128.

[2] Das, R. K., Huang, Y., Phillips, A. H., Kriwacki, R. W., & Pappu, R. V. (2016). Cryptic sequence features within the disordered protein p27Kip1 regulate cell cycle signaling. Proceedings of the National Academy of Sciences, 113(20), 5616-5621.

[3] Peran, I., Holehouse, A. S., Carrico, I. S., Pappu, R. V., Bilsel, O., & Raleigh, D. P. (2019). Unfolded states under folding conditions accommodate sequence-specific conformational preferences with random coil-like dimensions. Proceedings of the National Academy of Sciences, 116(25), 12301-12310.

[4] Lalmansingh, J. M., Keeley, A. T., Ruff, K. M., Pappu, R. V., & Holehouse, A. S. (2023). SOURSOP: A Python package for the analysis of simulations of intrinsically disordered proteins. Journal of Chemical Theory and Computation, 19(16), 5609-5620.

get_radius_of_gyration_distribution()

Return the radius of gyration (\(R_g\)) distribution.

This uses equation 3 of [Lhuillier1988] with the composition-weighted AFRC prefactor.

Returns:

(radii, probabilities), where radii are in Angstroms and the probabilities are a normalized probability mass function (they sum to 1).

Return type:

tuple of np.ndarray

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

Generate 3D conformations (one bead per residue) drawn from the AFRC.

The AFRC defines every inter-residue distance as Gaussian, with a mean-squared value set by the composition of the sequence between the two residues. Taken together, those distances define a Gaussian chain whose bead coordinates are jointly Gaussian, with covariance

\[C = -\frac{1}{6} J D J, \qquad D_{ij} = \langle r_{ij}^2 \rangle, \qquad J = I - \frac{1}{n}\mathbf{1}\mathbf{1}^T\]

for each of x, y and z. We factor \(C\) once per sequence and then draw conformations from it directly. Every inter-residue distance in the resulting ensemble has exactly the AFRC’s distance distribution, for every pair of residues at once - this is not an approximation or a fit.

There are a few things worth knowing about the ensembles:

• Adjacent beads are not at a fixed 3.8 A spacing. Their separation follows the AFRC’s own distribution for neighbouring residues (a mean of ~5.8 A, with most values between 2 and 10 A).

• The AFRC is an ideal chain, so there is no excluded volume and beads can overlap.

• The end-to-end distance of the ensemble is the distance between the first and last residues, which the AFRC treats as a chain of N-1 residues. It therefore matches get_mean_interresidue_distance(0, N-1), which is slightly smaller than get_mean_end_to_end_distance() (N residues).

• The radius of gyration is that of the Gaussian chain, which agrees with get_mean_radius_of_gyration() to within ~2% (the AFRC’s Rg was calibrated separately from its distances).

Each conformation is centred on the origin and randomly oriented.

Parameters:

• n (int) – Number of conformations to generate. Default is 1000.

• seed (int, np.random.Generator or None) – Seed (or generator) for the random numbers, 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, where N is the sequence length.

Return type:

np.ndarray

Raises:

AFRCException – If n is not a positive integer.

Examples

>>> protein = AnalyticalFRC('MEEPQSDPSVEPPLSQETFSDLWKLLPENNVLSPLPSQAMDDLMLSPDDI')
>>> xyz = protein.sample_conformations(n=5000, seed=1)
>>> xyz.shape
(5000, 50, 3)

sample_end_to_end_distribution(n=1000)

Draw random end-to-end distances from the end-to-end distribution.

Useful for building a size-matched, uncorrelated sample to compare against simulation data.

Parameters:

n (int) – Number of values to draw. Default is 1000.

Returns:

n independent end-to-end distances (in Angstroms).

Return type:

np.ndarray

sample_inter_residue_distance_distribution(R1, R2, n=1000)

Draw random distances between two residues from their distance distribution.

Useful for building a size-matched, uncorrelated sample to compare against simulation data.

Parameters:

• R1 (int) – Index of the first residue (zero-based).

• R2 (int) – Index of the second residue (zero-based). The order of R1 and R2 does not matter.

• n (int) – Number of values to draw. Default is 1000.

Returns:

n independent distances (in Angstroms). If R1 == R2 every value is 0.

Return type:

np.ndarray

Raises:

AFRCException – If either residue index is invalid.

sample_radius_of_gyration_distribution(n=1000)

Draw random radii of gyration from the \(R_g\) distribution.

Useful for building a size-matched, uncorrelated sample to compare against simulation data.

Parameters:

n (int) – Number of values to draw. Default is 1000.

Returns:

n independent radii of gyration (in Angstroms).

Return type:

np.ndarray

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

Generate an AFRC ensemble and write it to disk.

This draws n conformations with sample_conformations() and, by default, writes them as a PDB/XTC pair: <filename>.pdb (the topology, one CA bead per residue, holding the first conformation) and <filename>.xtc (every conformation). Load them together with, for example, mdtraj.load('ens.xtc', top='ens.pdb') or SOURSOP’s SSTrajectory('ens.xtc', 'ens.pdb'). Writing the XTC file needs mdtraj (pip install mdtraj).

With pdb_only=True every conformation is instead written to a single multi-model <filename>.pdb, which needs no extra dependencies but is larger and limited to 9999 conformations.

Parameters:

• filename (str) – Output path without an extension; .pdb (and .xtc) are added. If the name already ends in .pdb or .xtc that extension is dropped first.

• n (int) – Number of conformations to generate. Default is 1000.

• seed (int, np.random.Generator or None) – Seed (or generator) for the random numbers, for reproducible ensembles. If None (default) a fresh, unpredictable seed is used.

• pdb_only (bool) – If True, write a single multi-model PDB file instead of a PDB/XTC pair. Default is False.

Returns:

The conformations that were written, as an array of shape [n x N x 3] in Angstroms.

Return type:

np.ndarray

Raises:

AFRCException – If n is not a positive integer, mdtraj is not installed (for a PDB/XTC pair), the sequence is longer than the PDB format allows (9999 residues), there are more than 9999 conformations (with pdb_only=True), or a conformation that goes in the PDB file does not fit in the PDB coordinate format.

Ensemble utilities

The functions behind sample_conformations(), save_ensemble() and afrc-ensemble, which can also be used to write, or check, your own one-bead-per-residue coordinates.

ensemble.py

Tools for generating and writing 3D conformational ensembles with one bead per residue.

There are two kinds of generator here:

• gaussian_chain_factor() / sample_gaussian_chain() work for any Gaussian chain - one in which every inter-residue distance is Gaussian-distributed. Such a chain is fully specified by its matrix of mean-squared inter-residue distances, and conformations are drawn from it exactly: the bead coordinates are jointly Gaussian, with a covariance fixed by that matrix. The AFRC uses this, as do the self-avoiding walk models (as a Gaussian approximation with the right mean-squared distances).

• sample_freely_jointed_chain(), sample_freely_rotating_chain() and sample_worm_like_chain() build chains bond by bond from the model’s own rules, so they are exact (the worm-like chain up to a controlled discretization). Their exact mean-squared inter-residue distances are given by freely_rotating_chain_msd() and discrete_worm_like_chain_msd().

Ensembles are written as a PDB/XTC pair: the PDB holds the topology (one CA bead per residue) and the first conformation, and the XTC holds every conformation. Writing the XTC needs mdtraj, which is an optional dependency. Alternatively, write_multimodel_pdb() writes the whole ensemble as a single multi-model PDB file, with no extra dependencies.

Copyright Alex Holehouse 2018-2026 (holehouselab.com).

afrc.ensemble.discrete_worm_like_chain_msd(k: NDArray[float64], segment_length: float, subdivisions: int, tangent_correlation: float) → NDArray[float64]

Exact mean-squared distance between beads k residues apart on the discretized worm-like chain generated by sample_worm_like_chain().

Each residue is split into subdivisions straight sub-segments of length \(s\), and successive sub-segments have tangent correlation \(c\), so this is the freely rotating chain result with \(k m\) bonds of length \(s\) and \(\alpha = c\).

Parameters:

• k (np.ndarray) – Separation in residues.

• segment_length (float) – Contour length per residue (in Angstroms).

• subdivisions (int) – Sub-segments per residue, \(m\).

• tangent_correlation (float) – Mean cosine between successive sub-segments, \(c\).

Returns:

Mean-squared distances (in Angstroms squared).

Return type:

np.ndarray

afrc.ensemble.freely_rotating_chain_msd(k: NDArray[float64], bond_length: float, cos_angle: float) → NDArray[float64]

Exact mean-squared distance between beads k bonds apart on a freely rotating chain.

\[\langle r^2 \rangle = k b^2 \frac{1 + \alpha}{1 - \alpha} - 2 b^2 \frac{\alpha (1 - \alpha^k)}{(1 - \alpha)^2}\]

where \(\alpha\) is the cosine of the angle between successive bonds. With \(\alpha = 0\) this is \(k b^2\), which is also exact for the freely jointed chain.

Parameters:

• k (np.ndarray) – Number of bonds between the two beads.

• bond_length (float) – Bond length \(b\) (in Angstroms).

• cos_angle (float) – \(\alpha\), strictly between -1 and 1.

Returns:

Mean-squared distances (in Angstroms squared).

Return type:

np.ndarray

afrc.ensemble.gaussian_chain_factor(mean_squared_distances: NDArray[float64]) → NDArray[float64]

Factor the coordinate covariance of a Gaussian chain.

For a Gaussian chain whose conformations are centred on the origin, each Cartesian coordinate of the beads is jointly Gaussian with covariance

\[C = -\frac{1}{6} J D J, \qquad J = I - \frac{1}{n}\mathbf{1}\mathbf{1}^T\]

where \(D_{ij} = \langle r_{ij}^2 \rangle\) is the matrix of mean-squared inter-residue distances. This returns a matrix \(A\) with \(C = A A^T\), so that A @ z for a standard normal vector z is one coordinate (x, y or z) of one conformation.

Parameters:

mean_squared_distances (np.ndarray) – Symmetric [n x n] matrix of mean-squared inter-residue distances (in Angstroms squared), with zeros on the diagonal.

Returns:

The [n x n] factor \(A\).

Return type:

np.ndarray

Raises:

AFRCException – If the matrix is not square and symmetric, or if the distances cannot be realized by any Gaussian chain in three dimensions (the covariance has a genuinely negative eigenvalue).

afrc.ensemble.mean_squared_distance_map(n_res: int, msd_of_separation: Callable[[NDArray[float64]], NDArray[float64]]) → NDArray[float64]

Build an [n x n] mean-squared distance matrix from a function of separation.

For the homogeneous chain models (everything except the AFRC) the mean-squared distance between two residues depends only on how far apart they are in sequence, \(k = |i - j|\).

Parameters:

• n_res (int) – Number of residues.

• msd_of_separation (callable) – Function mapping an array of separations \(k \geq 1\) to the mean-squared distances (in Angstroms squared).

Returns:

Symmetric [n_res x n_res] matrix with zeros on the diagonal.

Return type:

np.ndarray

afrc.ensemble.sample_freely_jointed_chain(n_res: int, bond_length: float, n_conformations: int, rng: Generator) → NDArray[float64]

Draw freely jointed chains: bonds of fixed length pointing in independent, uniformly random directions.

Parameters:

• n_res (int) – Number of beads.

• bond_length (float) – Bond length (in Angstroms).

• n_conformations (int) – Number of conformations.

• rng (np.random.Generator) – Random number generator.

Returns:

Bead coordinates, shape [n_conformations x n_res x 3], each conformation centred on the origin.

Return type:

np.ndarray

afrc.ensemble.sample_freely_rotating_chain(n_res: int, bond_length: float, cos_angle: float, n_conformations: int, rng: Generator) → NDArray[float64]

Draw freely rotating chains: bonds of fixed length with a fixed angle between successive bonds and uniformly random torsions.

Parameters:

• n_res (int) – Number of beads.

• bond_length (float) – Bond length (in Angstroms).

• cos_angle (float) – Cosine of the angle between successive bond vectors.

• n_conformations (int) – Number of conformations.

• rng (np.random.Generator) – Random number generator.

Returns:

Bead coordinates, shape [n_conformations x n_res x 3], each conformation centred on the origin.

Return type:

np.ndarray

afrc.ensemble.sample_gaussian_chain(factor: NDArray[float64], n_conformations: int, rng: Generator) → NDArray[float64]

Draw conformations of a Gaussian chain from its covariance factor.

Parameters:

• factor (np.ndarray) – The [n x n] covariance factor from gaussian_chain_factor().

• n_conformations (int) – Number of conformations to draw.

• rng (np.random.Generator) – Random number generator to draw from.

Returns:

Array of shape [n_conformations x n x 3] with the bead coordinates. Each conformation is centred on the origin and randomly oriented.

Return type:

np.ndarray

afrc.ensemble.sample_worm_like_chain(n_res: int, segment_length: float, n_conformations: int, rng: Generator, subdivisions: int, tangent_correlation: float) → NDArray[float64]

Draw worm-like chains.

Each residue is split into subdivisions straight sub-segments of length \(s\). The direction of each sub-segment is drawn from a von Mises-Fisher distribution about the previous one, with the concentration set so the mean cosine between successive sub-segments is tangent_correlation, and only the bead at the end of each residue is kept. Use worm_like_chain_discretization() to choose the discretization; the exact mean-squared distances of the resulting chain are discrete_worm_like_chain_msd().

Parameters:

• n_res (int) – Number of beads.

• segment_length (float) – Contour length per residue (in Angstroms).

• n_conformations (int) – Number of conformations.

• rng (np.random.Generator) – Random number generator.

• subdivisions (int) – Sub-segments per residue.

• tangent_correlation (float) – Mean cosine between successive sub-segments, strictly between 0 and 1.

Returns:

Bead coordinates, shape [n_conformations x n_res x 3], each conformation centred on the origin.

Return type:

np.ndarray

Raises:

AFRCException – If tangent_correlation is not strictly between 0 and 1.

afrc.ensemble.save_conformations(conformations: NDArray[float64], sequence: str, filename: str, pdb_only: bool = False, remark: str = '') → None

Write conformations to disk as a PDB/XTC pair or a single multi-model PDB.

This is the shared back end of every model’s save_ensemble().

Parameters:

• conformations (np.ndarray) – Array of shape [n_conformations x n x 3] (in Angstroms).

• sequence (str) – Amino acid sequence of length n (case insensitive).

• filename (str) – Output path without an extension; .pdb (and .xtc) are added. A trailing .pdb or .xtc is dropped first.

• pdb_only (bool) – If True, write one multi-model PDB instead of a PDB/XTC pair.

• remark (str) – Text for the REMARK lines of the PDB file.

Raises:

AFRCException – If the inputs are inconsistent or the files cannot be written (see write_ensemble() and write_multimodel_pdb()).

afrc.ensemble.worm_like_chain_discretization(segment_length: float, persistence_length: float, n_res: int, tolerance: float = 0.0002, max_subdivisions: int = 4096) → tuple[int, float]

Choose how to discretize the worm-like chain.

Each residue is split into straight sub-segments, and we need the smallest number for which the discretized chain’s mean-squared distances match the continuous worm-like chain to within tolerance (relative) at every separation in the chain. Two choices of the correlation between successive sub-segments are tried:

• the worm-like chain’s own, \(e^{-s/L_p}\), which is best for stiff chains;

• \((2L_p - s)/(2L_p + s)\), which reproduces the continuous chain’s long-range size exactly and leaves only a short-range error of \(s^2/2\), so needs far fewer sub-segments when \(L_p\) is small compared with a residue (e.g. 221 rather than 807 per residue at \(L_p = 0.1\) A).

Whichever needs fewer sub-segments is used.

Parameters:

• segment_length (float) – Contour length per residue (in Angstroms).

• persistence_length (float) – Persistence length (in Angstroms).

• n_res (int) – Number of residues.

• tolerance (float) – Largest acceptable relative difference from the continuous chain. Default is 2e-4.

• max_subdivisions (int) – Upper limit on the number of sub-segments. Default is 4096.

Returns:

Sub-segments per residue, and the correlation between successive sub-segments.

Return type:

tuple of (int, float)

afrc.ensemble.worm_like_chain_msd(contour_length: NDArray[float64], persistence_length: float) → NDArray[float64]

Exact mean-squared end-to-end distance of a continuous worm-like chain.

\[\langle R^2 \rangle = 2 L_p L - 2 L_p^2 \left(1 - e^{-L/L_p}\right)\]

Parameters:

• contour_length (np.ndarray) – Contour length(s) \(L\) (in Angstroms).

• persistence_length (float) – Persistence length \(L_p\) (in Angstroms).

Returns:

Mean-squared distances (in Angstroms squared).

Return type:

np.ndarray

afrc.ensemble.write_ensemble(conformations: NDArray[float64], sequence: str, pdb_filename: str, xtc_filename: str, remark: str = '') → None

Write an ensemble as a PDB/XTC pair.

The PDB file holds the topology and the first conformation; the XTC file holds every conformation (including the first). Load them together with, for example, mdtraj.load(xtc_filename, top=pdb_filename) or SOURSOP’s SSTrajectory(xtc_filename, pdb_filename).

Parameters:

• conformations (np.ndarray) – Array of shape [n_conformations x n x 3] with the bead coordinates, in Angstroms.

• sequence (str) – One-letter amino acid sequence of length n (upper-case).

• pdb_filename (str) – Path of the PDB file to write.

• xtc_filename (str) – Path of the XTC file to write.

• remark (str) – Optional text written as a REMARK line at the top of the PDB file.

Raises:

AFRCException – If the conformations and sequence are inconsistent, a coordinate in the first conformation is too large for the PDB format, or mdtraj is not installed.

afrc.ensemble.write_multimodel_pdb(conformations: NDArray[float64], sequence: str, filename: str, remark: str = '') → None

Write a set of one-bead-per-residue conformations as a multi-model PDB file.

Each conformation is written as its own MODEL/ENDMDL block, laid out as in write_pdb(), so the whole ensemble lives in one file that mdtraj (and most viewers) read as a trajectory. This needs no extra dependencies, but the file is several times larger than a PDB/XTC pair, and the PDB format limits it to 9999 conformations.

Parameters:

• conformations (np.ndarray) – Array of shape [n_conformations x n x 3] with the bead coordinates, in Angstroms.

• sequence (str) – One-letter amino acid sequence of length n (upper-case).

• filename (str) – Path of the PDB file to write.

• remark (str) – Optional text written as REMARK lines at the top of the file.

Raises:

AFRCException – If the conformations and sequence are inconsistent, there are more than 9999 conformations, or a coordinate is too large to fit in the fixed-width PDB format.

afrc.ensemble.write_pdb(coordinates: NDArray[float64], sequence: str, filename: str, remark: str = '') → None

Write a single one-bead-per-residue conformation as a PDB file.

Each residue is written as a single CA atom, with the residue name taken from the sequence, in chain A. Consecutive beads are joined with CONECT records so that viewers draw the chain even though the beads are not at standard CA-CA spacing.

Parameters:

• coordinates (np.ndarray) – Array of shape [n x 3] with the bead coordinates, in Angstroms.

• sequence (str) – One-letter amino acid sequence of length n (upper-case).

• filename (str) – Path of the PDB file to write.

• remark (str) – Optional text written as a REMARK line at the top of the file.

Raises:

AFRCException – If the coordinates and sequence are inconsistent, or a coordinate is too large to fit in the fixed-width PDB format.

afrc.ensemble.write_xtc(conformations: NDArray[float64], filename: str) → None

Write a set of conformations as an XTC trajectory.

XTC files store coordinates in nanometres at a precision of 0.001 nm (0.01 A), so coordinates are converted from Angstroms on the way out.

Parameters:

• conformations (np.ndarray) – Array of shape [n_conformations x n x 3] with the bead coordinates, in Angstroms.

• filename (str) – Path of the XTC file to write.

Raises:

AFRCException – If mdtraj (needed to write XTC files) is not installed, or the array has the wrong shape.

ensemble_report.py

Checks that a generated 3D ensemble reproduces the statistics of the chain model it was drawn from.

Every model we generate ensembles for has exactly known mean-squared inter-residue distances \(D_{ij} = \langle r_{ij}^2 \rangle\), and these fix the root-mean-square radius of gyration and end-to-end distance too. For a Gaussian chain (the AFRC) every distance also follows a Maxwell distribution set by \(D_{ij}\), so the mean distances, the Kirkwood-Riseman hydrodynamic radius and the full distance distributions can be checked as well. Models with rigid geometry (fixed bond lengths or angles, or a maximum extension) add exact structural checks.

compare_ensemble_to_model() measures each of these in the ensemble, with a standard error estimated from the ensemble itself, and compares it with the model’s exact expectation, described by a ModelExpectations. A correct ensemble agrees with every check to within sampling error.

Copyright Alex Holehouse 2018-2026 (holehouselab.com).

class afrc.ensemble_report.Check(name: str, passed: bool | None, summary: str, values: dict[str, float]=<factory>, structural: bool = False)

One comparison between the ensemble and the model.

name

What is being compared.

Type:

str

passed

True if the ensemble agrees with the model, False if it does not, and None for a context-only comparison (or one that could not be assessed).

Type:

bool or None

summary

Human-readable result (may span several lines).

Type:

str

values

The numbers behind the summary, for programmatic use.

Type:

dict of str to float

structural

True for an exact geometric check (e.g. a fixed bond length), which does not depend on sampling and so is assessed for any number of conformations.

Type:

bool

class afrc.ensemble_report.EnsembleReport(model_name: str, n_residues: int, n_conformations: int, n_pair_conformations: int, checks: list[~afrc.ensemble_report.Check], context: list[~afrc.ensemble_report.Check], notes: list[str] = <factory>)

The result of comparing an ensemble with its model.

model_name

Name of the model the ensemble was drawn from.

Type:

str

n_residues

Number of residues (beads).

Type:

int

n_conformations

Number of conformations in the ensemble.

Type:

int

n_pair_conformations

Number of conformations used for the pair-based checks (all of them, unless that would exceed MAX_PAIR_EVALUATIONS distances).

Type:

int

checks

Comparisons the ensemble must pass to be consistent with the model.

Type:

list of Check

context

Comparisons that are informative but not expected to match exactly.

Type:

list of Check

notes

Model-specific caveats.

Type:

list of str

property assessed: bool

Whether the statistical checks could be assessed (enough conformations and residues).

property failed: list[str]

Names of the checks that failed.

format(header: list[str] | None = None) → str

Format the report as plain text.

Parameters:

header (list of str, optional) – Extra lines (e.g. the sequence and output files) shown under the title.

Returns:

The report, ready to print or save.

Return type:

str

property passed: bool

Whether the statistical checks were assessed and every check passed.

class afrc.ensemble_report.ModelExpectations(model_name: str, mean_squared_distances: NDArray[float64], gaussian: bool = False, bond_length: float | None = None, bond_angle_cosine: float | None = None, contour_length_per_residue: float | None = None, reference_rg: float | None = None, reference_rg_label: str = '', reference_re: float | None = None, reference_re_label: str = '', end_to_end_distribution: tuple[NDArray[float64], NDArray[float64]] | None = None, end_to_end_label: str = '', notes: tuple[str, ...] = ())

What a chain model says an ensemble drawn from it should look like.

model_name

Name of the model, for display.

Type:

str

mean_squared_distances

The model’s exact [N x N] matrix of mean-squared inter-residue distances (Angstroms squared).

Type:

np.ndarray

gaussian

True if every inter-residue distance is exactly Gaussian (Maxwell distributed), which enables the mean-distance, hydrodynamic-radius and distribution checks. Default is False.

Type:

bool

bond_length

Exact length of every bond between consecutive beads, if the model fixes it (Angstroms).

Type:

float, optional

bond_angle_cosine

Exact cosine of the angle between consecutive bond vectors, if the model fixes it.

Type:

float, optional

contour_length_per_residue

The contour length per residue, if the model is finitely extensible: no two beads k residues apart can be further apart than k times this.

Type:

float, optional

reference_rg

A radius of gyration the model reports itself, shown as context.

Type:

float, optional

reference_rg_label

What reference_rg is, for display.

Type:

str

reference_re

A whole-chain end-to-end distance the model reports itself, shown as context.

Type:

float, optional

reference_re_label

What reference_re is, for display.

Type:

str

end_to_end_distribution

The model’s analytical end-to-end distribution (distances, probabilities) for the first-to-last pair (N - 1 residues apart), shown as context.

Type:

tuple of np.ndarray, optional

end_to_end_label

What that distribution is, for display.

Type:

str

notes

Model-specific caveats printed with the report.

Type:

tuple of str

afrc.ensemble_report.compare_ensemble_to_gaussian_model(conformations: NDArray[float64], mean_squared_distances: NDArray[float64], model_name: str = 'AFRC', reference_mean_rg: float | None = None, reference_mean_re: float | None = None, max_pair_evaluations: int = 100000000) → EnsembleReport

Check that an ensemble reproduces the statistics of a Gaussian chain model.

A convenience wrapper around compare_ensemble_to_model() for a model in which every inter-residue distance is Gaussian, fully specified by its mean-squared inter-residue distances (as for the AFRC).

Parameters:

• conformations (np.ndarray) – Array of shape [n_conformations x N x 3] with the bead coordinates, in Angstroms.

• mean_squared_distances (np.ndarray) – The model’s symmetric [N x N] matrix of mean-squared inter-residue distances, in Angstroms squared.

• model_name (str) – Name of the model, for display. Default is ‘AFRC’.

• reference_mean_rg (float, optional) – The model’s own mean radius of gyration, if it defines one separately (the AFRC’s is calibrated independently of its distances). Reported as context.

• reference_mean_re (float, optional) – The model’s whole-chain mean end-to-end distance (the AFRC’s uses N rather than N-1 residues). Reported as context.

• max_pair_evaluations (int) – Upper limit on conformations x residue pairs for the pair-based checks.

Returns:

The checks and context comparisons.

Return type:

EnsembleReport

Raises:

AFRCException – If the conformations and distance matrix do not have matching shapes.

afrc.ensemble_report.compare_ensemble_to_model(conformations: NDArray[float64], expectations: ModelExpectations, max_pair_evaluations: int = 100000000) → EnsembleReport

Check that an ensemble reproduces the statistics of the model it came from.

For every model the ensemble is checked against the exact mean-squared inter-residue distances \(D_{ij}\):

• the root-mean-square radius of gyration, \(\langle R_g^2 \rangle = \frac{1}{N^2}\sum_{i<j} D_{ij}\);

• the first-to-last bead distance;

• the mean-squared distance of every residue pair, in bands of sequence separation.

For a Gaussian model it is also checked against the Maxwell distribution each distance follows: the mean distances (again in bands), the Kirkwood-Riseman hydrodynamic radius (\(\langle 1/r_{ij} \rangle = \sqrt{6/(\pi D_{ij})}\)), and the full distribution of three representative distances (Kolmogorov-Smirnov test). Models with fixed bond lengths, fixed bond angles or a maximum extension get exact structural checks too.

Parameters:

• conformations (np.ndarray) – Array of shape [n_conformations x N x 3] with the bead coordinates, in Angstroms.

• expectations (ModelExpectations) – What the model says the ensemble should look like.

• max_pair_evaluations (int) – Upper limit on conformations x residue pairs for the pair-based checks. If the ensemble exceeds it, those checks use the first conformations only.

Returns:

The checks, context comparisons and notes.

Return type:

EnsembleReport

Raises:

AFRCException – If the conformations and distance matrix do not have matching shapes.