Source code for arbdmodel.onck_polymer_model

# -*- coding: utf-8 -*-
## Test with `python -m arbdmodel.onck_polymer_model`

import numpy as np
import sys

## Local imports
from .logger import devlogger, logger, get_resource_path
from .core_objects import ParticleType, PointParticle
from .polymer import PolymerBeads, PolymerModel
from .interactions import AbstractPotential, HarmonicBond
from .version import Citation

"""Define particle types"""
_types = dict(
    A = ParticleType("ALA",
                     mass = 120,
                     charge = 0,
                     epsilon = 0.7, # dimensionless
                     lambda_ = 0.72973,
                 ),
    R = ParticleType("ARG",
                     mass = 120,
                     charge = 1,
                     epsilon = 0.0,
                     lambda_ = 0.0,
                 ),
    N = ParticleType("ASN",
                     mass = 120,
                     charge = 0,
                     epsilon = 0.33,
                     lambda_ = 0.432432,
                 ),
    D = ParticleType("ASP",
                     mass = 120,
                     charge = -1,
                     epsilon = 0.0005,
                 ),
    C = ParticleType("CYS",
                     mass = 120,
                     charge = 0,
                     epsilon = 0.68,
                 ),
    Q = ParticleType("GLN",
                     mass = 120,
                     charge = 0,
                     epsilon = 0.64,
                 ),
    E = ParticleType("GLU",
                     mass = 120,
                     charge = -1,
                     epsilon = 0.0005,
                 ),
    G = ParticleType("GLY",
                     mass = 120,
                     charge = 0,
                     epsilon = 0.41,
                 ),
    H = ParticleType("HIS",
                     mass = 120,
                     charge = 0.0,
                     epsilon = 0.53,
                 ),
    I = ParticleType("ILE",
                     mass = 120,
                     charge = 0,
                     epsilon = 0.98,
                 ),
    L = ParticleType("LEU",
                     mass = 120,
                     charge = 0,
                     epsilon = 1.0,
                 ),
    K = ParticleType("LYS",
                     mass = 120,
                     charge = 1,
                     epsilon = 0.0005,
                 ),
    M = ParticleType("MET",
                     mass = 120,
                     charge = 0,
                     epsilon = 0.78,
                 ),
    F = ParticleType("PHE",
                     mass = 120,
                     charge = 0,
                     epsilon = 1.0,
                 ),
    P = ParticleType("PRO",
                     mass = 120,
                     charge = 0,
                     epsilon = 0.65,
                 ),
    S = ParticleType("SER",
                     mass = 120,
                     charge = 0,
                     epsilon = 0.45,
                 ),
    T = ParticleType("THR",
                     mass = 120,
                     charge = 0,
                     epsilon = 0.51,
                 ),
    W = ParticleType("TRP",
                     mass = 120,
                     charge = 0,
                     epsilon = 0.96,
                 ),
    Y = ParticleType("TYR",
                     mass = 120,
                     charge = 0,
                     epsilon = 0.82,
                 ),
    V = ParticleType("VAL",
                     mass = 120,
                     charge = 0,
                     epsilon = 0.94,
                 )
)
for k,t in _types.items():
    t.resname = t.name
    t.version = 1.0

_types_versions = {1.0: _types}

## https://link.springer.com/article/10.1007/s12274-022-4647-1
_types_versions['1.0cp'] = {k:ParticleType(t.name,
                                           mass = t.mass,
                                           charge = t.charge,
                                           epsilon = t.epsilon,
                                           version = '1.0cp') for k,t in _types.items()}
for k in 'R D E K'.split(): _types_versions['1.0cp'][k].epsilon = 0.005

_types_versions[1.1] = dict()   # https://www.pnas.org/doi/10.1073/pnas.2221804120
__si_data = """A 71.07 0.7 0
L 113.16 1 0
R 156.19 0.005 +1
K 128.17 0.005 +1
N 114.10 0.41 0
M 131.20 0.78 0
D 115.09 0.005 -1
F 147.18 1 0
C 103.14 0.68 0
P 97.12 0.65 0
Q 128.13 0.33 0
S 87.08 0.45 0
E 129.11 0.005 -1
T 101.11 0.51 0
G 57.05 0.48 0
W 186.22 0.96 0
H 137.14 0.53 0
Y 163.18 0.82 0
I 113.16 0.98 0
V 99.13 0.94 0"""
for k, m, e, q in [l.split() for l in __si_data.split('\n')]:
    t0 = _types_versions[1.0][k]
    t  = _types_versions[1.1][k] = ParticleType(t0.name+'v1.1', mass=float(m),
                                                epsilon = float(e), charge = float(q), version=1.1)
    assert(t.charge == t0.charge)
    if t0.epsilon != t.epsilon:
        devlogger.info(f'Updating epsilon for 1BPA-1.1 {k}: {t0.epsilon} -> {t.epsilon}')

devlogger.info(f'Onck model version 1BPA-1.1 has {np.sum([t0.epsilon != _types_versions[1.1][k].epsilon for k,t0 in _types_versions[1.0].items()])} different epsilon parameters compared to 1BPA')

""" Cation-Pi interactions """
eps_cp_dict = {('ARG','PHE'): 4.3,
               ('ARG','TYR'): 5.0,
               ('ARG','TRP'): 6.7,
               ('LYS','PHE'): 1.79,
               ('LYS','TYR'): 3.13,
               ('LYS','TRP'): 4.26} # Table S2; kJ/mol
eps_cp_dict = {k:v*0.23900574 for k,v in eps_cp_dict.items()} # convert to kcal/mol
for k,v in list(eps_cp_dict.items()): eps_cp_dict[(k[1],k[0])] = v # add reversed keys to dictionary

version_name = {1.0:'1BPA', '1.0cp':'1BPA-CP', 1.1:'1BPA-1.1'}

__base_ref = Citation(
    author  = 'A. Fragasso, H.W. de Vries, J. Andersson, et al',
    title   = 'A designer FG-Nup that reconstitutes the selective transport barrier of the nuclear pore complex',
    journal = 'Nat Commun',
    volume  = 12,
    year    = 2021,
    doi     = '10.1038/s41467-021-22293-y'
)
__ref1 = Citation(
    author  = 'A. Fragasso, H.W. de Vries, J. Andersson, et al',
    title = 'Transport receptor occupancy in nuclear pore complex mimics',
    journal = 'Nano Res',
    volume = 15,
    pages = '9689–9703',
    year = 2022,
    doi = '10.1007/s12274-022-4647-1'
)
__ref2 = Citation(
    author="M. Dekker, E. Van der Giessen, and P.R. Onck",
    title="Phase separation of intrinsically disordered FG-Nups is driven by highly dynamic FG motifs",
    journal='PNAS',
    volume=120,
    number = 25,
    pages = 'e2221804120',
    year=2023,
    doi = '10.1073/pnas.2221804120'
)

version_refs = {1.0: (__base_ref,),
                '1.0cp': (__base_ref, __ref1),
                1.1: (__base_ref, __ref2),
                }

## Default to the latest
_types = _types_versions[1.1]

[docs] class OnckNonbonded(AbstractPotential): """ Nonbonded interaction potential for 1BPA and 1BPA-1.1. This class implements nonbonded interactions including electrostatics with distance-dependent dielectric and Lennard-Jones type potentials. Parameters ---------- debye_length : float, optional The Debye-Hückel screening length in angstroms. Default is 12.7 Å. resolution : float, optional The spatial resolution for potential calculation. Default is 0.1 Å. range_ : tuple, optional The range of distances for which to calculate the potential. Default is (0, None). max_force : float Maximum allowed force in the system, used as a cap for stability. Default is 50 Note ----- The potential includes: 1. Electrostatic interactions with distance-dependent dielectric constant 2. Modified Lennard-Jones interactions that depend on residue types: - Special 8-6 LJ potential for cationic-aromatic interactions in 1.0cp/1.1 versions - General LJ-type potential for other interactions The effective dielectric D decreases with distance according to a specific formula. """ def __init__(self, debye_length=12.7, resolution=0.1, range_=(0,None)): AbstractPotential.__init__(self, resolution=resolution, range_=range_) self.debye_length = debye_length self.max_force = 50
[docs] def potential(self, r, types): """ Electrostatics """ typeA, typeB = types ld = self.debye_length q1 = typeA.charge q2 = typeB.charge D = 80 # dielectric of water _z = 2.5 D = 80 * (1- (r/_z)**2 * np.exp(r/_z)/(np.exp(r/_z)-1)**2) ## units "e**2 / (4 * pi * epsilon0 AA)" kcal_mol A = 332.06371 u_elec = (A*q1*q2/D)*np.exp(-r/ld) / r """ LJ-type term """ try: """ Rather than interacting through the hydrophobic potential phi_hp, cationic residues interact with aromatic residues through an 8–6 Lennard Jones potential """ assert( all( t.version in ('1.0cp',1.1) for t in types ) ) key = tuple((typeA.name[:3],typeB.name[:3])) eps = eps_cp_dict[key] _r_m = 4.5 r6 = (_r_m/r)**6 r8 = (_r_m/r)**8 u_lj = eps * (3*r8 - 4*r6) except: sigma = 6.0 r6 = (sigma/r)**6 r8 = (sigma/r)**8 alpha = 0.27 epsilon_hp = 3.1070746 # units "13 kJ/N_A" kcal_mol epsilon_rep = 2.3900574 # units "10 kJ/N_A" kcal_mol epsilon = epsilon_hp*np.sqrt( (typeA.epsilon*typeB.epsilon)**alpha ) u_lj = (epsilon_rep-epsilon) * r8 s = r<=sigma u_lj[s] = epsilon_rep*r8[s] - epsilon*(4*r6[s]-1)/3 u_lj[r>25] = 0 u = u_elec + u_lj return u
[docs] class OnckBeads(PolymerBeads): """ A class that represents polymer beads in the Onck model. The Onck model is a coarse-grained model for proteins where each amino acid is represented by a single bead. The model includes specific bonded interactions (bonds, angles, dihedrals) based on the amino acid types in the sequence. Parameters ---------- polymer : Polymer The parent polymer object to which these beads belong. sequence : list, optional The sequence of amino acid types for each bead, must match polymer length. spring_constant : float, default=38.422562 Spring constant for peptide bonds in units of 8038 kJ/(N_A nm^2) or 0.5 kcal_mol/AA^2. rest_length : float, default=3.8 Equilibrium length of peptide bonds. version : float, optional Version of the Onck model to use. Defaults to 1.1 if not specified. **kwargs Additional keyword arguments passed to the parent PolymerBeads constructor. types_dict : dict Mapping between amino acid codes and bead types. peptide_bond : HarmonicBond Bond potential used for peptide bonds between adjacent beads. Raises ------ ValueError If the version is not supported or if the sequence length doesn't match the polymer length. NotImplementedError If sequence is None (random sequence generation not implemented). Note ----- The Onck model distinguishes between proline (P), glycine (G), and all other amino acids (X) for certain bonded interactions. Special potentials are loaded from files for angle and dihedral interactions based on this classification. """ def __init__(self, polymer, sequence=None, spring_constant = 38.422562, # units "8038 kJ / (N_A nm**2)" "0.5 * kcal_mol/AA**2" rest_length=3.8, version=None, **kwargs): if version == None: logger.warning(f'No Onck model version specified; using version 1BPA-1.1 from 2023') version = 1.1 if version not in _types_versions: raise ValueError(f'Unkown Onck model version "{version}"') self.version = version self.types_dict = _types_versions[version] _types = _types_versions[version] # Update global _types convenience variable self.peptide_bond = HarmonicBond(k = spring_constant, r0 = rest_length, range_ = (0,100), resolution = 0.01, max_force = 10) if sequence is None: raise NotImplementedError # ... set random sequence self.polymer = polymer self.sequence = sequence self.spring_constant = spring_constant PolymerBeads.__init__(self, polymer, sequence, rest_length=rest_length, **kwargs) assert(self.monomers_per_bead_group == 1) if len(sequence) != polymer.num_monomers: raise ValueError("Length of sequence does not match length of polymer") def _generate_ith_bead_group(self, i, r, o): s = self.sequence[i] return PointParticle(self.types_dict[s], r, name = s, resid = i+1) def _join_adjacent_bead_groups(self, ids): def bead_to_type(bead): if bead.type_.name == 'PRO': return 'P' elif bead.type_.name == 'GLY': return 'G' else: return 'X' ## Two consecutive groups if len(ids) == 2: b1,b2 = [self.children[i] for i in ids] """ units "10 kJ/N_A" kcal_mol """ bond = self.peptide_bond self.add_bond( i=b1, j=b2, bond = bond, exclude=True ) elif len(ids) == 3: b1,b2,b3 = [self.children[i] for i in ids] t1,t2,t3 = [bead_to_type(b) for b in (b1,b2,b3)] filename = 'onck_model_potentials/bend_O{}{}.txt'.format( t2, 'P' if t3 == 'P' else 'Y' ) self.add_angle( i=b1, j=b2, k=b3, angle = get_resource_path(filename) ) self.add_exclusion( i=b1, j=b3 ) elif len(ids) == 4: ## Four consecutive monomers b1,b2,b3,b4 = [self.children[i] for i in ids] t1,t2,t3,t4 = [bead_to_type(b) for b in (b1,b2,b3,b4)] filename = 'onck_model_potentials/dih_{}{}.txt'.format(t2,t3) self.add_dihedral( i=b1, j=b2, k=b3, l=b4, dihedral = get_resource_path(filename) ) self.add_exclusion( i=b1, j=b4 ) else: raise Exception('Programming error!')
[docs] class OnckModel(PolymerModel): """ A polymer model based on Onck's coarse-grained model for disordered FG-Nup peptides. The OnckModel class implements a coarse-grained model of polymers based on the work by Onck et al. This model is specifically designed for disordered FG-Nup peptides and incorporates electrostatic interactions through a Debye-Hückel potential. Parameters ---------- polymers : list List of polymers to be modeled. sequences : list, optional List of amino acid sequences corresponding to each polymer. Must be provided. Default is None. debye_length : float, optional Debye screening length in angstroms. Default is 10. damping_coefficient : float, optional Damping coefficient in ns. Default is 100. version : float, optional Version of the Onck model to use. Default is 1.1 (corresponding to 1BPA-1.1 from 2023). DEBUG : bool, optional Whether to enable debug mode. Default is False. **kwargs : dict Additional keyword arguments to pass to the parent PolymerModel class. Common options include: - timestep: Simulation timestep (default: 20e-6) - cutoff: Cutoff distance for non-bonded interactions (default: max(5*debye_length, 25)) - decomp_period: Decomposition period (default: 1000) types_dict : dict Dictionary mapping from type keys to type objects for the given version. types : list List of all types used in the model. Note ----- The model differs from the published Onck model if a Debye length other than 12.7 Å is chosen. In such cases, the non-bonded cutoff is set to 5 * debye_length. References ---------- Please refer to the citations displayed during initialization. """ def __init__(self, polymers, sequences = None, debye_length = 10, # damping_coefficient = 50e3, damping_coefficient = 100, version = None, DEBUG=False, **kwargs): if version is None: logger.warning(f'No Onck model version specified; using version 1BPA-1.1 from 2023') version = 1.1 if version not in _types_versions: raise ValueError(f'Unkown Onck model version "{version}"') self.version = version self.types_dict = _types_versions[version] logger.info(f'Using an implementation of the Onck model {version_name[self.version]} for disordered FG-Nup peptides.') _msg = 'Please cite all appropriate articles, including:\n' _msg = _msg + "\n and\n".join(ref.display() for ref in version_refs[self.version]) print(_msg) if debye_length != 12.7: logger.warning("""Deviating from the model published by Onck by choosing a Debye length differing from 1.27 nm. Be advised that the non-bonded cutoff is simply set to 5 * debye_length, but this is not necessarily prescribed by the model.""") if 'timestep' not in kwargs: kwargs['timestep'] = 20e-6 if 'cutoff' not in kwargs: kwargs['cutoff'] = max(5*debye_length,25) if 'decomp_period' not in kwargs: kwargs['decomp_period'] = 1000 """ Assign sequences """ if sequences is None: raise NotImplementedError("OnckModel must be provided a sequences argument") PolymerModel.__init__(self, polymers, sequences, monomers_per_bead_group=1, **kwargs) """ Update type diffusion coefficients """ self.types = all_types = [t for key,t in self.types_dict.items()] self.set_damping_coefficient( damping_coefficient ) """ Set up nonbonded interactions """ nonbonded = OnckNonbonded(debye_length) for t in all_types: self._add_nonbonded_interaction(nonbonded, t) def _add_nonbonded_interaction(self, interaction, type_): i = self.types.index(type_) if type_ in self.types else 0 for t in self.types[i:]: # self.add_nonbonded_interaction( interaction, typeA=type_, typeB=t ) self.add_nonbonded_interaction( interaction, typeA=type_, typeB=t ) def _generate_polymer_beads(self, polymer, sequence, polymer_index=None): return OnckBeads(polymer, sequence, monomers_per_bead_group = self.monomers_per_bead_group, poymer_index = polymer_index, version = self.version )
[docs] def set_damping_coefficient(self, damping_coefficient): for t in self.types: t.damping_coefficient = damping_coefficient
# t.diffusivity = 831447.2 * temperature / (t.mass * damping_coefficient) if __name__ == "__main__": logger.info("TYPES") for n,t in _types.items(): logger.info("{}\t{}\t{}\t{}\t{}".format(n, t.name, t.mass, t.charge, t.epsilon))