This tutorial is thought to illustrate how it is possible to compute SAXS intensities from single PDB files or trajectories using PLUMED. In particular, we will show how to compute scattering intensities with the hybrid coarse-grained/atomistic approach that is described in [86] . The tutorial will provide basic instructions to prepare files for the back-calculation of SAXS intensities using the hybrid approach (this process is simplified by the possibility to use an ad-hoc python script). Further, it is explained how to adjust the plumed input file for specific purposes, e.g. to use SAXS data as restraints in MD simulations or to compare SAXS intensities computed with both the atomistic and the hybrid approach.
Once this tutorial is completed users will be able to:
The TARBALL for this project contains the following files:
This tutorial has been tested on PLUMED version 2.5.1
Calculations of Small-angle X-ray scattering (SAXS) intensities from a structure of N atoms could be extremely demanding from a computational perspective as it is an \(O(N^2)\) problem. This issue has limited the applicability of SAXS in numerous situations, including their use as restraints in combination with MD simulations. A strategy to reduce the cost of computing SAXS from atomic structure consists in using a coarse grain representation of the structure, represented as a collection of M beads with M<N, each comprising a variable number of atoms. Niebling et al [84] have previously derived the Martini beads form factors for proteins, showing how this approach can be almost 50 times faster than the standard SAXS calculation. More recently, Martini beads form factors for nucleic acids have been also derived and, together with the ones for proteins, have been implemented in PLUMED. Importantly, it has been shown that the computation of scattering intensities using these Martini form factors achieves a good accuracy, as compared to the atomistic ones, for đť‘ž values up to 0.45 Ă…-1.
The Martini form factors can be exploited within a hybrid multi-resolution strategy to speed up SAXS profiles calculations, both for the back-calculation of scattering intensities from atomistic PDB or trajectory files and for the inclusion of SAXS data as restraints in experimental-driven all-atom simulations (e.g. METAINFERENCE). This can be achieved using PLUMED, which computes on the fly the virtual position of the Martini beads from the atomistic 3D-structure and then uses the centres of mass of these beads, along with Martini form factors, to quickly back-calculate SAXS intensities.
Given a PDB file, PLUMED is able to compute SAXS intensities for the molecule in the PDB and to compare these intensities with the experimental ones stored in a data file. It is possible to apply the same procedure to all the conformations visited during a MD trajectory. This is achieved by using the PLUMED driver utility and the SAXS variable of the ISDB module. While computing scattering intensities with a full atomistic representation is quite easy (see the SAXS keyword and later in this tutorial), in order to adopt the hybrid coarse-grained/atomistic approach a more elaborated procedure is needed, requiring the generation of few specific files to be used as input by driver. To facilitate this step, it is possible to use this script (martiniFormFactor_p3.py for python 3) and type in the bash shell:
python martiniFormFactor_p3.py -f start.pdb -dat SASDAB7.dat [-unit Ang/nm -nq 15]
The input files are:
The files generated in this step, to be used later as input for driver, are the following:
The default files generated should be sufficient to post-process single PDB files or trajectories, by typing:
plumed driver --plumed plumed.dat --mf_pdb start.pdb
or
plumed driver --plumed plumed.dat --mf_xtc samplextc.xtc
The output files generated by this step are:
Let's see now the meaning of each line in plumed.dat and how it is possible to modify it for other purposes, e.g. to compare atomistic/coarse-grained scattering intensities and to perform Metainference simulations where SAXS data are used as restraints. Here is a sample of plumed.dat produced by the python script martiniFormFactor_p3.py:
MOLINFOSTRUCTURE=aacg_template.pdb # BEADS DEFINITION INCLUDEcompulsory keyword a file in pdb format containing a reference structure.FILE=plumed_beads.dat # SAXS saxsdata: SAXS ...compulsory keyword file to be includedATOMS=martiniThe atoms to be included in the calculation, e.g.MARTINI# You can use SCALEINT keyword to set appropriate scaling factor. # SCALEINT is expected to correspond to the intensity in q=0 # SCALEINT=( default=off ) calculate SAXS for a Martini modelQVALUE1=0.0111721000Selected scattering lengths in Angstrom are given as QVALUE1, QVALUE2, ...EXPINT1=0.0527250000Add an experimental value for each q value.QVALUE2=0.0368675000Selected scattering lengths in Angstrom are given as QVALUE1, QVALUE2, ...EXPINT2=0.0327126000Add an experimental value for each q value.QVALUE3=0.0625615000Selected scattering lengths in Angstrom are given as QVALUE1, QVALUE2, ...EXPINT3=0.0128316000Add an experimental value for each q value.QVALUE4=0.0882529000Selected scattering lengths in Angstrom are given as QVALUE1, QVALUE2, ...EXPINT4=0.0045545200Add an experimental value for each q value.QVALUE5=0.1139410000Selected scattering lengths in Angstrom are given as QVALUE1, QVALUE2, ...EXPINT5=0.0022799600Add an experimental value for each q value.QVALUE6=0.1396240000Selected scattering lengths in Angstrom are given as QVALUE1, QVALUE2, ...EXPINT6=0.0013048600Add an experimental value for each q value.QVALUE7=0.1653020000Selected scattering lengths in Angstrom are given as QVALUE1, QVALUE2, ...EXPINT7=0.0007215740Add an experimental value for each q value.QVALUE8=0.1909730000Selected scattering lengths in Angstrom are given as QVALUE1, QVALUE2, ...EXPINT8=0.0004340930Add an experimental value for each q value.QVALUE9=0.2166360000Selected scattering lengths in Angstrom are given as QVALUE1, QVALUE2, ...EXPINT9=0.0002717150Add an experimental value for each q value.QVALUE10=0.2422910000Selected scattering lengths in Angstrom are given as QVALUE1, QVALUE2, ...EXPINT10=0.0002574160Add an experimental value for each q value.QVALUE11=0.2679360000Selected scattering lengths in Angstrom are given as QVALUE1, QVALUE2, ...EXPINT11=0.0001878030Add an experimental value for each q value.QVALUE12=0.2935700000Selected scattering lengths in Angstrom are given as QVALUE1, QVALUE2, ...EXPINT12=0.0001592670Add an experimental value for each q value.QVALUE13=0.3191930000Selected scattering lengths in Angstrom are given as QVALUE1, QVALUE2, ...EXPINT13=0.0000811279Add an experimental value for each q value.QVALUE14=0.3448030000Selected scattering lengths in Angstrom are given as QVALUE1, QVALUE2, ...EXPINT14=0.0001110630Add an experimental value for each q value.QVALUE15=0.3703990000Selected scattering lengths in Angstrom are given as QVALUE1, QVALUE2, ...EXPINT15=0.0001264680 # METAINFERENCE # Uncomment the following keywords and adjust parameters to activate METAINFERENCE # DOSCORE NOENSEMBLE SIGMA_MEAN0=0 # REGRES_ZERO=500 # SIGMA0=5 SIGMA_MIN=0.001 SIGMA_MAX=5.00 # NOISETYPE=MGAUSS ... # METAINFERENCE # Uncomment the following keyword to activate METAINF # saxsbias: BIASVALUE ARG=(saxsdata\.score) STRIDE=10 # STATISTICS statcg: STATSAdd an experimental value for each q value.ARG=(saxsdata\.q_.*)the input for this action is the scalar output from one or more other actions.PARARG=(saxsdata\.exp_.*) # PRINT # Uncomment the following line to print METAINFERENCE output # PRINT ARG=(saxsdata\.score),(saxsdata\.biasDer),(saxsdata\.weight),(saxsdata\.scale), (saxsdata\.offset),(saxsdata\.acceptSigma),(saxsdata\.sigma.*) STRIDE=500 FILE=BAYES.SAXS # change stride if you are using METAINFERENCE PRINTthe input for this action is the scalar output from one or more other actions without derivatives.ARG=(saxsdata\.q_.*)the input for this action is the scalar output from one or more other actions.STRIDE=1compulsory keyword ( default=1 ) the frequency with which the quantities of interest should be outputFILE=SAXSINT PRINTthe name of the file on which to output these quantitiesARG=statcg.corrthe input for this action is the scalar output from one or more other actions.STRIDE=1compulsory keyword ( default=1 ) the frequency with which the quantities of interest should be outputFILE=ST.SAXSCGthe name of the file on which to output these quantities
The first two lines tell PLUMED to use aacg_template.pdb as a template and to read the file plumed_beads.dat (which rules how to compute the CENTER of the beads and defines the “martini” group, containing all the beads). Note that you do not need to prepare these two files, since they both are produced by the python script martiniFormFactor_p[2,3].py. These lines are followed by the SAXS keyword, labelled “saxsdata”. Here are defined the atoms to be used for SAXS calculations (in our case these are all the beads within the group “martini”) and the structure factors to adopt: in this case we are using the Martini form factors (flag MARTINI), alternatively it is possible to use the atomistic ones (flag ATOMISTIC) or define custom form factors using a polynomial expansion to any order (flag PARAMETERS). By default, PLUMED computes a scaling factor to fit experimental and back-calculated intensities, however it is possible to set manually this scaling factor using the flag SCALEINT, which is expected to correspond to the intensity in q=0. The following lines indicate the q-values at which the calculation of scattering intensities is required (QVALUE) and the corresponding intensities (EXPINT). These are the ones selected by default by the python script, but it is possible to manually adjust them and/or to add new values. Lastly, within the SAXS keyword, there are few lines that are needed to activate METAINFERENCE. These are commented since they are not necessary for the back-calculation of SAXS intensities, however you have to uncomment and adjust them if you want to perform Metainference simulations in which SAXS data are used as restraints. The same is true for the line containing the keyword BIASVALUE. In the last part of the file, we ask PLUMED to compute some statistics comparing experimental and back-calculated intensities; finally, we print the computed intensities and statistics into the SAXSINT and ST.SAXSCG files, respectively.
It is easy to modify the plumed.dat for your own purposes, for instance:
MOLINFOSTRUCTURE=aacg_template.pdb # BEADS DEFINITION INCLUDEcompulsory keyword a file in pdb format containing a reference structure.FILE=plumed_beads.dat # SAXS saxscg: SAXS ...compulsory keyword file to be includedATOMS=martiniThe atoms to be included in the calculation, e.g.MARTINI( default=off ) calculate SAXS for a Martini modelQVALUE1=0.010Selected scattering lengths in Angstrom are given as QVALUE1, QVALUE2, ...QVALUE2=0.020 # etc.. ... saxsaa: SAXS ...Selected scattering lengths in Angstrom are given as QVALUE1, QVALUE2, ...ATOMS=1-11104The atoms to be included in the calculation, e.g.ATOMISTIC( default=off ) calculate SAXS for an atomistic modelQVALUE1=0.010Selected scattering lengths in Angstrom are given as QVALUE1, QVALUE2, ...QVALUE2=0.020 # etc.. ... #PRINT PRINTSelected scattering lengths in Angstrom are given as QVALUE1, QVALUE2, ...ARG=(saxscg\.q_.*)the input for this action is the scalar output from one or more other actions.FILE=QVAL_CG PRINTthe name of the file on which to output these quantitiesARG=(saxsaa\.q_.*)the input for this action is the scalar output from one or more other actions.FILE=QVAL_AAthe name of the file on which to output these quantities