Sequence-Dependent Collapse of a Block Copolymer (OpenMM)

Same beads, different arrangement¶

Take $N = 50$ beads, half A and half B. A homopolymer of all-A (strong self-attraction) collapses into a compact globule. An alternating ABABAB sequence behaves like a weak homopolymer — the average pair interaction is the arithmetic mean of $\varepsilon_{AA}, \varepsilon_{BB}, \varepsilon_{AB}$ . A block copolymer AAA...BBB... forms a two-lobed structure: two collapsed globules joined at the interface. Random sequences interpolate between these.

This is the coarse-grained picture behind stickers-and-spacers models of intrinsically disordered proteins and biomolecular condensates: sequence is encoded in the pattern of attractive beads.

# Only install conda + OpenMM when running on Google Colab.
# On a local/CI kernel this cell is a no-op.
try:
    import google.colab  # noqa: F401
    _IN_COLAB = True
except ImportError:
    _IN_COLAB = False

if _IN_COLAB:
    get_ipython().system('pip install -q condacolab')
    import condacolab
    condacolab.install()

%%capture
!conda install -c conda-forge openmm mdtraj parmed
!pip install py3dmol

import numpy as np
import matplotlib.pyplot as plt
import openmm as mm
import openmm.app as app
import openmm.unit as unit
from sys import stdout

Build a polymer with a specified AB sequence¶

We use CustomNonbondedForce with per-particle type labels so the same object handles all three pair types. Bonded neighbors are excluded from the nonbonded force and connected by a harmonic spring.

def make_polymer_system(sequence, eps_AA=1.0, eps_BB=1.0, eps_AB=0.1,
                        sigma=1.0, k_bond=100.0, r0=1.0, box=20.0):
    # sequence is a string of 'A'/'B' characters; its length sets the chain length
    N = len(sequence)
    system = mm.System()
    system.setDefaultPeriodicBoxVectors([box, 0, 0], [0, box, 0], [0, 0, box])
    for _ in range(N):
        system.addParticle(1.0)

    # Custom nonbonded: type-dependent epsilon via per-particle parameter a (1=A, 0=B)
    expr = (
        "4*eps*((sig/r)^12 - (sig/r)^6); "
        "eps = eps_AA*a1*a2 + eps_BB*(1-a1)*(1-a2) + eps_AB*(a1*(1-a2) + (1-a1)*a2)"
    )
    nb = mm.CustomNonbondedForce(expr)
    nb.addGlobalParameter('eps_AA', eps_AA)
    nb.addGlobalParameter('eps_BB', eps_BB)
    nb.addGlobalParameter('eps_AB', eps_AB)
    nb.addGlobalParameter('sig', sigma)
    nb.addPerParticleParameter('a')
    for s in sequence:
        nb.addParticle([1.0 if s == 'A' else 0.0])
    nb.setCutoffDistance(3.0 * sigma)
    nb.setNonbondedMethod(mm.CustomNonbondedForce.CutoffNonPeriodic)

    # Harmonic bonds + exclusions
    bond = mm.HarmonicBondForce()
    for i in range(N - 1):
        bond.addBond(i, i + 1, r0, k_bond)
        nb.addExclusion(i, i + 1)
    system.addForce(nb)
    system.addForce(bond)

    # Initial positions: stretched chain along x
    pos = np.zeros((N, 3))
    pos[:, 0] = np.arange(N) * r0
    pos -= pos.mean(axis=0)
    return system, pos


def run_polymer(sequence, T=0.7, gamma=1.0, dt=0.005, n_steps=200_000,
                save_every=500, **kwargs):
    system, pos = make_polymer_system(sequence, **kwargs)
    integrator = mm.LangevinIntegrator(T, gamma, dt)
    context = mm.Context(system, integrator)
    context.setPositions(pos)
    context.setVelocitiesToTemperature(T)
    n_save = n_steps // save_every
    traj = np.empty((n_save, len(sequence), 3))
    for k in range(n_save):
        integrator.step(save_every)
        traj[k] = context.getState(getPositions=True).getPositions(asNumpy=True)
    return traj

Run three sequences at the same composition¶

50% A, 50% B in three arrangements. Same code, same $T$ , same interaction parameters — only the sequence changes.

sequences = {
    'alternating': 'AB' * 25,
    'diblock':     'A' * 25 + 'B' * 25,
    'random':      ''.join(np.random.default_rng(3).choice(['A', 'B'], size=50)),
}

trajs = {name: run_polymer(seq, T=0.7, eps_AA=1.0, eps_BB=1.0, eps_AB=0.1,
                           n_steps=100_000) for name, seq in sequences.items()}


def R_g(coords):
    com = coords.mean(axis=-2, keepdims=True)
    return np.sqrt(((coords - com) ** 2).sum(axis=-1).mean(axis=-1))


fig, ax = plt.subplots(figsize=(8, 4))
for name, traj in trajs.items():
    rg = R_g(traj)
    ax.plot(rg, lw=1, label=f'{name}: $\\langle R_g\\rangle = {rg[len(rg)//2:].mean():.2f}$')
ax.set(xlabel='frame', ylabel='$R_g$',
       title='Radius of gyration vs. time for three sequences')
ax.legend()
ax.grid(True, ls=':', alpha=0.5)
fig.tight_layout()

Problems¶

Collapse transition vs. sequence. For each of the three sequences above, run at $T \in \{0.3, 0.5, 0.7, 1.0, 1.5\}$ and plot $\langle R_g(T) \rangle$ . Identify the approximate collapse temperature $T_\theta$ for each sequence. Why does the diblock collapse at higher $T$ than the alternating chain?
Contact maps. Define a contact at $r < 1.5\,\sigma$ . Compute the average contact map $C_{ij} = \langle \Theta(1.5\sigma - |r_i - r_j|)\rangle$ for the diblock at $T = 0.5$ . You should see two compact square blocks on the diagonal (intra-block contacts) and a sparse off-diagonal strip (interface contacts).
Patterned sequence. Design a patterned sequence of your choice (e.g. (ABBA)*12 or stickers-and-spacers (AXXXXA)...). Compare its $R_g(T)$ and contact map to the three reference sequences. Can you design a sequence that collapses more compactly than the homopolymer of the same composition? If not, why not?

Reference¶

A. R. Khokhlov, P. G. Khalatur, Conformation-dependent sequence design (engineering) of AB copolymers, Phys. Rev. Lett. 82, 3456 (1999). Khokhlov & Khalatur (1999)

References¶

Khokhlov, A. R., & Khalatur, P. G. (1999). Conformation-Dependent Sequence Design (Engineering) of<mml:math xmlns:mml=“http://www.w3.org/1998/Math/MathML” display=“inline”><mml:mi mathvariant=“italic”>AB</mml:mi></mml:math>Copolymers. Physical Review Letters, 82(17), 3456–3459. 10.1103/physrevlett.82.3456