Collagen & Peptide NutritionNutrition and collagen guides

Nutrition guide

redesigned λ N peptide binds boxB RNA - PMC

. Author manuscript; available in PMC: 2017 Oct 15. Published in final edited form as: J Comput Chem. 2016 Aug 4;37(27):2423–2435. doi: 10.1002/jcc.24466 Abstract Our previously-developed peptide-design algorithm was improved by adding an energy minimization s

By Collagen & Peptide Nutrition Editorial TeamPublished

. Author manuscript; available in PMC: 2017 Oct 15.

Published in final edited form as: J Comput Chem. 2016 Aug 4;37(27):2423–2435. doi: 10.1002/jcc.24466

Abstract

Our previously-developed peptide-design algorithm was improved by adding an energy minimization strategy which allows the amino acid sidechains to move in a broad configuration space during sequence evolution. In this work, the new algorithm was used to generate a library of 21-mer peptides which could substitute for λ N peptide in binding to boxB RNA. Six potential peptides were obtained from the algorithm, all of which exhibited good binding capability with boxB RNA. Atomistic molecular dynamics simulations were then conducted to examine the ability of the λ N peptide and three best evolved peptides, viz. Pept01, Pept26 and Pept28, to bind to boxB RNA. Simulation results demonstrated that our evolved peptides are better at binding to boxB RNA than the λ N peptide. Sequence searches using the old (without energy minimization strategy) and new (with energy minimization strategy) algorithms confirm that the new algorithm is more effective at finding good RNA-binding peptides than the old algorithm.

Keywords: Peptide-design algorithm, λ N peptide-boxB RNA, RNA-binding peptide, RNA tethering, Energy minimization strategy

Graphical abstract

Peptide-design algorithm with a new addition of energy minimization strategy was used to search novel, high-affinity peptides in substitution of the native λ N peptide to bind to boxB RNA. Atomistic molecular dynamics simulations were then employed to examine the binding ability of the native and the evolved peptides to the RNA target. Simulation results suggested that the improved algorithm enables better design of peptides to tether boxB RNA.

Introduction

Proteins that bind to specific RNA motifs comprise useful tools for basic research and biotechnology. To date, the protein-RNA pairs in widest use are derived from bacteriophages, since bacteriophage coat proteins have evolved to bind their cognate RNA genomes with high affinity, and in general, the RNA motifs recognized by such coat proteins are distinct from those found in either eukaryotic or prokaryotic cells [1]. Such bacteriophage systems have been harnessed to tether proteins to RNA in vivo in order to study the function of RNA-binding proteins, dissect RNA-protein or RNA-RNA interactions, and visualize RNA localization, including single-molecule imaging [2]. Such bacteriophage systems have been also used in biotechnology applications ranging from drug delivery via virus-like particles [3] to building spatially-defined protein assemblies using RNA scaffolds [4].

Currently, the potential usage of bacteriophage protein-RNA tethering systems is limited by the properties of naturally-evolved “parts”. The most commonly used system is the MS2 bacteriophage coat protein, which is a dimeric protein that binds to its cognate MS2 RNA stem loop [5-7]. The MS2 coat protein dimer binds to the MS2 RNA stem loop with a Kd of (1-3)×10−9 M, and most applications employ a high affinity mutant of the MS2 RNA stem loop which binds the coat protein with a Kd of (2-6)×10−10 M [8, 9]. Most applications employ a tandem dimer of the MS2 coat protein, which is more stable than the natural monomer. An alternative system comprises the coat protein of the λ bacteriophage (λ N protein), which binds its cognate RNA stem loop (termed boxB RNA) with a Kd of 1.3×10−9 M [10-12]. An attractive feature of the λ N-boxB system is that only the first 22 amino acids of the λ N protein are required to mediate boxB RNA binding, and this λ N peptide may be useful for applications where the large tandem MS2 dimer is too bulky, e.g., disrupting the interaction under investigation. However, λ N peptide binds boxB RNA with an affinity that is too low for some applications. Thus, developing λ N peptide variants that bind boxB RNA with a higher affinity, and potentially with a range of affinities, would comprise unique and powerful tools for biotechnology and basic research. An attractive strategy for developing such variants is the use of computational protein/peptide design.

During the past decade, a number of researchers have strived to develop computational algorithms [13, 14] to de novo design or redesign peptides and proteins for particular purposes, e.g. strengthening the thermodynamic stability of enzymes [15, 16], modifying immune epitopes with structural specificity to antibodies [17, 18], etc. The side-chain grafting method, a molecular modeling technique, was used by Correia et al. [18] to redesign the scaffolds for the HIV epitope 4E10. They found that their designed epitope not only prevented undesired immune reactions but also exhibited a high binding affinity to monoclonal antibody (mAb) 4E10. Using a Monte Carlo (MC) approach, Ambroggio et al. [19] successfully designed a 32-mer amino acid sequence that can be stabilized simultaneously in two distinct protein conformations: a coiled coil and a 2C-2H zinc finger, depending upon the pH of solution. This study explored the underlying conformational switch mechanism in response to the changes of biological environment. Cochran et al. [20] employed a self-consistent mean field (SCMF) approach to reengineer natural protein scaffolds in a completely de novo design so that the redesigned proteins could bind selectively to a specific non-biological cofactor, e.g. DPP-Fe cofactor. In our previous work [21], a search algorithm combining MC and SCMF techniques was developed to evolve peptide sequences that have good binding capability to the anticodon stem and loop (ASL) of human lysine tRNA species, viz. ASLLys3, with the ultimate purpose of breaking the replication cycle of human immunodeficiency virus-1. This search algorithm has been validated experimentally by synthesizing and testing peptide sequences selected in this manner against ASLLys3 [22].

In this paper, we describe a major improvement to our peptide design algorithm and use it to generate a library of peptides predicted to bind boxB RNA with affinities greater than does the wild-type λ N(2-22) peptide. The major improvement is the addition of an energy minimization step to the search algorithm that can move the amino acid sidechains in a broad conformational space during peptide sequence evolution. Multiple searches are conducted on several conformations of the λ N(2-22) peptide-boxB RNA complex, and many potential peptides that have low score (high affinity) for binding with boxB RNA are found. Atomistic molecular dynamics (MD) simulations are employed to examine the binding ability of the original λ N(2-22) and our evolved peptides to boxB RNA. Finally, we compare the effectiveness of the old and new search algorithms and find that the addition of an energy minimization step greatly improves the performance of our new algorithm in the search for novel, high-affinity RNA-binding peptides.

Structures and Methods

Structures of λ N(2-22) peptide and boxB RNA

The tertiary structure of a 15-mer nutR boxB RNA hairpin complexed with the 36-mer N-terminal peptide of λ N protein is available in the Protein Data Bank (PDB ID: 1QFQ), which includes 29 possible binding conformations [23]. We extract a 21-mer peptide fragment, i.e. residues 2-22 from the λ N peptide, because these residues are in close contact with the 15-mer boxB RNA. Henceforth, this 21-mer peptide fragment will be referred to as λ N(2-22) peptide. Figure 1 shows the binding conformation of the complex formed by this peptide fragment with boxB RNA along with the sequences and secondary structures of the λ N(2-22) peptide and boxB RNA. All of the 29 complexes exhibit a bent α helix for the λ N(2-22) peptide which binds exclusively to the 5′ strand of boxB RNA [24]. According to convention, we will number the residues on the 21-mer λ N peptide fragment 2 to 22. The nucleotides on the 15-mer boxB RNA are also numbered for identification. We chose six of the 29 complexes randomly as candidates for conformational optimizations. These complexes are then used as initial conformations in the search algorithm to find other candidate peptides that would have high affinity to boxB RNA. Details regarding how we optimized the conformations of the six complexes using AMBER package will be given in section 3.1.

Figure 1.

Structure of the complex formed by λ N(2-22) peptide and boxB RNA hairpin based on PDB ID: 1QFQ showing the RNA (green ribbon), the peptide (orange ribbon), the positively charged amino acids (blue) and the negatively charged amino acids (red). The other nucleotides and amino acids are not shown for clarity. The secondary structure of the 15-mer nutR boxB RNA hairpin is shown in the box on the left. The primary sequence of λ N peptide from aspartate (D) at site 2 to asparagine (N) at site 22 is shown in the box on the top.

Search algorithm procedure

Figure 2 shows a flow sheet that illustrates the basic steps in the new algorithm. The conformation of boxB RNA and the backbone conformation of the peptide are fixed during sequence evolutions. Unlike our previous algorithm [21], the new algorithm does not need the self-consistent mean field technique to determine appropriate rotamer combinations from a library of fixed sidechain configurations. Instead we add an energy minimization step to determine optimal sidechain configuration for the amino acid repacking. This enables the sidechains to move more freely in configuration space. A detailed description of the energy minimization strategy will be given later in Figure 3. In the search algorithm, there are two types of trial “moves” to change the peptide sequence (Figure 2). The first type of trial “move” is a random substitution of a new residue for an old one. The new residue should be of the same residue type as the old one to maintain the peptide's hydration properties (see below). All possible rotamers on the new residue are then energy minimized, after which the optimal sidechain configuration with the lowest score (see below) is chosen to repack the side chain. The second type of trial “move” is a random exchange of two chosen residues, regardless of their residue types. All possible rotamers on the two exchanged residues are then energy minimized, after which the optimal sidechain configuration for each residue is chosen to repack the side chains.

Figure 2.

The flow sheet for the search algorithm.

Figure 3.

Block diagram showing the energy minimization strategy to determine optimal sidechain configuration for one amino acid repacking.

The procedure for the search algorithm is the following: a random initial sequence, S0, is generated and draped over the fixed backbone conformation; the score of this initial peptide-boxB RNA complex, Γscore0, is evaluated [25]. Next, a random number is generated to determine whether to mutate one random residue or to exchange two random residues. Once new residue(s) is(are) determined to use in substitution for old one(s), energy minimization is performed on its(their) sidechain(s) to get the optimal configuration(s), thereby resulting in the generation of a new peptide sequence, Si. The new peptide sequence is evaluated further by calculating its score, Γscorei. Finally, the new peptide sequence is accepted or rejected according to the Metropolis sampling method based on the evaluation of the score. After a total of 10,000 evolution steps, with each step containing either 21 mutation or 21 exchange attempts, the best peptides are identified from the sequence evolution. The code for the peptide-design algorithm was developed in our group and written in Fortran 90.

Setting constraints on the peptide's hydration properties in the search algorithm allows us to find candidate peptides with a range of affinities to boxB RNA. As the peptide's hydration property varies, its ability to bind to boxB RNA changes accordingly. We classified the 20 standard amino acids into six residue types according to their hydrophobicity, polarity, size and charge [26]. The first column in Table 1 gives the residue type and the second column lists the amino acids of that type. The 21-mer λ N(2-22) peptide that is of interest to us contains one hydrophobic residue (Nhydrophobic=1), three negatively charged residues (Nnegative charge=3), seven positively charged residues (Npositive charge=7), five hydrophilic residues (Nhydrophilic=5), five other residues (Nother=5) and no glycine (Nglycine=0). In this work, these numbers (Nhydrophobic, Nnegative charge, Npositive charge, Nhydrophilic, Nother, and Nglycine) will be kept the same during the process of sequence evolution. By doing so, we can maintain the hydration properties of all the evolved peptides the same as that in the λ N(2-22) peptide. In this way, we can avoid introducing excess hydrophilicity or positively charges in the evolved peptides due to strong negatively charges on boxB RNA.

Table 1.

The six residue types for the 20 standard amino acids.

Hydrophobic Leu, Val, Ile
Met
Phe
Tyr, Trp
Negatively charged Glu, Asp
Positively charged Arg, Lys
Hydrophilic Ser, Thr
Asn, Gln
His
Other Ala
Cys
Pro
Glycine Gly

The score function used in this search algorithm is:

Γscore=ΔGbinding+λ⋅(UVDWligand−bound+UELEligand−bound+GEGBligand−bound). (1)

The right hand side of the equation (1) contains two energy terms: one for the binding of the complex, as evaluated by equation (2), and the other for the stability of the ligand in the bound state within the complex, as evaluated by equation (3). The parameter λ that multiplies the second energy term is a weighting factor which is used to adjust the importance of the stability term of the ligand in the score function. The value of λ was set to 0.010 in this work; for discussion of how to set the value of λ see our previous work [25].

The binding free energy ΔGbinding for a ligand and a receptor is defined to be the difference between the free energy of the complex, and the free energies of the ligand and the receptor prior to binding [27]. It can be calculated according to:

ΔGbinding=GTOTcomplex−GTOTligand−GTOTreceptor, (2)

where GTOTcomplex, GTOTligand and GTOTrecepton represent the total free energies of the complex, the ligand and the receptor in solution, respectively. The total free energy GTOT (here, we neglect the entropy contribution in the search algorithm) of the molecule is calculated as follows:

GTOT=UINT+UVDW+UELE+GEGB+GGBSUR, (3)

where UINT, UVDW, UELE, GEGB and GGBSUR indicate the internal energy (INT), van der Waals energy (VDW), electrostatic energy (ELE), the polar solvation energy (EGB) and the non-polar solvation energy (GBSUR). The internal energy UINT is associated with the vibration of the bonds, bond angles, and torsional or dihedral angles. The van der Waals energy UVDW is modeled using the typical 12-6 Lennard-Jones equation. The electrostatic energy UELE follows Coulomb's law. The polar solvation energy GEGB is calculated based on the Generalized Born method with the Onufriev-Bashford-Case model [28]. The non-polar solvation energy GGBSUR is calculated by multiplying the solvent-accessible surface area of the solute molecules by the surface tension (in this work, the surface tension is set to 0.0072 kcal/mol/Å2). More details concerning the calculation of each energy term can be found in our previous work [22, 29]. It is worth noting that the internal energy ΔUINT is always zero in the calculation of the binding free energy since UINT is associated with the individual internal motions of ligand and receptor, and this term does not refer to the binding between ligand and receptor. We neglect the non-polar solvation energy ΔGGBSUR throughout the entire evolution process, because the calculation of ΔGGBSUR is time-consuming, and its value is quite small compared to the other energy terms, e.g. ΔUVDW, ΔUELE and ΔGEGB.

To accurately evaluate UELE and GEGB between the charged (or polar) molecules, we employ a variable internal dielectric constant model to capture electrostatic shielding effects. The expressions for the UELE and GEGB are:

UELE=12∑i∑j≠iqiqjεin(ij)rij, (4)
GEGB=−12∑i(1εin(ii)−1εsol)[qi2rii2+αi2exp(−rii24αi2)]−12∑i∑j≠i(1εin(ij)−1εsol)(qiqjrij2+αiαjexp(−rij24αiαj)), (5)

where qi and qj represent the charges on atoms i and j, respectively; rij is the distance between atoms i and j (thus, rii=0.00); εin(ij) is the internal dielectric constant whose value depends on the residue types for atoms i and j; εsol is the dielectric constant of the solvent (in this work, εsol=80.0 for water); and αi and αj are the effective Born radii of atoms i and j, respectively. In the residue-based dielectric model, we use distinct internal dielectric constants εin(i) for the polar and charged groups; these are listed in Table 2. For a pair of interacting atoms i and j, we always choose the larger of its two dielectric constants εin(i) and εin(j) as the residue-based internal dielectric constant εin(ij), viz. εin(ij)=max{εin(i), εin(j)}. Detailed descriptions regarding this point can be found in references [30, 31]. The values of the internal dielectric constants listed in Table 2 were obtained from reference [32]. All of the force field parameters in this work originated from the AMBER ff99SB force field library [33].

Table 2.

Internal dielectric constants (εin(i)) used in the variable dielectric model [32].

Residue type Internal Dielectric Constants (εin(i))
Lysine (K) 4
Glutamic Acid (E) 3
Aspartic Acid (D) 2
Arginine (R) 2
Histidine (H) 2
Phosphate linkage 4

An energy minimization strategy is used to determine optimal sidechain configuration for each amino acid that is to be repacked during the process of sequence evolution. We use the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method [34-37]. Details of the BFGS method for adjusting the configuration of an amino acid sidechain can be found in Supplemental Material. The energy minimization strategy consists of three stages, as shown in Figure 3. The first two stages eliminate bad sidechain conformers and select one potential sidechain conformer, and the last stage starts with this potential sidechain conformer to determine the optimal (final) sidechain configuration for repacking the amino acid. A detailed description of how the optimal sidechain configuration is obtained during the repacking of an amino acid sidechain is below.

We utilize Lovell's rotamer library [38] to obtain a series of initial rotamers for a particular amino acid. Consider an amino acid containing m initial rotamers, each of which is labeled by i (i.e. i=1, ..., m). All rotamers of the amino acid are subjected to sidechain configuration optimization. We begin by selecting one rotamer i. In the first stage of the strategy, the non-bonded VDW energy between this rotamer and the other amino acids on the peptide and nucleic acids on boxB RNA is calculated. The BFGS method is then used to generate a new “trial” sidechain conformer i(1) and the non-bonded VDW energy of this “trial” sidechain conformer is evaluated. If the value of the VDW energy is greater than 100.00 kcal/mol, the “trial” sidechain conformer i(1) and its corresponding rotamer i are discarded due to atom overlaps with other groups and another rotamer i+1 is then considered. If the value of the VDW energy is smaller than 100.00 kcal/mol, the “trial” sidechain conformer i(1) is eligible to continue on to the second stage. In the second stage of the strategy, the non-bonded VDW+ELE energies between this sidechain conformer i(1) and the other amino acids on the peptide and nucleic acids on boxB RNA are calculated. The BFGS method is used to further generate the second “trial” sidechain conformer i(2). As a result, a “temporary” peptide sequence is obtained by repacking the second “trial” sidechain conformer i(2) on the chosen amino acid. We next employ equation (1) to compute the score for the “temporary” peptide-boxB RNA complex, and record the relevant energy and configuration information. In order to obtain all possible sidechain conformers, we repeat the first two stages on all of the initial rotamers (viz. i=1, ..., m) for the chosen amino acid. This is feasible because the first two stages of the strategy are time efficient. However, since the last (third) stage of the strategy is very time-consuming, we select only one potential sidechain conformer obtained from the second stage to enter the last round for a refined sidechain configuration optimization. Herein, the sidechain conformer with the lowest score is considered as potential conformer, j(2). (Justification for this point is given in Supplemental Materials.) The last stage of the strategy starts with this potential sidechain conformer j(2) and then uses the BFGS method to determine a “final” sidechain conformer based on calculating the score of the peptide-boxB RNA complex. Eventually, the “final” sidechain conformer is used as the optimal sidechain configuration of the amino acid to repack the peptide backbone.

Procedure for atomistic molecular dynamics simulation

We carried out atomistic MD simulations using the AMBER 14 package to examine the dynamics of the binding process of the λ N(2-22) peptide and three evolved peptides (viz. Pept01, Pept26 and Pept28 that are obtained via the search algorithm) with boxB RNA. All simulations were run with explicit solvent in the canonical (NVT) ensemble. The AMBER ff99SB force field [33] was used to describe the nucleosides and amino acids. Each complex consisting of the peptide and boxB RNA was solvated in a periodic truncated octahedral box containing an 8 Ångstrom buffer of TIP3P water (~3000 TIP3P water molecules) surrounding the complex in each direction. Sodium counterions were added to neutralize the system. More simulation details can be found in our previous work [25, 29]. To reach system equilibrium, we ran an 120-ns simulation for the λ N(2-22) peptide-boxB RNA complex, and a 150-ns simulation for the three complexes formed by the evolved peptides and boxB RNA. The root mean square deviation (RMSD) from the initial conformation was calculated as a function of time to monitor whether the binding process reached a stable state or not. Clustering analysis was performed on the last 5 ns of the simulation trajectory to obtain a representative structure for the complex in solution.

Results

Generating initial conformations of the λ N(2-22) peptide-boxB RNA complex for input to search algorithm

The original PDB file (PDB ID: 1QFQ) for the λ N peptide-boxB RNA complex provides structural information on the 29 complexes [23]. We identify the 29 complexes as 1qfq1, 1qfq2, ... 1qfq29, where the “1qfq” indicates the PDB entry (1QFQ) and the last number reflects the index within the 29 complexes. Since hydrogen atoms are absent from the original PDB file, the xLeap program was used to automatically add the missing hydrogen atoms on each amino acid. The binding free energies of these 29 complexes were positive, which indicated that the addition of the missing hydrogen atoms caused some bad contacts between the λ N(2-22) peptide and boxB RNA. (The energies of the 29 λ N(2-22) peptide-boxB RNA complexes are given in supplement table 1.) In order to remove the bad contacts from the λ N(2-22) peptide-boxB RNA complex, the AMBER 14 package was run for a total of 500 energy optimization cycles: 250 steepest descent cycles followed by 250 conjugate gradient cycles. Six of the 29 original complexes were randomly chosen and subjected to an AMBER energy optimization: 1qfq1, 1qfq4, 1qfq15, 1qfq21, 1qfq26 and 1qfq28.

Table 3 lists the various contributions to the binding free energy, i.e. the VDW, ELE+EGB, and GBSUR energies, for the six optimized complexes. By comparing Table 3 and Supplemental Table 1, we see that after this simple energy optimization to remove atomic overlaps, the binding free energies of the six complexes dropped from high positive values to very low negative values, which indicates that the optimized λ N(2-22) peptides are in good contact with boxB RNA. The decrease in the binding free energy is mainly attributed to the change of the VDW and ELE+EGB energies. Since the VDW and ELE+EGB interactions are quite sensitive to the distance between atoms, elimination of the atomic overlaps helps avoid high VDW and ELE+EGB energies and strengthens the binding of the λ N(2-22) peptide and boxB RNA. It is worth noting that the GBSUR energy did not change a lot as a result of the energy optimization. For this reason, we decided to neglect the GBSUR contribution when evaluating the score function in the evolution process (more justification for this choice will be given later in the text). The scores of the six optimized λ N(2-22) peptide-boxB RNA complexes are shown in Table 3.

Table 3.

Contributions to the binding free energy for the six λ N(2-22) peptide-boxB RNA complexes after energy optimization using the AMBER 14 package.

VDW ELE+EGB GBSUR Binding Energy without GBSUR ΔGbinding (kcal/mol) Γ score
1qfq1 −56.69 11.08 −10.21 −45.62 −55.82 −57.18
1qfq4 −55.70 15.77 −9.98 −39.93 −49.91 −51.27
1qfq15 −59.72 14.74 −10.19 −44.98 −55.17 −56.52
1qfq21 −60.21 14.66 −9.80 −45.55 −55.36 −57.05
1qfq26 −66.31 9.96 −10.13 −56.35 −66.48 −67.87
1qfq28 −64.71 15.65 −10.66 −49.06 −59.72 −60.58

Sequence evolution starting with the backbone conformation of the six optimized complexes

The new search algorithm was performed on six different random sequences draped on the backbone conformations of the six optimized λ N(2-22) peptide-boxB RNA complexes. In the sequence evolution, the conformations of the complexes are kept fixed, but the peptide sequences are altered. Figure 4 shows the score profiles versus the number of search steps. After a quick drop in the scores at an early stage, all profiles fluctuate about certain levels, indicating that the amino acids with high energy contributions (negative binding free energy) quickly occupy significant sites along the chain, while occupancy of insignificant sites along the chain varies frequently as other amino acids sample these sites. A detailed discussion of this point will be given in the section on energy analysis. By examining the score profiles over the process of evolution, we identified the lowest score in each profile which corresponds to the best peptide sequence for each search. As a result, a total of six best peptide sequences were found; their backbone conformations are the same as the backbone scaffolds within the six optimized λ N(2-22) peptide-boxB RNA complexes.

Figure 4.

The profiles of score vs. number of evolution steps during the sequence evolution for the six λ N(2-22) peptide-boxB RNA complexes.

Table 4 lists the sequences of the λ N(2-22) peptide and the six best evolved peptides that exhibit a high affinity to boxB RNA. For simplicity, we named each sequence, viz. Pept01 (from 1qfq1), Pept04 (from 1qfq4), Pept15 (from 1qfq15), Pept21 (from 1qfq21), Pept26 (from 1qfq26) and Pept28 (from 1qfq28). A comparison between Table 3 and Table 4 clearly demonstrates that our evolved peptides have lower scores than the λ N(2-22) peptide when accessing boxB RNA, implying that these evolved peptides have better capability to bind to boxB RNA. For instance, for 1qfq28, the score for the complex between λ N(2-22) peptide and boxB RNA is −60.58 kcal/mol (Table 3), while the score for the complex formed by Pept28 and boxB RNA is −103.85 kcal/mol (Table 4). Thus, the evolved Pept28 has a score that is 71.43% higher than the λ N(2-22) peptide. (Detailed discussion of this point is given in Supplement materials.) All of the evolution results confirm that the new algorithm is effective at finding other good RNA-binding peptides.

Table 4.

The best peptides evolved based on the backbone scaffolds of the six λ N(2-22) peptide-boxB RNA complexes. The residues that are identical at the same sites in the evolved and λ N(2-22) peptides are highlighted in bold. The (random) starting sequences for each search are shown in the footnote. The score and binding energy without GBSUR of the six λ N(2-22) peptide-boxB RNA complexes can be found in table 3.

site 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 Binding Energy without GBSUR Score (kcal/mol)
λ N(2-22) peptide D A Q T R R R E R R A E K Q A Q W K A A N
Pept01 (1qfq1)a Y A R R R Q Q C R Q C D R Q E R Q R A E C −91.27 −99.57
Pept04 (1qfq4)b P A R N R R R T A N E A R Q D R N W C E R −83.69 −90.46
Pept15 (1qfq15)c E A R C R R N A R N Q D R R C W N H A D R −83.78 −94.02
Pept21 (1qfq21)d F A R R R R Q A R Q A C Q Q D R Q R E A E −82.66 −94.06
Pept26 (1qfq26)e Q A R R R R N Q A Q C A R R C R W Q E E E −91.36 −101.44
Pept28 (1qfq28)f D A R R R Q Q A A R C Q R Q C R Q R D D M −94.01 −103.85

a

starting sequence: ECHQKKARKQPDDARSARNWR, score: 199.38

b

starting sequence: HSQKKCKARDRKEAKDTCSWC, score: 18.84

c

starting sequence: RARAENCKDKPQRKSQKCTDY; score: 33.04

d

starting sequence: KADRCEKRTAAIKCDRQKNQS, score: −18.94

e

starting sequence: PCKRCEARANDRSSQKEYQKK; score: 35.79

f

starting sequence: KCRYAQDAESDCNKKNKRQAR, score: 247.70.

We also examined the similarity of the evolved peptides with the λ N(2-22) peptide. The residues that are identical at particular sites in the evolved and λ N(2-22) peptides are highlighted in bold (Table 4). Among the evolved peptides, Pept21 shares the most (seven) residues in common with the λ N(2-22) peptide with three alanines (A) at sites 3, 12 and 21, three arginines (R) at sites 6, 7, and 10, and one glutamine (Q) at site 15. Also, Pept28, which has the lowest score (highest affinity), has five residues identical to the λ N(2-22) peptide: one alanine (A) at site 3, two arginines (R) at sites 6 and 11, one aspartic acid (D) at site 2 and one glutamine (Q) at site 15.

Analysis of the energy for the three top evolved peptide-RNA complexes

To explore the roles of the peptide residues in binding to boxB RNA, we analyze the binding free energies of the three top evolved peptides, i.e. Pept01 (score: −99.57 kcal/mol), Pept26 (score: −101.44 kcal/mol) and Pept28 (score: −103.85 kcal/mol). Figures (5a, 6a and 7a) show snapshots of the three top peptide sequences draped on the backbone conformations of 1qfq1, 1qfq26 and 1qfq28, respectively. The side chain configurations of the key residues and nucleotides are exhibited in distinct colors. The contributions to the score, binding free energy as well as the VDW, ELE+EGB, and GBSUR energies of the three complexes are shown. Additionally, the individual contributions of each residue and each nucleotide within the three complexes to the inter-chain VDW, inter-chain ELE+EGB and GBSUR energies are shown in Figures (5-b, c), (6-b, c) and (7-b, c). We next analyze the structures of the three peptide-boxB RNA complexes, compare them with each other and with their corresponding λ N(2-22) peptide-boxB RNA complexes, and discuss three aspects of the energy contributions below.

Figure 5.

Pept01 binds to boxB RNA with high affinity. (a) The binding structure of Pept01 and boxB RNA obtained from search algorithm showing the peptide backbone (orange ribbon) and the ribose-phosphodiester backbone (green ribbon). Several key residues and nucleotides are specified in blue and red colors, respectively. The table on the right characterizes the contributions of the different binding modes for the Pept01-boxB RNA complex. (b) Individual contributions of each nucleotide to the inter-chain VDW, inter-chain ELE+EGB, and GBSUR. (c) Individual contributions of each amino acid to the inter-chain VDW, inter-chain ELE+EGB, and GBSUR.

Figure 6.

Pept26 binds to boxB RNA with high affinity. (a) The binding structure of Pept26 and boxB RNA obtained from search algorithm. The table on the right characterizes the contributions of the different binding modes for the Pept26-boxB RNA complex. (b) Individual contributions of each nucleotide to the inter-chain VDW, inter-chain ELE+EGB, and GBSUR. (c) Individual contributions of each amino acid to the inter-chain VDW, inter-chain ELE+EGB, and GBSUR.

Figure 7.

Pept28 binds to boxB RNA with high affinity. (a) The binding structure of Pept28 and boxB RNA obtained from search algorithm. The table on the right characterizes the contributions of the different binding modes for the Pept28-boxB RNA complex. (b) Individual contributions of each nucleotide to the inter-chain VDW, inter-chain ELE+EGB, and GBSUR. (c) Individual contributions of each amino acid to the inter-chain VDW, inter-chain ELE+EGB, and GBSUR.

Firstly, we see from Figures (5, 6 and 7) and Table 3 that our three top evolved peptides exhibit a higher binding capability to boxB RNA (viz. a lower binding free energy without GBSUR) and a lower score than the λ N(2-22) peptide due to a sizeable decrease in the ELE+EGB energy. The negative VDW energy of the evolved peptides-boxB RNA complexes also improves the peptides’ binding capability to a certain extent. There is hardly any difference in the GBSUR energy for the evolved peptide and that for the λ N(2-22) peptide, which indicates that the GBSUR energy has no impact on sequence evolution. Taking for example Pept28 (Figure 7a) and the λ N(2-22) peptide in the 1qfq28 configuration (Table 3), the ELE+EGB energy drops from 15.65 kcal/mol for λ N(2-22) peptide to −13.03 kcal/mol for Pept28, resulting in an almost 2-fold improvement in the binding affinity. The decrease in the VDW energy from −64.71 kcal/mol for λ N(2-22) peptide to −80.98 kcal/mol for Pept28 indicates that the VDW interaction also plays an important role in improving the peptides’ binding capability. Our neglect of the GBSUR contribution to the score function during sequence evolution is supported by the fact that the change in the GBSUR energy from the λ N(2-22) peptide (−10.66 kcal/mol) to Pept28 (−11.31 kcal/mol) is negligible.

Secondly, by analyzing the various contributions of each residue along the peptides to the binding free energy, we can determine which residues contribute high (negative) energies to the binding. In the three complexes between our evolved peptides and boxB RNA, the residues with high energy contributions (negative binding free energy) are particularly prone to be hydrophilic or positively charged. They strongly interact with the nucleotides on boxB RNA via both inter-chain ELE+EGB and inter-chain VDW interactions. The inter-chain ELE+EGB energy is mainly attributed to the hydrophilic or positively charged amino acids which attract the negatively charged phosphate linkages. The inter-chain VDW energy involves a wide distribution of pairs of residues and nucleotides. Take for example the data presented on Figure 5(b, c), tyrosine (Y) at site 2 and alanine (A) at site 3 on Pept01 have inter-chain VDW interactions with the base and sugar ring on the nucleotides. Nevertheless, the majority of the binding contribution involving Pept01 and boxB RNA is still attributed to seven positively charged arginines (R) at sites 4, 5, 6, 10, 14, 17 and 19, and five hydrophilic glutamines (Q) at sites 7, 8, 11, 15 and 18 (Figure 5c). These twelve significant sites on Pept01 also provide good steric spaces for residues R and Q to interact with the nucleotides of boxB RNA, thereby resulting in high energy contributions to the binding. Since the seven R residues and five Q residues occupy 12 out of a total of 21 sites on the peptide chain, it is relatively easy for them to wind up at significant sites when we exchange the amino acids in the search algorithm. This also explains why the score profiles in Figure 4 drop very quickly at an early stage of the search process. Besides the twelve significant sites (including seven arginine sites and five glutamine sites), other sites along the chain are not as significant. As the amino acids sample the insignificant sites during the sequence evolution, the occupancy of these sites varies quite frequently, causing a persistent medium-sized fluctuation in the score profiles (Figure 4).

Thirdly, analysis of the overall energy distribution in the complexes reveals that the majority of the binding energy comes from the first eight nucleotides on boxB RNA and from certain hydrophilic and positively charged residues on the peptides. Examination of the various energy contributions for the three best complexes shows that the inter-chain ELE+EGB energy is attributed to several individual nucleotides and residues; whereas the inter-chain VDW energy is distributed relatively uniformly along the sequences of boxB RNA and the peptides. Taking for example the Pept26-boxB RNA complex (Figure 6-b, c), most of the inter-chain ELE+EGB energies for boxB RNA are attributed to the nucleotides U5, G6, A7 and A8. However, the inter-chain VDW energy associated with these four nucleotides contributes only 54.2% of the total inter-chain VDW energy; 45.8% of the total inter-chain VDW energy comes from the other nucleotides (see Figure 6b). As for Pep26, a large portion of the inter-chain ELE+EGB energy is attributed to the hydrophilic asparagine (N) at site 8, glutamine (Q) at site 11 and the positively charged arginines (R) at sites 7, 14 and 15. However, these residues contribute only 40.4% of the inter-chain VDW energy. Other residues like alanine (A) at site 3 and tryptophan (W) at site 18 on Pept26 also contribute to the inter-chain VDW energy (see Figure 6c). According to the above energy analysis of Pept26 and boxB RNA, we concluded that the residues (N, Q and R) on the peptide chain contribute to both the inter-chain ELE+EGB and inter-chain VDW energies for the binding to boxB RNA (Figure 6c), and thus are responsible for both the binding affinity and specificity [29]. Other residues (A, R, Q and W) on the peptide chain have a strong preference for the bases of the nucleotides on boxB RNA via the inter-chain VDW interaction (Figure 6c), and thus contribute to the binding specificity as well [29]. A similar situation also holds for the Pept01-boxB RNA and Pept28-boxB RNA complexes.

Verification of the ability of the λ N(2-22) peptide and the evolved peptides to bind to boxB RNA by atomistic molecular dynamics simulations

Atomistic simulations were performed to examine and rank the binding ability of the designed peptides. Motivation for these calculations is the following. Low scores do not guarantee that the evolved peptides from the algorithm have a better ability to bind to boxB RNA than the λ N peptide; the search algorithm is merely a scheme for identifying promising peptide candidates with high affinity to boxB RNA. Atomistic MD simulations were performed to examine the binding abilities of the λ N(2-22) peptide, Pept01, Pept26 and Pept28 to boxB RNA. The starting structure of the λ N(2-22) peptide-boxB RNA complex was 1qfq28, and the starting structures of boxB RNA complexed with Pept01, Pept26 and Pept28 were the structures shown in Figures 5(a), 6(a) and 7(a), respectively. Note that the four complexes have the same peptide backbone conformation but different peptide sequences.

Figure 8 shows the RMSD profiles and the equilibrated binding conformations for the four complexes. The RMSD profiles in Figure 8 undergo only slight fluctuations within acceptable range at the end of simulations, implying that the binding conformations of the peptides and boxB RNA reached a stable equilibration state. By comparing the conformations of the λ N(2-22) peptide-boxB RNA complex before (Figure 3) and after (Figure 8a) the MD simulations, we can see that when the λ N(2-22) peptide binds to boxB RNA in solution, its C-terminal domain undergoes a conformational change from a helix (Figure 1) to a β-strand (Figure 8a). This is not surprising since short peptides in solution have more conformational flexibility than when they are part of a larger peptide, due to the absence of restrictions from surrounding protein residues. The λ N(2-22) peptide is a short peptide fragment in the λ N protein, thus it exhibits a larger conformational flexibility in binding to boxB RNA in our MD simulations than it would when it is part of the λ N protein. This is consistent with the experimental observations by Xia et al.[39] and Zhang et al.[40]. They found that the complex formed by the peptide fragment λ N(2-22) and boxB RNA exists in a dynamic equilibrium, with the N-terminal domain of the λ N(2-22) peptide forming a helix that stacks with the RNA, but the C-terminal domain adopting multiple discrete conformations within the complex. Similar behaviors also take place in our atomistic MD simulations for the other short binding peptides, Pept01, Pept26 and Pept28 as demonstrated in Figure 8. Thus, the simulated binding structures are not exactly the same as the structures used in the design algorithm.

Figure 8.

Atomistic MD simulation studies of the binding of (a) the λ N(2-22) peptide, (b) Pept01, (c) Pept26 and (d) Pept28 with boxB RNA. The profile of RMSD vs. time was shown on the left side of each sub-picture to assure that our simulations reach conformational equilibrium. The binding structure of the peptide and boxB RNA was shown on the right side of each sub-picture, which was extracted as a structural representative by performing cluster analysis for the last 5-ns trajectory. The RNA is represented by the green ribbon. The peptide is represented by the orange ribbon, with the positively charged amino acids in blue and the negatively charged amino acids in red.

The enthalpy (ΔHbindingavg) and entropy (TΔSbindingavg) for the λ N(2-22) peptide and the three evolved peptides bound to boxB RNA over the last 5 ns of the simulation trajectories were calculated using the Molecular Mechanics/Generalized Born Surface Area (MM/GBSA) method [27] with the residue-based dielectric constant model [32] and normal mode analysis approach [41], as listed in Table 5. Since ΔG = ΔH – TΔS, we computed the mean binding free energy (ΔGbindingavg) of the four peptide-boxB RNA complexes in their last 5-ns simulations. It is clear in Table 5 that the abilities of the three evolved peptides bound to boxB RNA, as measured by the values of (ΔGbindingavg), are equivalent (−15.85 kcal/mol for Pept01, −16.01 kcal/mol for Pept26 and −15.92 kcal/mol for Pept28), and are better than the λ N(2-22) peptide (−9.60 kcal/mol). The good result for Pept28 is, however, somewhat misleading as Pept28 is a good binder mainly because its entropy (−41.91 kcal/mol) is much higher than the λ N(2-22) peptide's (−52.09 kcal/mol), not because its enthalpy is lower. (Since our search algorithm essentially looks for peptides with the lowest energy, one would expect that the designed peptides would have lower enthalpy.) Pept01 and Pept26 are good binders because of their low enthalpy; their entropies are somewhat higher than the λ N(2-22) peptide's (see Table 5). Thus the simulation results suggest that the new search algorithm shows promise for design of RNA-binding peptides with higher affinity to boxB RNA than the λ N(2-22) peptide. The ultimate test is, of course, experimental verification.

Table 5.

The mean values of the binding free energy, enthalpy and entropy for the λ N(2-22) peptide and the three evolved peptides bound to boxB RNA in the last 5 ns of the atomistic explicit-solvent MD simulation trajectories. Unit: kcal/mol

ΔGbindingavg ΔHbindingavg TΔSbindingavg
λ N(2-22) peptide-boxB RNA −9.60±0.17a −61.69±0.14 −52.09±0.10
Pept01-boxB RNA −15.85±0.17 −62.21±0.15 −46.36±0.10
Pept26-boxB RNA −16.01±0.19 −61.93±0.17 −45.92±0.11
Pept28-boxB RNA −15.92±0.15 −57.83±0.13 −41.91±0.09

a

Calculation errors are estimated by the standard errors of the mean.

Comparison of the efficiencies of the old and the new search algorithms

In order to compare the merits of the energy minimization, we conducted sequence searches using the old and the new algorithms starting from the same random sequence draped on the 1qfq28 backbone conformation. The flow sheet for the old search algorithm is shown in Supplemental Figure 1. In that algorithm, the mutation of one amino acid or the exchange of two amino acids involves selection of fixed amino acid sidechain configurations from a rotamer library. However, in the new algorithm (Figure 2), the sidechain configuration of each trial amino acid is adjusted using an energy minimization strategy, thereby enhancing the ability of the configuration to fit into the available space.

Figure 9 shows a comparison of the profiles of the score versus the number of steps in the sequence evolutions using the old and the new search algorithms. After a quick drop in the first 50 steps, the score for the old algorithm remains roughly constant, leveling off at an average value of approximate −45.00 kcal/mol. In contrast, the score profile for the new algorithm drops quickly in the early stage but then enters a period of medium-sized fluctuations with an average value of approximate −80.00 kcal/mol. Through determining the lowest score for each profile, we identified the two best peptide sequences from the old and new algorithms, respectively. The best peptide sequence (IARARAHARRQNRRENRQDEA) obtained from the old algorithm has a score of −57.52 kcal/mol, which is a little higher than −60.58 kcal/mol for the λ N(2-22) peptide in the 1qfq28 conformation (Table 3). It suggests that this peptide would not be as effective as the λ N(2-22) peptide when approaching boxB RNA. However, the best peptide sequence Pept28 (DARRRQQAARCQRQCRQRDDM) obtained from the new algorithm has a score of −103.85 kcal/mol, which is much lower than −60.58 kcal/mol for the λ N(2-22) peptide, indicating that this peptide has a higher affinity to boxB RNA than the λ N(2-22) peptide. Thus the addition of an energy minimization step to the algorithm greatly improves the performance of the new algorithm in the search for good peptide candidates with low scores.

Figure 9.

Comparison of sequence evolution between the old and the new search algorithms. The profiles of the score vs. the number of evolution steps are plotted.

Conclusion

In this work, we described the development and implementation of an algorithm to search for RNA-binding peptides that serve as useful tools in biotechnology and basic research. The RNA-binding peptides considered here were designed to have improved binding affinity for the boxB RNA compared to the native λ N peptide. Such improved peptides could ultimately enable new approaches for characterizing and manipulating RNA localization through the use of high affinity, RNA-peptide interactions. To achieve this goal, we improved our search algorithm so that it could generate a library of RNA-binding peptides with the same or even higher affinity than the λ N peptide when bound to boxB RNA.

By incorporating an energy minimization step in the algorithm, we can search for peptide sequences in a broad sidechain configuration space. At the core of the energy minimization strategy is the Broyden-Fletcher-Goldfarb-Shanno (BFGS) optimization method which adjusts the amino acid sidechain configurations during the sequence evolution process. All the original rotamers of the amino acid are energy minimized to optimize their sidechain configurations. The energy minimization strategy consists of three stages: the first two stages eliminate bad sidechain conformers and select one potential sidechain conformer, and the last stage starts with this potential sidechain conformer to determine an optimal sidechain configuration for the amino acid repacking. When employing the BFGS method to determine a new sidechain configuration in this strategy, we adopt different energy functions: the non-bonded VDW energy between a single rotamer and all the other groups on the peptide and boxB RNA which is calculated in the first stage; the non-bonded VDW+ELE energy between the amino acid sidechain and all the other groups which is calculated in the second stage; and the score of the peptide-boxB RNA complex which is calculated in the last stage.

We performed the new algorithm to evolve potential candidate peptides based on the fixed backbone conformations of six λ N(2-22) peptide-boxB RNA complexes. This led to six best peptides which exhibit higher binding capability with boxB RNA than the λ N(2-22) peptide. Subsequent analysis of the energy of the complexes formed by the three top peptides, viz. Pept01, Pept26 and Pept28, with boxB RNA revealed that the first eight nucleotides on boxB RNA and the hydrophilic and positively charged residues on the peptides account for the majority of the contributions to the binding energy.

Atomistic MD simulations were carried out to study the ability of the λ N(2-22) peptide and the three top evolved peptides to bind to boxB RNA. The binding free energy of the four complexes was calculated over the last 5-ns of at least 120-ns simulation trajectory using the MM/GBSA method with the residue-based dielectric constant model. Simulation results revealed that the computed binding free energy of the evolved peptide-boxB RNA complexes are lower than that for the λ N(2-22) peptide-boxB RNA complex, which suggests that our new search algorithm enables design of good RNA-binding peptides to recognize and load boxB RNA.

To highlight the merits of the new energy-minimization based search algorithm, we compared the results of the old and the new search algorithms on the λ N(2-22) peptide-boxB RNA complex in the fixed backbone conformation, 1qfq28. The evolution results confirmed that the new algorithm is effective at finding other good RNA-binding peptides, because the peptides evolved from the new algorithm have a higher affinity to boxB RNA than those from the old algorithm. Adding an energy minimization step to the algorithm evidently improves its performance in the search for good RNA-binding peptides.

Supplementary Material

Supp Info

Acknowledgements

Financial support for this work was awarded to C.K.H. by the National Institutes of Health USA (EB006006). This work was also supported in part by the NSF's Research Triangle MRSEC, DMR-1121107. J.N.L would like to acknowledge support from a 3M Non-tenured Faculty Award and the Northwestern University Prostate Cancer Specialized Program of Research Excellence (SPORE) through NIH award P50 CA090386. This work was also supported to M.E.H. by the National Science Foundation's Graduate Research Fellowship Program (NSF GRFP) award DGE-0824162. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1053573. We thank Texas Advanced Computing Center (TACC) and San Diego Supercomputer Center (SDSC) for providing us computing time.

References

Associated Data

This section collects any data citations, data availability statements, or supplementary materials included in this article.

Supplementary Materials

Supp Info

Continue reading

Related guides

· Collagen & Peptide Nutrition Editorial Team

PEP-FOLD Peptide Structure Prediction Server

Welcome to the PEP-FOLD 2011 improved service! PEP-FOLD is a de novo approach aimed at predicting peptide structures from amino acid sequences. This method, based on structural al…

Read guide →

· Collagen & Peptide Nutrition Editorial Team

Peptide Library Design Tools

Enter Protein Sequence Amino Acid 3 letter code 1 letter code Alanine Ala A Cysteine Cys C Aspartic acid Asp D Glutamic acid Glu E Phenylalanine Phe F Glycine Gly G Histidine His …

Read guide →