Introducing FRET restraints

Introducing FRET restraints#

import numpy as np
from matplotlib import pyplot as plt
import pandas as pd
import mdtraj as md
import fretraj as ft

import importlib

Accessible contact volume#

Calculate ACVs for donor and acceptor dye on an solvated and energy minimized biomolecule.

Note: the generated .gro file must first be converted to a .pdb file using gmx trjconv

labels = ft.cloud.labeling_params('gmx/mn_riboswitch_labels.json', verbose=False)

struct = md.load_pdb('gmx/mn_riboswitch.pdb')
struct_nosolv = struct.remove_solvent()

don = ft.cloud.Volume(struct_nosolv, '1-B-U1-O5\'', labels)
acc = ft.cloud.Volume(struct_nosolv, '1-A-A1-O5\'', labels)

Compute a FRET trajectory from the ACVs

fret_traj = ft.cloud.FRET(don, acc, 'A1-U1', labels)

fret_traj.values
R0 (A) <R_DA> (A) <E_DA> <R_DA_E> (A) R_attach (A) R_mp (A)
value 54.0 96.1 0.04 91.9 86.6 94.4
std NaN 10.7 0.03 13.0 NaN NaN

Create a Plumed object from the loaded structure and the ACVs. Look for neighboring phosphorus atoms within a range of 15 A.

plumed = ft.restraints.Plumed(struct_nosolv, [don, acc], selection='name P', cutoff=15)

Set your desired FRET value that you want the simulation converge to. The corresponding \(\langle R_{DA,E}\rangle\) will be computed from the \(\langle E\rangle\). The \(\langle R_{mp}\rangle\) will be calculated from \(\langle E\rangle\) using a third order polynomial (Kalinin et al., Nat. Methods, 2012)

targetFRET = 0.6
mean_R_DA_E = ft.fret.mean_dist_DA_fromFRET(None,None,targetFRET,54)
targetRmp = ft.fret.R_DAE_to_Rmp(mean_R_DA_E)

Compile a plumed input file for the MD simulation

plumed.write_plumed('gmx/out/plumed_A1-U1.dat', targetRmp, 100, 100)

True

To illustrate the restaints save a visualization script for VMD or PyMOL which can be run with either vmd_vis -c filename.gro -x trajectory.xtc -v vis.vmd or pymol_vis -c filename.gro -x trajectory.xtc -v vis.py

plumed.write_vmd('gmx/out/vis.vmd')

plumed.write_pymol('gmx/out/vis.py')

Save force field and topology files files of the mean position pseudoatoms

plumed.write_pseudo('gmx/out/MP.pdb', 'gmx/out/MP.itp')

(1) Copy the "MP.itp" file into your force field directory

(2) Add the following lines to your "topology.top" file:
; Include topology for dye mean position
#include "<path/to/forcefield>/MP.itp"

#ifdef POSRES
[ position_restraints ]
1     1  1000  1000  1000
#endif

(3) Add the following line to your "atomtypes.atp" file in the force field directory:
MP        0.00000 ;

(4) Add the following lines to your "ffnonbonded.itp" file in the force field directory:
; Dummy mass for dye mean position
MP          35      79.90    0.0000  A   4.64693e-01  2.45414e-01

don.save_mp('gmx/out/mn_riboswitch_D.dat', units='nm', format='plain')
acc.save_mp('gmx/out/mn_riboswitch_A.dat', units='nm', format='plain')

First, backup the original .gro file. Then, insert the mean position coordinates given in MP.pdb at positions specified by riboswitch_D.dat and riboswitch_A.dat

mv em/mn_riboswitch.gro em/mn_riboswitch_original.gro
gmx insert-molecules -ci MP.pdb -f mn_riboswitch_original.gro -o mn_riboswitch.gro -ip out/mn_riboswitch_D.dat -replace water |& tee tmp_insert.dat
gmx insert-molecules -ci MP.pdb -f mn_riboswitch.gro -o mn_riboswitch.gro -ip out/mn_riboswitch_A.dat -replace water |& tee -a tmp_insert.dat

Run the following command will tell you how many solvent molecules have been replaced by the MP pseudoatoms and how to update the topology.

awk '$1 == "Replaced" {sum += $2}; END {print "\nIn your topology.top file under the section [ molecules ] do the following:\n(1) decrease the number of solvent molecules by " sum "\n(2) add the following line:\nMP        2"}' tmp_insert.dat && rm tmp_insert.dat
Contents