This tutorial demonstrates how to use perturbation response scanning (PRS) to determine sensors and effectors, which are important for allosteric signal transduction. The PRS approach is derived from linear response theory where perturbation forces are applied via a covariance matrix, which can be derived from elastic network models or MD simulations.
The example used in this tutorial is the AMPA-type ionotropic glutamate receptor (AMPAR; PDB 3kg2), which we recently studied using this method (see Figure 6 of Dutta et al. 2015, Structure 23(9):1692-1704).
First we need to import required packages.
from prody import * import numpy as np from pylab import * %matplotlib inline confProDy(auto_show=False)
@> ProDy is configured: auto_show=False
First, we parse a structure that we want to analyse with PRS. For this tutorial we will fetch the near-intact full length AMPAR structure 3kg2 from the PDB. In the same step, we can create a selection containing the calpha atoms, which we will use in downstream steps.
ampar_ca = parsePDB('3kg2', compressed=False).calpha
@> Connecting wwPDB FTP server RCSB PDB (USA). @> 3kg2 downloaded (3kg2.pdb) @> PDB download via FTP completed (1 downloaded, 0 failed). @> 22686 atoms and 1 coordinate set(s) were parsed in 0.11s.
Next, create an ANM instance and calculate modes from which the covariance matrix can be calculated. By default, n_modes=20, giving us the first 20 modes. By setting n_atoms to None, we ask for all modes rather than a subset. We could alternatively apply the PRS to another model from which a covariance matrix could be derived such as PCA, GNM or an MD simulation.
anm_ampar = ANM('AMPAR 3kg2') anm_ampar.buildHessian(ampar_ca) anm_ampar.calcModes(n_modes=None) saveModel(anm_ampar,'3kg2',matrices=True) writeNMD('anm_3kg2.nmd',anm_ampar,ampar_ca)
@> Hessian was built in 8.12s. @> 9342 modes were calculated in 235.52s.
If you have previously run a block like the one above, the you can run one like that below to load the ANM back in again.
anm_ampar = loadModel('3kg2.anm.npz')
The PRS matrix is then calculated from the covariance matrix from the ANM of the AMPAR. The raw PRS matrix is symmetric and does not allow differentiation of sensors and effectors. We therefore normalize it by dividing each row by its diagonal element. This is achieved by specifying normMatrix=True.
prs_mat_ampar = calcPerturbResponse(anm_ampar, normMatrix=True)
@> Calculating covariance matrix @> Covariance matrix calculated in 42.8s. @> Perturbation response matrix calculated in 7.3s. @> Perturbation response scanning completed in 55.5s.
These profiles are the averages over the rows and columns of the normalized PRS matrix, respectively. We calculate them as follows:
eff_ampar, sen_ampar = \ calcPerturbResponseProfiles(prs_mat_ampar)
They can also be calculated during analysis steps as shown at the end.
Effectors are the most effective residues whose perturbation has large effects on the structure and dynamics. Conversely, sensors are the most sensitive residues, which respond most strongly to perturbations of effectors and themselves undergo structural changes. We can identify these two sets of residues by (a) viewing the profiles as graphs and deciding upon a reasonable cutoff or (b) coloring residues by effectiveness and sensitivity on the structure in a molecular graphics program. For the latter approach, we write the profiles into new PDB files in place of the b-factor field.
The functions associated with the two kinds of methods are able to calculate and return the effectiveness and sensitivity profiles if they aren’t provided but we pre-calculated them, which makes it easier to apply them for both analyses.
This approach is implemented in the showPerturbResponse function, which makes use of the matplotlib pyplot package. This function can calculate the PRS matrix (which will be normalized by default) as well as the effectiveness and sensitivity profiles if they are not provided. You then would need to provide a model from which the covariance matrix can be calculated instead.
showPerturbResponse(prs_matrix=prs_mat_ampar, effectiveness=eff_ampar, sensitivity=sen_ampar, atoms=ampar_ca, cmap=cm.Greys, norm=Normalize(0,np.max(prs_mat_ampar)/5.))
The last two options make the matrix greyscale and normalize the scale to make weaker signals more apparent. There are usually a few very strong signals, which otherwise drown out everything else. The current choice of capping at 1/5 of the max value looks reasonable for seeing the structure of the matrix.
This approach is implemented in the writePerturbResponsePDB function, which can also write out individual rows and columns of the PRS matrix to look at the perturbation and response properties of an individual residue as discussed later.
writePerturbResponsePDB(prs_mat_ampar, '3kg2.pdb', effectiveness=eff_ampar, sensitivity=sen_ampar)
@> 22686 atoms and 1 coordinate set(s) were parsed in 0.14s. @> The effectiveness and sensitivity profiles were written to 3kg2_effectiveness.pdb and 3kg2_sensitivity.pdb.
To visualize this data, load the files into the graphics program of your choice and color by b-factor. In VMD, you would do this through the Graphical Representation window (from Graphics > Representations menu). The window that comes up gives various Color Method options from which you would pick Beta. The residues with high b-factor are shown in blue followed by white and ending at red for low b-factor. You can change this via the Color Controls window (Graphics > Color); this has a Color Scale tab with a Method dropdown from which you can pick other options.
To look at the effectiveness that perturbing a residue has in eliciting a response in individual residues (instead of its overall effectiveness) or to look at the sensitivity of a residue to perturbations of individual residues (instead of its overall sensitivity), we read out rows or columns from the perturbation response matrix.
This can again be shown in a plot as in Figure 6D or on a structure as in Figure 7.
For this purpose, use the showPerturbResponseProfiles function, which creates line graphs (this can also be used for average effectiveness and sensitivity profiles).
showPerturbResponseProfiles(prs_matrix=prs_mat_ampar, atoms=ampar_ca, resnum=84, chain='B', overlay=True)
The overlay option can be removed to show the chains one after another. You can also use this for plotting overall effectiveness and sensitivity profiles, which can be provided or calculated. As with writePerturbResponsePDB and showPerturbResponse, you can also ask for the profiles to be returned for further analysis.
To this aim, we use writePerturbResponsePDB again but we provide a residue number and optionally a chain ID or a concatenated subset of chain IDs (e.g. 'AC'). A PDB will be produced for each chain requested (if no chain IDs are provided then all chains with that residue number will be assumed).
The default behaviour in this usage mode is to read out a row from the PRS matrix, which is a residue-specific effectiveness profile. For example, to generate Figure 7B, use the following command:
writePerturbResponsePDB(prs_mat_ampar, '3kg2.pdb', chain='A', resnum=475)
You can also ask for columns rather than rows to get residue-specific sensititivity profiles as in Figure 7A and C:
writePerturbResponsePDB(prs_mat_ampar, '3kg2.pdb', chain='B', resnum=84, direction='sensing')
writePerturbResponsePDB(prs_mat_ampar, '3kg2.pdb', chain='A', resnum=475, direction='sensing')
The first two arguments have default order so don’t need to be named: (1) the PRS matrix and (2) the input PDB file, which will be modified with data from the PRS matrix or profiles. Any other arguments are optional. As we saw above, if you don't enter a residue number, it writes out the effectiveness and sensitivity profiles.
Another example is that you can ask for the profiles to be returned for further analysis using the returnData option as below. This could also be used for calculating the effectors and sensors if you want to straight away write them to PDB before doing anything else as shown.
eff_ampar, sen_ampar = writePerturbResponsePDB(prs_mat_ampar, '3kg2.pdb', returnData=True)
A similar command form could be used with showPerturbResponse or showPerturbResponseProfiles.