Recent site activity

APBS‎ > ‎Examples‎ > ‎Binding energies‎ > ‎

Ionic strength dependence of peptide-RNA interactions

This is example is taken directly from
García-García C, Draper DE. Electrostatic interactions in a peptide-RNA complex. Journal of Molecular Biology. 331 (1), 75-88, 2003. http://dx.doi.org/10.1016/S0022-2836(03)00615-6
with the minimized PDB files kindly provided by David Draper. Files associated with this example are linked below and provided with APBS in the examples/protein-rna directory.

Introduction

The example explores the electrostatic contributions to the binding interaction between a 22-residue α-helical peptide of protein λ with the "box B" RNA hairpin structure. In particular, this example uses nonlinear Poisson-Boltzmann equation calculations to look at the non-specific screening effects of monovalent salt on the peptide-RNA complex.
In the experimental portion of the García-García and Draper manuscript, they isolate the contribution of KCl concentration to the binding of the folded peptide with the folded RNA hairpin and determine a fairly linear relationship between the binding free energy \Delta_{\rm bind} G and the logarithm of the KCl concentration which yields

 
This slope can be used to determine the number  of KCl ions linked to the binding equilibrium through the expression
where RT is the thermal energy, to determine n = -4.4 \pm 0.2 for the RNA-peptide binding equilibrium.
García-García and Draper also used nonlinear Poisson-Boltzmann equation calculations to estimate the electrostatic contributions to the binding free energy as a function of the monovalent salt concentration. As discussed elsewhere, the Poisson-Boltzmann equation is only able to describe non-specific interactions of ions with solutes, including the effects of ion size and charge but otherwise ignoring the important differences between ionic species. Interestingly (and perhaps surprisingly), they find excellent agreement between the experimental binding energy dependence on [KCl] and their Poisson-Boltzmann calculations with equivalent concentrations of monovalent ions. This agreement strongly suggests that the binding of RNA and the peptide is primarily determined by electrostatic interactions. It also suggests that the primary interaction of the KCl with this system is through non-specific screening interactions. The García-García and Draper nonlinear Poisson-Boltzmann equation calculations gave


and n = -4.3 \pm 0.2 for KCl linkage to the RNA-peptide binding equilibrium.

The binding energy calculation

Structures

As mentioned above, a minimized structure for the peptide-RNA complex was provided by David Draper. The PQR file for this complex is available for download here and shown in the following figure:


The complex between a 22 residue α-helical peptide from the N protein of phage λ and an RNA hairpin.

PQR files for the isolated RNA and peptide components are also available or can simply be extracted from the PQR file for the complex.

The calculation

The basic aspects of the electrostatics calculation will be highlighted by calculating the electrostatic contribution to the binding energy at a 0.225 M monovalent salt concentration. In what follows, we'll assume you've named the peptide-RNA complex PQR file model_outNB.pqr, the RNA PQR file model_outBoxB19.pqr, and the peptide PQR file model_outNpep.pqr.
After reading in our PQR files with the READ command, our calculation will have three parts:
  1. Calculation of the total electrostatic energy (including self-interaction energies) of the peptide-RNA complex. This calculation is named complex in the input file.
  2. Calculation of the total electrostatic energy (including self-interaction energies) of the peptide. This calculation is named peptide in the input file.
  3. Calculation of the total electrostatic energy (including self-interaction energies) of the RNA. This calculation is named rna in the input file.
Each of these parts will be contained in a separate ELEC section of the input file which specifies the settings for our nonlinear Poisson-Boltzmann equation calculation. The calculations themselves will not be overly demanding, since we will use relatively coarse grids. This grid coarseness has a significant impact on the absolute electrostatic binding energy we obtain from this particular calculation: the calculated energy isn't converged with respect to grid spacing. However, the overall slope of binding energy with respect to monovalent ion concentration is rather insensitive with respect to the grid spacing, allowing us to save computational time and effort during the calculations. Finally, the calculation will conclude with a PRINT command which will combine the total energies from the three parts to obtain our approximate absolute electrostatic binding energy for the complex at 0.225 M monovalent salt concentration. It is very important to note that this absolute energy no meaning in isolation for several reasons:
  • It is not converged with respect to grid spacing
  • It does not contain other very important non-electrostatic aspects of the binding energy which are important for the measured affinity
For purposes of comparison, APBS 1.0.0 yields an energy of 113.81 kJ/mol for this particular calculation using the following input file:
   read
       mol pqr model_outNB.pqr
       mol pqr model_outNpep.pqr
       mol pqr model_outBoxB19.pqr
   end
   elec name complex
       mg-auto
       dime 65 97 129
       cglen 45.3322 54.9498 82.2633
       fglen 45.3322 52.3234 68.3902
       cgcent mol 1
       fgcent mol 1
       mol 1
       npbe
       bcfl sdh
       pdie 4.0
       ion charge 1 conc 0.225 radius 2.0
       ion charge -1 conc 0.225 radius 2.0
       sdie 80.0
       srfm mol
       chgm spl2
       sdens 10.00
       srad 1.40
       swin 0.30
       temp 298.15
       calcenergy total
       calcforce no
   end
   elec name peptide
       mg-auto
       dime 65 97 129
       cglen 45.3322 54.9498 82.2633
       fglen 45.3322 52.3234 68.3902
       cgcent mol 1
       fgcent mol 1
       mol 2
       npbe
       bcfl sdh
       pdie 4.0
       sdie 80.0
       ion charge 1 conc 0.225 radius 2.0
       ion charge -1 conc 0.225 radius 2.0
       srfm mol
       chgm spl2
       sdens 10.00
       srad 1.40
       swin 0.30
       temp 298.15
       calcenergy total
       calcforce no
   end
   elec name rna
       mg-auto
       dime 65 97 129
       cglen 45.3322 54.9498 82.2633
       fglen 45.3322 52.3234 68.3902
       cgcent mol 1
       fgcent mol 1
       mol 3
       npbe
       bcfl sdh
       pdie 4.0
       sdie 80.0
       ion charge 1 conc 0.225 radius 2.0
       ion charge -1 conc 0.225 radius 2.0
       srfm mol
       chgm spl2
       sdens 10.00
       srad 1.40
       swin 0.30
       temp 298.15
       calcenergy total
       calcforce no
   end
   print energy complex - peptide - rna end
   quit

Visualizing the results

This is also a very good system for visualization. You can visualize the electrostatic properties of the system using the methods described in the VMD tutorial and the PyMOL tutorial. Note that the input file above did not include write statements to write out electrostatic potentials or other properties. Therefore, you will either need to add this to the input file directly or calculate the electrostatic potentials through separate calculations within VMD or PyMOL. As an example, the following image illustrates some aspects of the electrostatic driving force for peptide binding to the negatively-charged RNA molecule:


The surface potential of the peptide helix binding to the negatively-charged RNA molecule. Surface potentials are shown from −5 kT/e (red) to +5 kT/e (blue)

As discussed above, the primary purpose of these calculations is to examine the dependence of binding energy on ionic strength. To gain some insight into this dependence, we can visualize the change in ionic environment around the RNA and peptide upon binding. First, we need to modify the APBS input file for this calculation to include a write qdens command in each of the ELEC statements. If we want to examine the change in ion charge density, then we need to add statements of the form
 write qdens dx qdens-...
to each ELEC input file section. If we want to examine the change in ion number density, then we need to add statements of the form
   write ndens dx ndens-...
to each ELEC input file section. In both cases, the
...
should be replaced with some string denoting the particular ELEC statement; e.g., complex, peptide, or rna. An example APBS input file is available here.

After running this modified calculation to output ion densities, you will have OpenDX maps of ion densities for the complex, RNA, and peptide components of the system. These can be visualized individually, but it's often more interesting to see how the densities change upon binding. APBS provides a tool tools/mesh/dxmath to subtract OpenDX maps of scalar quantities calculated on the same grids (e.g., same grid dimensions, spacings, etc.). The dxmath program is invoked as
   dxmath dxmath.in
where dxmath.in is an input file with the stack-based arithmetic operations to be performed on the OpenDX maps. Running dxmath with no arguments provides more information about its usage. For example, the input file
   qdens-complex-0.225.dx
   qdens-pep-0.225.dx -
   qdens-rna-0.225.dx -
   qdens-diff-0.225.dx = 
will create the OpenDX map qdens-diff-0.225.dx by subtracting the maps qdens-pep-0.225.dx and qdens-rna-0.225.dx from qdens-complex-0.225.dx. A similar input file can be used to generate a difference map for number densities.

The resulting number or charge density or density-difference OpenDX maps can be read into VMD through the File → New Molecule... menu. They are generally best visualized through isocontour or "field line" representations which are described (in the context of electrostatic potentials) in the VMD section. As an example, here is the charge density profile for the peptide:


Mobile ion charge density isocontours around a peptide alpha helix. Orange isocontours are positive charge density and green contours are negative charge density.

and the complex:


Mobile ion charge density isocontours around the peptide-RNA complex. Orange isocontours are positive charge density and green contours are negative charge density.

(the RNA charge density isocontours are rather boring). This picture shows isocontours of the change in charge density upon complex formation as calculated from the dxmath-generated difference map:



Isocontours of the change in mobile ion charge density around the peptide-RNA complex. Blue isocontours denote positive changes in charge density and red isocontours denote negative changes.

Finally, this picture shows isocontours of the change in number density upon complex formation as calculated from the dxmath-generated difference map:


Isocontours of the change in mobile ion number density around the peptide-RNA complex. Blue isocontours denote positive changes in number density and red isocontours denote negative changes.

The predominance of the negative isocontour suggests that a negative overall change in ion number density upon binding. As discussed above, experimental data suggests that this qualitative observation is correct. The next section will demonstrate how to generate a quantitative estimate of the linkage between monovalent ion density and the protein-RNA binding.

Examining ionic strength dependence

After understanding the binding energy calculation above, it should be straightforward to generalize this to perform binding energy calculations for a series of monovalent ion concentrations. García-García and Draper examine binding affinities over a range of ion concentrations starting from close to 0 mM and up to 800 mM. Although Poisson-Boltzmann theory is most appropriate for ionic strengths under 500 mM, we will cover the same range in our calculations. In fact, if you perform this same analysis using only concentrations under 500 mM, you obtain a very slightly different answer (try it!).

It would be exceedingly tedious to run the calculation described above for several ionic strengths, so I recommend creating an input file template which replaces all of the specific ionic strength values with a placeholder; e.g.,
 
  read
       mol pqr model_outNB.pqr
       mol pqr model_outNpep.pqr
       mol pqr model_outBoxB19.pqr
   end
   elec name complex
       mg-auto
       dime 65 97 129
       cglen 45.3322 54.9498 82.2633
       fglen 45.3322 52.3234 68.3902
       cgcent mol 1
       fgcent mol 1
       mol 1
       npbe
       bcfl sdh
       pdie 4.0
       ion charge 1 conc IONSTR radius 2.0
       ion charge -1 conc IONSTR radius 2.0
       sdie 80.0
       srfm mol
       chgm spl2
       sdens 10.00
       srad 1.40
       swin 0.30
       temp 298.15
       calcenergy total
       calcforce no
   end
   elec name peptide
       mg-auto
       dime 65 97 129
       cglen 45.3322 54.9498 82.2633
       fglen 45.3322 52.3234 68.3902
       cgcent mol 1
       fgcent mol 1
       mol 2
       npbe
       bcfl sdh
       pdie 4.0
       sdie 80.0
       ion charge 1 conc IONSTR radius 2.0
       ion charge -1 conc IONSTR radius 2.0
       srfm mol
       chgm spl2
       sdens 10.00
       srad 1.40
       swin 0.30
       temp 298.15
       calcenergy total
       calcforce no
   end
   elec name rna
       mg-auto
       dime 65 97 129
       cglen 45.3322 54.9498 82.2633
       fglen 45.3322 52.3234 68.3902
       cgcent mol 1
       fgcent mol 1
       mol 3
       npbe
       bcfl sdh
       pdie 4.0
       sdie 80.0
       ion charge 1 conc IONSTR radius 2.0
       ion charge -1 conc IONSTR radius 2.0
       srfm mol
       chgm spl2
       sdens 10.00
       srad 1.40
       swin 0.30
       temp 298.15
       calcenergy total
       calcforce no
   end
   print energy complex - peptide - rna end
   quit
You can then generate and use a series of input files simply through replacement of the IONSTR placeholder. For example, I ran my series of calculations using the following bash script:
for conc in 0.025 0.050 0.075 0.100 0.125 0.150 0.175 0.225 0.250 0.275 0.300 0.325 0.400 0.500 0.600 0.700 0.800; do
    infile=apbs-${conc}.in
    outfile=${infile%.in}.out
    tempfile=template.in
    cat ${tempfile} | sed -e "s/IONSTR/${conc}/g" > ${infile}
    echo "Starting calculation for ${conc} M ionic strength..."
    apbs ${infile} > ${outfile}
    echo "Done."
done
This will produce a series of APBS output files containing the calculated binding energies at the desired ionic strengths. Note that these energies are in kJ/mol and will need to be appropriately scaled for comparison with the data from García-García and Draper. The raw results from this calculation using APBS 1.0.0 are shown below:

Further analysis requires the ability to perform linear regression. There are many programs that do this; I have provided a simple Python script for this purpose here. This script is used as:
   cat file.dat | python fit.py
where file.dat is a file with whitespace-separated data in two columns. With linear regression, we can obtain



from



Likewise, we can also obtain


from

Lnconc-energykcal.png/400px-Lnconc-energykcal.png

Comparison with the values from the García-García and Draper paper (see above) shows good agreement between the calculated number of monovalent ions linked to the binding and the experimental number of KCl ions linked to the equilibrium.