Source code for arbdmodel.mpipi_polymer

## Test with `python -m arbdmodel.mpipi_polymer_model`

import numpy as np
import sys
import pandas as pd
## Local imports
from .logger import logger, get_resource_path
from . import ParticleType, PointParticle

from .polymer import PolymerBeads, PolymerModel
from .interactions import AbstractPotential, HarmonicBond

"""Define particle types"""
_types = dict(
    A = ParticleType("ALA",
                     mass = 71.08,
                     charge = 0,
                     sigma = 5.04,
                     freq=0.091,
                 ),
    R = ParticleType("ARG",
                     mass = 156.2,
                     charge = 1,
                     sigma = 6.56,
                     freq=0.552
                 ),
    N = ParticleType("ASN",
                     mass = 114.1,
                     charge = 0,
                     sigma = 5.68,
                     freq=0.353
                 ),
    D = ParticleType("ASP",
                     mass = 115.1,
                     charge = -1,
                     sigma = 5.58,
                     freq=0.195,
                 ),
    C = ParticleType("CYS",
                     mass = 103.1,
                     charge = 0,
                     sigma = 5.48,
                     freq=0.127,
                 ),
    Q = ParticleType("GLN",
                     mass = 128.1,
                     charge = 0,
                     sigma = 6.02,
                     freq=0.365,
                 ),
    E = ParticleType("GLU",
                     mass = 129.1,
                     charge = -1,
                     sigma = 5.92,
                     freq=0.211,
                 ),
    G = ParticleType("GLY",
                     mass = 57.05,
                     charge = 0,
                     sigma = 4.5,
                     freq=0.22
                 ),
    H = ParticleType("HIS",
                     mass = 137.1,
                     charge = 0.5,
                     sigma = 6.08,
                     freq=0.668
                 ),
    I = ParticleType("ILE",
                     mass = 113.2,
                     charge = 0,
                     sigma = 6.18,
                     freq=0.005
                 ),
    L = ParticleType("LEU",
                     mass = 113.2,
                     charge = 0,
                     sigma = 6.18,
                     freq=0.021
                 ),
    K = ParticleType("LYS",
                     mass = 128.2,
                     charge = 1,
                     sigma = 6.36,
                     freq=0.048
                 ),
    M = ParticleType("MET",
                     mass = 131.2,
                     charge = 0,
                     sigma = 6.18,
                     freq=0.073
                 ),
    F = ParticleType("PHE",
                     mass = 147.2,
                     charge = 0,
                     sigma = 6.36,
                     freq=0.712
                 ),
    P = ParticleType("PRO",
                     mass = 97.12,
                     charge = 0,
                     sigma = 5.56,
                     freq=0.144
                 ),
    S = ParticleType("SER",
                     mass = 87.08,
                     charge = 0,
                     sigma = 5.18,
                     freq=0.113
                 ),
    T = ParticleType("THR",
                     mass = 101.1,
                     charge = 0,
                     sigma = 5.62,
                     freq=0.057,
                 ),
    W = ParticleType("TRP",
                     mass = 186.2,
                     charge = 0,
                     sigma = 6.78,
                     freq=1.0,
                 ),
    Y = ParticleType("TYR",
                     mass = 163.2,
                     charge = 0,
                     sigma = 6.46,
                     freq=0.762
                 ),
    V = ParticleType("VAL",
                     mass = 99.07,
                     charge = 0,
                     sigma = 5.86,
                     freq=0.011
                 )
)

for k,t in list(_types.items()):
    t.resname = t.name
    t.is_idp = False
    
    ## Add types for IDPs

    
[docs] class MpipiNonbonded(AbstractPotential): def __init__(self, debye_length=10, 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*0.75 q2 = typeB.charge*0.75 D = 80 # dielectric of water ## units "e**2 / (4 * pi * epsilon0 AA)" kcal_mol A = 332.06371 u_elec = (A*q1*q2/D)*np.exp(-r/ld) / r """Mpipi Wang-Frenkel Potential""" WF=pd.read_csv(get_resource_path('mpipi_params/mpipi_protein_resname.csv'),index_col=0) indices=list(WF.columns) if indices.index(typeA.resname)<indices.index(typeB.resname): sigma_ij=WF[typeB.resname][typeA.resname+"_sigma"]*10 eps_ij=WF[typeB.resname][typeA.resname+"_eps"]*0.23900574 else: sigma_ij=WF[typeA.resname][typeB.resname+"_sigma"]*10 eps_ij=WF[typeA.resname][typeB.resname+"_eps"]*0.23900574 vij=1 muij=2 if typeA.resname=="ILE": if typeB.resname=="ILE": muij=11 elif typeB.resname=="VAL": muij=4 if typeB.resname=="ILE" and typeA.resname=="VAL": muij=4 Rij=3*sigma_ij #alpha=2*(3**(2*muij))*((2*vij+1/(2*vij*(3**(2*muij)-1))))**(2*vij+1) alpha=2*vij*(Rij/sigma_ij)**(2*muij)*((2*vij+1)/(2*vij*((Rij/sigma_ij)**(2*muij)-1)))**(2*vij+1) u_wang=eps_ij*alpha*((sigma_ij/r)**(2*muij)-1)*((Rij/r)**(2*muij)-1)**(2*vij) """ Mpipi scale model """ """ A_is_idp = B_is_idp = False try: A_is_idp = typeA.is_idp except: pass try: B_is_idp = typeB.is_idp except: pass _idp_scale = (int(A_is_idp)*int(B_is_idp)) alpha = 0.159 + _idp_scale * (0.228 - 0.159) epsilon0 = -1.36 + _idp_scale * (1.36 - 1.0) e_mj = "ERR"[(typeA.resname,typeB.resname)] epsilon = alpha * np.abs( e_mj - epsilon0 ) lambda_ = -1 if epsilon0 > e_mj else 1 sigma = 0.5 * (typeA.sigma + typeB.sigma) r6 = (sigma/r)**6 r12 = r6**2 u_lj = 4 * epsilon * (r12-r6) u_hps = lambda_ * np.array(u_lj) s = r<=sigma*2**(1/6) u_hps[s] = u_lj[s] + (1-lambda_) * epsilon u = u_elec + u_hps """ u = u_elec + u_wang return u
[docs] class MpipiBeads(PolymerBeads): def __init__(self, polymer, sequence=None, spring_constant = 19.19, rest_length = 3.8, **kwargs): if sequence is None: raise NotImplementedError # ... set random 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(_types[s], r, name = s, resid = i+1) def _join_adjacent_bead_groups(self, ids): ## Two consecutive groups if len(ids) == 2: b1,b2 = [self.children[i] for i in ids] """ units "10 kJ/N_A" kcal_mol """ bond = HarmonicBond(k = self.spring_constant, r0 = self.rest_length, range_ = (0,100), resolution = 0.01, max_force = 10) self.add_bond( i=b1, j=b2, bond = bond, exclude=True ) elif len(ids) == 3: ... else: pass
[docs] class MpipiModel(PolymerModel): def __init__(self, polymers, sequences = None, rest_length = 3.8, spring_constant = 19.19, debye_length = 7.95, damping_coefficient = 10, DEBUG=False, **kwargs): """ [debye_length]: angstroms [damping_coefficient]: ns """ logger.info("""You are using an implementation of the Mpipi Polymer model as described for proteins and RNA based on: Physics-driven coarse-grained model for biomolecular phase separation with near-quantitative accuracy. Published in final edited form as: Nat Comput Sci. 2021 Nov; 1(11): 732–743. Published online 2021 Nov 22. doi: 10.1038/s43588-021-00155-3 Please cite all appropriate articles!""") if 'timestep' not in kwargs: kwargs['timestep'] = 10e-6 if 'cutoff' not in kwargs: kwargs['cutoff'] = max(4*debye_length,20) if 'decomp_period' not in kwargs: kwargs['decomp_period'] = 1000 self.rest_length = rest_length self.spring_constant = spring_constant """ Assign sequences """ if sequences is None: raise NotImplementedError("MpipiModel 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 _types.items()] self.set_damping_coefficient( damping_coefficient ) """ Set up nonbonded interactions """ nonbonded = MpipiNonbonded(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 j in range(i,len(self.types)): t = self.types[j] self.add_nonbonded_interaction( interaction, typeA=type_, typeB=t ) def _generate_polymer_beads(self, polymer, sequence, polymer_index = None): return MpipiBeads(polymer, sequence, rest_length = self.rest_length, spring_constant = self.spring_constant, monomers_per_bead_group = self.monomers_per_bead_group, polymer_index = polymer_index )
[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__": pass """ from matplotlib import pyplot as plt nt = len(_types) # print("TYPES") # for n,t in _types.items(): # print("{}\t{}\t{}\t{}\t{}".format(t.name, t.mass, t.charge, t.sigma, t.lambda_)) type_string = 'WYFMLIVAPGCQNTSEDKHR' d = np.zeros([nt,nt]) for i in range(nt): n1 = type_string[i] t1 = _types[n1] for j in range(nt): n2 = type_string[j] t2 = _types[n2] d[nt-i-1,j] = "ERR"[(t1.name,t2.name)] plt.imshow(d.T) plt.show() """