API reference

This module provides an implementation of the Rouse model.

Its content is mostly the Model class, so see there for more description. Otherwise, there is the functional form for the MSD of two loci on an infinite continuous Rouse polymer in twoLocusMSD. This might be used as phenomenological expression for fitting to data.

See also

Model

class rouse.Model(N, D=1.0, k=1.0, d=3, setup_dynamics=True, add_bonds=None)

The Rouse model

This model follows the linear Langevin equation

\[\dot{x}(t) = -kAx(t) + F + \xi(t)\,,\]

where \(x\) is a vector of monomer positions and \(\xi\) is Gaussian white noise of strength proportional to D driving the system. The connectivity matrix \(A\) describes which monomers are connected to which by harmonic bonds with spring constant \(k\). The connectivity matrix is conveniently set up using setup_free_chain and add_bonds. Finally, \(F\) can be used to incorporate an external force on individual monomers.

This model implementation provides the following functionality:

  • setting up a spring network with arbitrary connectivity matrix

  • calculating the steady state ensemble.

  • propagating an arbitrary Gaussian ensemble forward in time

  • evolving a single conformation \(x\) forward in time

  • analytically calculating MSD & ACF, also for arbitrary linear combinations of loci (e.g. MSD of the vector between two loci)

  • other auxiliary stuff, like giving equilibration time scales and contact probabilities

N

number of monomers

Type:

int

D

diffusivity of an individual monomer (i.e. strength of the driving white noise)

Type:

float

k

default spring constant for harmonic bonds

Type:

float

d

spatial dimension

Type:

int

A

connectivity matrix. Initialized to a linear polymer by default

Type:

(N, N) array, dtype=float, internal

F

deterministic external force on each monomer. Zero by default

Type:

(N, d) array, dtype=float, internal

Parameters:

  • setup_dynamics (bool) – whether to pre-calculate matrices needed for dynamics (propagation of ensembles / evolution of single conformation).

  • add_bonds (list of bonds) – shortcut for adding bonds to the system before initializing dynamics. Corresponds to the links argument of add_bonds, see there for more details.

setup_free_chain(d=3)

Set up connectivity matrix of a linear chain and zero external force

Mostly used internally upon initialization

Parameters:

d (int) – spatial dimension

See also

add_bonds, add_tether

add_bonds(links, k_rel=1.0)

Add additional bonds to connectivity matrix

Parameters:

  • links (list of (monA, monB, k_rel) or (monA, monB) tuples) – bonds to add. Introduces an additional bond between monA and monB (integer indices < N) of strength k_rel*k. Giving k_rel as third entry in the tuple overwrites use of the function argument k_rel for this bond.

  • k_rel (float > 0, optional) – strength of the new bonds, as multiple of the class attribute k.

See also

add_tether

add_tether(mon=0, k_rel=1.0, point=None)

Tether a monomer to a fixed point in space

This simply introduces an additional harmonic bond between this monomer and the specified point. Note that you can give a large k_rel to actually pin the monomer in place

Parameters:

  • mon (int) – the index of the monomer that should get the tether

  • k_rel (float > 0, optional) – strength of the tether, in multiples of the default bond strength k.

  • point ((d,) array-like, optional) – the point to tether to; defaults to the origin

See also

add_bonds

update_dynamics(dt=1.0)

Pre-calculate dynamic matrices

This sets the model up for evaluation of model dynamics, such as evolving a conformation or propagating an ensemble. It should be called after any modification to model parameters, which is handled automatically for internal modifications (such as when add_bonds is called). If you modify anything externally, recommended behavior is to simply set model._dynamics['needs_updating'] = True, which will ensure that stuff is recalculated as needed.

Ultimately, this is just a wrapper for the update_* functions. If you know what you’re doing, you can also just call these individually as needed.

Parameters:

dt (float) – time step to prepare for

See also

check_dynamics

update_F_only(override_full_update=False)

Update specifically the stuff depending on the external force F.

Parameters:

override_full_update (bool) – pretend that this was a full update. Useful when you changed only F since the last update, to avoid recalculating everything else

See also

update_dynamics

check_dynamics(dt=None, run_if_necessary=True)

Check that dynamics are set up properly

Parameters:

  • dt (float or None) – the time step we should be set up for. Specify None to accept whatever was the value set up

  • run_if_necessary (bool) – controls the behavior when dynamics are not set up. If set to True (default), update_dynamics is simply run with the settings we have. If set to False raise a RuntimeError if dynamics are not (correctly) set up for the current model.

Raises:

RuntimeError – if something is not set up properly

See also

update_dynamics

steady_state()

Calculate the steady state ensemble

Since our model is linear and driven by Gaussian noise, the steady state ensemble will always be Gaussian, and as such given by a mean vector and covariance matrix over all monomers.

We use the Moore-Penrose pseudo-inverse of the connectivity matrix A to calculate the steady state distribution. This means that any degrees of freedom (such as center of mass for a free chain) will be pinned to the origin instead of diverging to infinity.

Returns:

  • M ((N, d) np.ndarray, dtype=float) – mean position of each monomer

  • C ((N, N) np.ndarray, dtype=float) – covariance between monomer coordinates in each dimension. Note that we assume independence of spatial dimensions, such that the full covariance C_ijab := <x_ia*x_jb>_c (with monomer indices i, j and spatial indices a, b) is written as the tensor product C_ijab = C_ij*δ_ab with Kronecker’s δ. The matrix returned by this function is C_ij.

propagate_M(M, dt=None, check_dynamics=True)

Propagate the mean of a Gaussian ensemble

Parameters:

  • M ((N, d) np.ndarray, dtype=float) – mean conformation

  • dt (float or None, optional) – the time step to propagate for. When None uses the time step set through update_dynamics (default, recommended)

  • check_dynamics (bool) – whether to check for correct setup of precalculated matrices. Can be disabled for performance when its otherwise clear that the setup is correct

Returns:

M ((N, d) np.ndarray, dtype=float) – like input M, but propagated by dt.

See also

propagate, propagate_C, steady_state

propagate_C(C, dt=None, check_dynamics=True)

Propagate the covariance of a Gaussian ensemble

Parameters:

  • C ((N, N) np.ndarray, dtype=float) – current covariance

  • dt (float or None, optional) – the time step to propagate for. When None uses the time step set through update_dynamics (default, recommended)

  • check_dynamics (bool) – whether to check for correct setup of precalculated matrices. Can be disabled for performance when its otherwise clear that the setup is correct

Returns:

C ((N, N) np.ndarray, dtype=float) – like input C, but propagated by dt.

See also

propagate, propagate_M, steady_state

propagate(M, C, dt=None, check_dynamics=True)

Propagate a Gaussian ensemble one time step

This is simply a wrapper for propagate_M and propagate_C.

Parameters:

  • M ((N, d) np.ndarray, dtype=float) – current mean conformation

  • C ((N, N) np.ndarray, dtype=float) – current covariance

  • dt (float or None, optional) – the time step to propagate for. When None uses the time step set through update_dynamics (default, recommended)

  • check_dynamics (bool) – whether to check for correct setup of precalculated matrices. Can be disabled for performance when its otherwise clear that the setup is correct

Returns:

M, C – like input, but propagated by dt.

See also

propagate_M, propagate_C, steady_state

conf_ss()

Draw a conformation from steady state

Returns:

(N, d) np.ndarray, dtype=float – the drawn conformation

See also

steady_state

evolve(conf, dt=None, check_dynamics=True)

Evolve a conformation forward one time step

Parameters:

  • conf ((N, d) np.ndarray, dtype=float) – the conformation to start from

  • dt (float or None, optional) – the time step to propagate for. When None uses the time step set through update_dynamics (default, recommended)

  • check_dynamics (bool) – whether to check for correct setup of precalculated matrices. Can be disabled for performance when its otherwise clear that the setup is correct

Returns:

conf – like input, but evolved for dt.

See also

conf_ss

MSD(dts, w=None)

Calculate MSD for given degrees of freedom

This is simply an evaluation of analytical expressions. Which path exactly we take depends on the system, see Notes.

Parameters:

  • dts ((T,) array-like, dtype=float) – the time lags to evaluate the MSD at. Entries should be greater or equal to zero

  • w ((N,) np.ndarray, dtype=float) – measurement vector. Use this to specify the linear combination of monomers whose MSD you are interested in. See Examples. If unspecified, the function returns the full covariance matrix at lag Δt.

Returns:

(T,) np.ndarray or (T, N, N) np.array – the MSD, either as scalar function evaluated at dts or matrix valued (if no measurement vector w was specified).

Notes

The analytical solution for the MSD at time lag \(\Delta t\) is given by

\[MSD(\Delta t) \equiv \left\langle\left(x(t+\Delta t) - x(t)\right)^{\otimes 2}\right\rangle = V\left[ 2D E + (\Delta t)^2 SV^TF F^TVS \right] V^T\]

with \(E = \int_0^{\Delta t}\mathrm{d}\tau\exp(-kA\tau)\), S_ij = 1 if i == j and λ_i = 0 else 0, and \(A = V\Lambda V.T\) is the eigendecomposition of A. Note that the Δt^2 terms just contributes ballistic motion for force components acting on singular degrees of freedom (e.g. center of mass).

Examples

>>> model = Model(N=25)
... dts = np.array([0, 1, 2, 10]) # some sample lag times
...
... # Calculate full matrices. This is usually not as relevant
... msd = model.MSD(dts)
...
... # MSD of the middle monomer
... w = np.zeros(model.N)
... w[12] = 1
... msd = model.MSD(dts, w)
...
... # MSD of the end-to-end vector
... w = np.zeros(model.N)
... w[0]  =  1
... w[-1] = -1
... msd = model.MSD(dts, w)
...
... # The end-to-end dof equilibrates, so we can calculate the limit
... # value
... msd_inf = model.MSD([np.inf], w)
ACF(dts, w=None)

Calculate autocovariances

Exists mostly for completeness. Calculates the autocovariance as γ = exp(-k*A*dt) @ C, where C is the steady state covariance as returned by steady_state(). Note that for any equilibrating degrees of freedom the ACF can also be calculated from the MSD as

ACF(Δt) = 0.5*( MSD(∞) - MSD(Δt) )

so sticking to the MSD is preferred. Any non-equilibrating dofs are pinned to zero by the Moore-Penrose inverse.

The one benefit of this function over MSD is that it is faster (roughly 1.5x).

Parameters:

  • dts ((T,) array-like, dtype=float) – the time lags to evaluate the ACF at. Entries should be greater or equal to zero

  • w ((N,) np.ndarray, dtype=float) – measurement vector. Use this to specify the linear combination of monomers whose ACF you are interested in. If unspecified, the function returns the full covariance matrix at lag Δt.

Returns:

(T,) np.ndarray or (T, N, N) np.array – the ACF, either as scalar function evaluated at dts or matrix valued (if no measurement vector w was specified).

See also

MSD

timescales()

Give time scales relevant to the model

Returns:

dict – t_microscopic is where Rouse scaling (i.e. the model) starts to apply; t_Rouse is the classic Rouse time (relaxation of the slowest mode); t_equilibration is the actual crossover from Rouse to equilibrium of the end-to-end distance (for a free chain)

Gamma()

MSD prefactor of a single locus on the polymer

Returns:

float

rms_Ree(L=None)

RMS end-to-end distance for chain of length L

Parameters:

L (float, optional) – number of (effective) bonds in the chain. Defaults to self.N-1.

Returns:

float

contact_frequency()

Calculate a contact frequency matrix for the system

This is intended to produce HiC-like maps for the system. It is based on the mean squared distance between each pair of loci, and then converts to contact frequency by a scaling exponent of -self.d/2. Note that these are not actual probabilities, but relative frequencies

Returns:

(N, N) np.ndarray, dtype=float

rouse.twoLocusMSD(dts, Gamma, J)

MSD for relative position of two loci on infinite continuous Rouse polymer

The analytical expression is given by

\[MSD(\Delta t) = 2\Gamma\sqrt{\Delta t}\left(1 - \exp\left(-\frac{\tau}{\Delta t}\right)\right) + 2J\mathrm{erfc}\sqrt{\frac{\tau}{\Delta t}}\,,\]

with \(\tau\equiv \frac{1}{\pi}\left(\frac{J}{\Gamma}\right)^2\).

Parameters:

  • dts (array-like, dtype=float) – the times at which to evaluate the MSD

  • Gamma (float) – the prefactor for the polymer part (0.5 scaling) of the MSD. Note that this parameter should be the prefactor for tracking of one locus, i.e. the MSD produced by this function (which is the distance between two loci) will have a prefactor of 2*Gamma.

  • J (float) – the mean squared distance of the two loci at steady state. This is half the asymptotic value of the MSD: MSD(Δt --> inf) = 2*J.

Returns:

(T,) np.array – the MSD evaluated at times dts