Predicting the Initial Steps of Salt-Stable Cowpea ... - ACS Publications

Subscriber access provided by Stockholm University Library

Article

Predicting the Initial Steps of Salt Stable CCMV Capsid Assembly with Atomistic Force Fields Zoltan Antal, Janos Szoverfi, and Szilard N. Fejer J. Chem. Inf. Model., Just Accepted Manuscript • DOI: 10.1021/acs.jcim.7b00078 • Publication Date (Web): 06 Apr 2017 Downloaded from http://pubs.acs.org on April 6, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Information and Modeling is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 24

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Predicting the Initial Steps of Salt Stable CCMV Capsid Assembly with Atomistic Force Fields Zoltan Antal,† Janos Szoverfi,†,‡ and Szilard N. Fejer∗,† Provitam Foundation, 16 Caisului Street, Cluj-Napoca, Romania, and University Politehnica of Bucharest, Faculty of Applied Chemistry and Materials Science, 1-7 Gh. Polizu Street, Bucharest, Romania E-mail: [email protected]

Abstract Computational prediction of native protein-protein interfaces still remains a challenging task. In virus capsid proteins, each protein unit is in contact with copies of itself through several interfaces. The relative strength of the different contacts affects the dynamics of the assembly, especially if the process is hierarchical. We investigate the dimerization of the salt stable CCMV capsid protein using a combination of different computational tools. The best predictions of dimer configurations provided by blind docking with ZDOCK are rescored using geometry optimization with the Amber and Rosetta force fields. We also evaluate the relative stabilities of the three main interfaces present in the icosahedral capsid, using locally restricted docking with Rosetta. Both the rescoring and locally restricted docking results report a particularly stable protein-protein interface, which is the most likely intermediate during the first stage of ∗

To whom correspondence should be addressed Provitam Foundation, 16 Caisului Street, Cluj-Napoca, Romania ‡ University Politehnica of Bucharest, Faculty of Applied Chemistry and Materials Science, 1-7 Gh. Polizu Street, Bucharest, Romania †

1

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

the hierarchical capsid assembly. The blind docking results rescored with both Amber and Rosetta yield docking funnels, i.e. three or more near-native structures among the top five predictions. The result supports experimental observations for the in vitro assembly for CCMV capsids. The cross-validation of the results suggests that energy landscape-based methods with biomolecular force fields have the potential to improve existing docking procedures.

Introduction Computational studies on protein-protein interfaces have become recently increasingly important for understanding various biological processes, due to advances in accurate modelling, accessibility of computing power, and increase of available high-resolution 3D protein structures. Accurately predicting protein-protein binding modes is still a challenge, 1 and even assessing the quality of the predicted structures is not trivial. 2,3 Several main issues influence the quality of structure prediction for protein-protein complexes, among which the most important is accurate description of interactions, and sampling the relevant configuration space. Accuracy can be increased by using computationally expensive force fields, 4,5 and various efficient algorithms have been developed recently to enhance sampling of candidate dimer structures. 6,7 Combining these approaches can potentially enhance the quality of predictions, there are several recent examples of such method developments. 8–10 For the better understanding of self-assembling systems, 11 virus capsid proteins capable of spontaneous aggregation into well-defined structures are of particular importance. These proteins are special in the way that they probably evolved to have more than one specific interaction site. The blind prediction of preferred binding modes is therefore an interesting challenge. The cowpea chlorotic mottle virus (CCMV) capsid protein is capable to spontaneously aggregate into empty icosahedral capsids in vitro, under appropriate environmental conditions (pH, ionic strength). The experimentally determined mechanism for assembly is hierarchi2


Page 2 of 24

Page 3 of 24

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


cal, starting with dimerization, followed by formation of a pentamer of dimers, then possibly either cooperative addition of dimers to form an icosahedral T=3 capsid, or aggregation of pentamers to form pseudo-T=2 capsids for higher protein concentration or different pH. 12 The association free energy of the capsid proteins was determined to be about −2.2 kcal/mol, using data fitting on experimental light scattering measurements. 13 In this contribution we investigate the association of two salt stable CCMV capsid protein units, using a combination of monomer relaxation with Rosetta, 14 blind sampling of possible dimer structures with ZDOCK, and rescoring the predictions with both Amber and Rosetta. We also carry out locally restricted docking in order to assess the relative stabilities of the various experimentally relevant protein-protein interfaces. Amber force fields are physics-based, while Rosetta contains statistics-based potentials as well. We chose the two, conceptually different force fields in order to assess the robustness of our blind docking approach and that of relative stability estimates. Among the current, widely used blind docking methods, we selected ZDOCK for initial sampling of candidate docked structures because it carries out rigid body docking returning a large number of wellpacked dimers. We include a degree of flexibility through relaxation of initial structures, and subsequent relaxation of ZDOCK predictions, by using the two force fields, given that they are more physically realistic than simplistic energy functions used for taking into account some flexibility in fast blind docking methods. 15

The model system The best available experimental structure for the 190 residue CCMV capsid protein is that of a single-residue mutant in the N-terminal arm of wild type CCMV (K42R), with a resolution of 2.7 ˚ A (PDB id 1za7). This mutation stabilizes the virus structure compared to the wild type by forming extra hydrogen bonds between neighbouring residues along the quasi-sixfold vertices. 16 The structure of salt stable CCMV (ss-CCMV) is very similar to that of the wild

3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

type protein (RMSD < 0.6 ˚ A for the Cα atoms for residues 42-190). The asymmetric unit in the crystal structure contains three, nearly identical chains, labelled A, B and C. Among these, the position of residues 1-39 of chain A and that of residues 26-39 for chains B and C could not be determined, probably due to the flexibility of this chain segment. The disordered N-terminal region is involved in RNA packaging, 17 being located towards the interior of the capsid. This region is not important for capsid assembly, since a deletion mutant of the protein (residues 1-36 deleted) is also able to form icosahedral capsids. 18 However, the pathway for capsid assembly is likely different for empty capsids than that for the complete virus, with RNA (and the long N-terminal region) present. Our study is focused on investigating the association of capsid proteins which lead to the formation of empty capsids in experiments.

Methods Docking with Rosetta. The input dimer was first relaxed with the Rosetta Relax protocol, generating 10 decoys, out of which the lowest energy structure was used for the docking protocol. 20000 decoys were generated during the locally restricted docking and each new structure was compared against the starting structure. The protocol was applied to representative dimers of each protein-protein interface (CC, BC1 and AB). The docking protocol is described in more detail in the Supporting Information. Generating dimer structures with ZDOCK. Each separate chain from the asymmetric unit was relaxed with Rosetta to remove any clashes. The Relax protocol uses all-atom refinement with all-heavy-atom constraints in the Rosetta force field. 5,19 The most important conformational changes during the relaxation occurred on the flexible N-terminal part of the protein (residues 26 to 39), when relaxing chains that contained this sequence (B and C). We therefore used all possible pairwise permutations of the three chains as input into the ZDOCK server, 6,20 version 3.0.2, in order to maximize the diversity of input structures.

4


Page 4 of 24

Page 5 of 24

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


We did not complete any chains with the missing residues. ZDOCK uses a fast Fourier transformation-based algorithm to produce candidate docked structures from rigid protein units. 21 The best 2000 predictions for dimer structures of each pair (AA, AB, AC, BB, BC, CC) were saved and refined further with Amber and Rosetta. Binding energy calculation. The structures generated with ZDOCK were prepared for Amber energy minimization with tleap from AMBERTOOLS: hydrogen atoms were added and input coordinate and topology files were generated using the ff03 force field. 22 For solvation, we used a modified version of the generalized Born solvation model (igb=2), 23 and the mbondi2 radii set. Other parameters used: rgbmax = 8.23 ˚ A, salt concentration of 0.1 M and no cutoffs for the nonbonded interactions. The energy minimization of the dimer structures was performed using the L-BFGS minimizer, 24 as implemented in GMIN 25 with the Amber 12 software interfaced, 4 compiled for use on graphics cards. 26 The 6000 dimer structures which contain the fewest missing residues (BB, BC and CC) were minimized to an RMS force of 10−4 kcal/mol. Each optimized dimer structure was then separated and reoptimized, to estimate the total energy for the individual monomers. A typical geometry optimization took between 5-7 minutes on an nVidia Tesla K40 GPU, running an L-BFGS implementation for GPUs. 26 The binding energy for each dimer was calculated as the difference between the minimized potential energies for the associated and separated form of the dimers, respectively. We used an analogous scheme to rescore the ZDOCK predictions with Rosetta. In this case the binding energy for all 12000 structures was determined by optimizing the ZDOCK structures with the Relax protocol, and using the Prepack protocol on the relaxed structures in order to calculate the interface score. In principle, each dimer could be further optimized through Rosetta’s own docking protocol, as in Figure 1. However, the computational cost for such an analysis on all 12000 structures is rather large, given that it would involve generating and optimizing more than 200 million decoys in total. Structural alignment. The best alignment to the reference dimer structures was cal-

5



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

culated while taking into account the permutation of the two chains as well. Euclidean distances (D) and the corresponding RMSDs were calculated with the PERMOPT routine of GMIN, 25 using the alpha carbons of residues 40 to 190 in each protein chain. The total √ Cα RMSD was calculated as D/ N, where N = 302, the number of Cα atoms used in the analysis. For identifying interface residues, we used a surface area-based approach, as implemented in the InterfaceResidues Pymol script. 27 A residue is defined as being on the interface, if the difference between its solvent accessible area in the dissociated and associated forms is larger 2

than 1 ˚ A . Since each of the three main interface types had two or more variants (8 in total), we redefined the interface residues for each pair. For example, in case of the type I interface, the CC dimer had 68 residues at the interface, while the BA dimer had 69, according to the above definition.

Evaluating the relative stability of the main experimental protein-protein interfaces Starting from the full virus capsid structure generated by ViperDB, 28 one can identify three main interfaces for protein dimers in an icosahedral shell, labelled according to their interface areas (type I - CC, BA; type II - BC1, BC2, AA1, AA2; type III - AB, AC). The type I interface has the largest area, therefore it is logical to assume that this interface plays a role in primary dimer formation. 12,29 The type I interface is the only one among these three types which contains residues of the C-terminal part of the protein. Experimental data shows that a deletion mutant with the last 15 residues removed from the C-terminal is unable to form dimers, 18 and the capsid does not form either. However, it was not clear from the experiments whether dimers cannot form due to the missing interface residues or the deletion causes misfolding of the protein. The CC and BA type dimers are in a slightly different conformation, due to the position 6


Page 6 of 24

Page 7 of 24

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


they take in the T=3 icosahedral capsid (chain A is closest to the fivefold symmetry axis). This interface is described as a ‘hinge’, the angle between the two protein chains can take different values depending on the position in the capsid (and the type of capsid). For example, an N-terminal deletion mutant of the protein was found to form T=1, pseudo-T=2 and T=3 capsids as well. 30 The angle values around this interface were determined as 38◦ (BA), 42◦ (CC), and 45◦ (dimer in the T=1 capsid). In order to assess the relative stabilities of the three interface types, we subjected one dimer from each type to locally restricted docking, using the Rosetta Docking protocol. This protocol starts from a local perturbation of the input structure and carries out a lowresolution Monte Carlo search with adaptive rotational and translational steps, 31,32 followed by a high-resolution refinement of the generated lowest energy structure. The relative orientation of the two docking partners was kept fixed, so that we can explore the local conformational space around the input structure. 33 Figure 1 shows the interface score of each decoy, compared to the interface RMSD (iRMSD) with the native structure. The graphs for type I and type II interfaces are characteristic of docking funnels, 34 while that for the type III interface does not show a preference for a local binding mode. Type I interface has the lowest docking score, -21.09 Rosetta Energy Units (REU), followed by type II (-14.71 REU) and type III (-6.16 REU). The relative stability of the three interfaces is in a good agreement with the number of contacts in each interface. The docking funnels indicate a more flexible binding mode for the type I interface (wider funnel), compared to that of the type II interface. The observed increased flexibility of the binding mode might be necessary to accommodate the ‘bend’ of the B and A monomers at the pentameric site, relative to the CC dimer. Details of the applied protocol can be found in the Supporting Information. In order for the capsid assembly to be hierarchical, the formation of various oligomers (dimers, pentamers of dimers) must take place on different time scales, therefore the dimerization interaction has to be much stronger than that for subsequent association of dimers.

7



ts 0

−5 Interface score (REU)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 24

−10

−15

−20

−25

0

5

15

10

20

25

30

iRMSD (˚ A)

Figure 1: Locally restricted docking results starting from the three main interfaces in the complete 1za7 capsid. Red: CC (type I interface), green: BC (type II interface), blue: AB (type III interface). The three types of protein dimers are also shown. If we assume that the next step of assembly is the formation of pentamers, that step involves forming two extra interfaces per dimer, a type II and a type III interface. Given that we could not find a specific docking pose for the type III interface, it is possible that the formation of pentamers is driven exclusively by interactions along the type II interface. This conclusion is further supported by the locally restricted docking results for the dimerization of a type I dimer, shown in Figure 2. The starting structure (two neighbouring type I dimers) was also obtained from the full capsid, as generated by ViperDB. The results show the typical energy funnel related to a successful docking run, with the lowest interface energies ranging in the -14.0 REU area, with a docking funnel very similar to that for dimerization along the type II interface (Figure 1). We therefore observe no further stabilization of the tetramer due to

8


Page 9 of 24

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the formation of the type III interface. A similar idea of finding the most likely capsid assembly pathways was first proposed by Horton and Lewis, 35 through estimating free energies of binding for subsequent addition of protein units onto an icosahedral surface. More stable configurations were proposed to be likely intermediates for capsid assembly. However, the estimation was done by treating each protein as a rigid unit, and scaling the various components of the interaction free energy (nonpolar and polar solvation free energies, as well as association entropies) in order to obtain a good fit to experimental values. A more refined version of this approach was developed by Reddy et al., 36 in which association energies were estimated using the CHARMM force field. In this case, each subunit was relaxed, but with harmonic restraints on all heavy atoms. The deviation of each protein unit from the crystal structure was therefore minimal, so in this case the protein was treated essentially as a rigid unit as well. Clearly, association of capsid proteins does not occur on an icosahedral surface, therefore each unit has much more flexibility than in the final assembly. The parameterization of modern force fields is nowadays accurate enough to allow relaxation of each protein without any heavy atom restraints, and observed deviations from the crystal structure are not necessarily an artefact of the force field, especially for virus capsid proteins.

Blind docking of monomers In this section we present the protocol applied to the ss-CCMV monomer for estimating the best binding modes from using only structural information of the monomer chains. We treat computational protein-protein docking as being essentially a global optimization problem on a hierarchically organized energy landscape. 37 As such, blind protein-protein docking has two essential parts: sampling and ranking of structures. For sampling, we use the ZDOCK server, and for structure ranking we employ two atomistic force fields: Amber ff03 and Rosetta. The best 2000 structures from each ZDOCK

9



ts 0

Interface score (REU)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 24

−5

−10

−15

0

5

15

10

20

25

30

iRMSD (˚ A)

Figure 2: Locally restricted docking results starting from two adjacent type I dimers from the pentameric site of the 1za7 capsid, docking them around their common interfaces (a type II and type III interface). run are relaxed to the nearest local energy minimum, using these two force fields, and ranking is done according to the ‘binding energy’ terms. Table 1 summarizes the number of distinct dimers among the best 2000 predictions generated with ZDOCK for each pair. For clustering the structures, we used the MaxCluster software 38 with the hierarchical ‘average linkage’ scheme. 39 We also took into account the possible permutation of the two chains of the dimer when calculating the best Cα RMSDs between any two pairs of dimers. The flexible N-terminal part (residues 26 to 39) was excluded from the backbone RMSD calculation. The distance threshold used for clustering was 4 ˚ A. In total, out of 12000 best predictions, 2229 distinct clusters were identified with this clustering scheme.

10


Page 11 of 24

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 1: Clustering results for the structures generated with ZDOCK. Dimer Clusters

AA AB AC BB BC CC All 403 579 619 488 540 351 2229

Calculating binding energies It is important to note that we cannot assign any meaningful physical value to the binding energies calculated with our protocol either using Rosetta or Amber. The binding energy values in the arbitrary Rosetta Energy Units (Rosetta) or kcal/mol (Amber) are not related directly to the real binding affinities of the protein. Predicting real binding affinities is still a challenge, given the complexity of factors influencing the thermodynamics of proteinprotein aggregation, 1 such as the importance of using conformational ensembles, solvent effects, entropy contribution, crowding in cells etc. The binding energies we calculate can be compared for various dimer configurations of the same system, the force fields are acting as enhanced scoring functions with more physically realistic energy terms.

Structural alignment with experimentally relevant dimers We compared each structure optimized with Amber and Rosetta with each of the three main interfaces, calculating their best alignment to these. We also calculated the iRMSD for each of the optimized structures, relative to the protein-protein interfaces obtained from the icosahedral capsid. Euclidean distances were determined between the interface Cα atoms of the eight experimental dimer structures and all structures from ZDOCK reoptimized with both Amber and Rosetta. Figures 3a-c and 4a-c show the calculated binding energies for each structure as a function of the best iRMSD for the three interface types with these two force fields. The structure with the lowest Amber binding energy has an iRMSD of only 1.67 ˚ A to the experimental type I interface, and the three structures with the lowest binding energies have iRMSD values within 3 ˚ A to this interface. These values would account for a medium

11



b

0

0

−20

−20

−40

−40

Binding energy (kcal/mol)


a

−60 −80 −100

−60 −80 −100

−120

−120

−140

−140

−160

0

5

10

20

15

25

30

35

−160

40

0

5

10

15

iRMSD (˚ A)

d

0

35

−40

30

−60

25

−80 −100

−140

5

10

20

15

35

40

25

30

35

40

15 10

5

30

20

−120

0

25

40

−20

−160

20 iRMSD (˚ A)

total RMSD (˚ A)

c


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 24

25

30

35

40

0

0

5

iRMSD (˚ A)

10

15

20 iRMSD (˚ A)

Figure 3: a, Binding energies calculated with Amber, as a function of iRMSD of the structure from the experimental type I interface, b, type II interface, and c, type III interface, respectively. d, Correlation between iRMSD and total Cα RMSD. quality prediction, according to the rules of the Critical Assessment of PRedicted Interactions (CAPRI) experiment. 3,40 The total RMSD values for the top three structures are 1.77 ˚ A, A, respectively. The two structures with the lowest binding energies were A, and 2.86 ˚ 1.50 ˚ also the two best ranked in the ZDOCK run using chains B and C as input, while the other two sets of runs (BB and CC) did not produce highly ranked structures which would be experimentally relevant. In addition to the low-energy structures identified as type I interface, we also find a symmetric structure with a low binding energy (having a C2 axis), which does not align well with any of the three experimentally relevant interface types, suggesting a possible bias

12


Page 13 of 24

of the Amber force field for such symmetrical structures, but this bias is not significant enough to affect the quality of our predictions in this case. Although ZDOCK produced a few dimer structures similar to the type II interface, the Amber binding energies of these are not especially low, as shown in Figure 3b.

−5

−5

−10

−10


b0


a0

−15

−20

−25

−30

−15

−20

−25

0

5

10

15

20 iRMSD (˚ A)

25

30

35

−30

40

d30

−5

20 Rosetta interface score (REU)

c0


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


−10

−15

−20

5

10

15

20 iRMSD (˚ A)

25

30

35

40

10

0

−10

−20

−25

−30

0

0

5

10

15

20 iRMSD (˚ A)

25

30

35

40

−30

−150

−100 −50 Amber binding energy (kcal/mol)

0

50

Figure 4: a, Interface score calculated with Rosetta, as a function of iRMSD of the structure from the experimental type I interface, b, type II interface, and c, type III interface, respectively. d, Amber binding energies vs Rosetta interface scores. The two sets of lowest 50 structures for Amber and Rosetta have in common seven structures, coloured in orange. We found no clear correlation between binding energies calculated with Amber and interface scores calculated with Rosetta (Figure 4d). However, among the best 50 predictions for both force fields, seven structures are common, six of which are of type I interface. The dimer with the lowest Rosetta interface score is type I as well, and corresponds to the third best ranked ZDOCK prediction from chains B and C as input. 13



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 24

Rescoring the ZDOCK predictions with both Amber and Rosetta achieves a docking funnel, if we set the type I interface as the native interface of the protein dimer. The docking procedure has a funnel if there are three or more near-native structures among the top five predictions (N5 ≥ 3). 31 For the Amber predictions, N5 = 4, while the Rosetta scoring gives us N5 = 3. We emphasize that these are the results of global docking, and not locally restricted docking as in ref. 31 The interface scores as a function of iRMSD calculated to the three main experimental interface types (relaxed with Rosetta) are shown in Figure 4a-c. Any rescoring method can fail if the sampling used to generate the starting complex structures is not adequate. In our case, we could increase the number of predicted dimer configurations by using different combinations of nearly identical chains as input to ZDOCK, with the main structural difference between the chains being the orientation of the flexible N-terminal region relative to the core of the protein. If we consider only the ZDOCK scoring results, the performance of this method is very heterogenous: N5 = 0 for runs with chain combinations AA, AB, AC, BB, CC, while N5 = 5 for chain combination BC. Rescoring such results with different force fields and cross-matching the top predictions can narrow down the number of candidate structures, and also exclude sampling or force field artefacts. For blind docking, one has to strike a balance between extensive sampling of the most relevant configurational space and keeping the computational cost at a reasonable level. Approaches that involve initial sampling with fast Fourier transform-based methods therefore attempt to improve the energetics of the relevant basins of attraction for each predicted binding mode, for example by MM-PBSA/GBSA calculations, 41 or by structure clustering and subsequent rescoring with restricted local Monte Carlo searches. 42 There is no systemindependent, generally agreed protocol for rescoring docking predictions. When using MMPBSA/GBSA for rescoring, such calculations are not started from properly relaxed structures, usually involving only a fixed number of minimization steps and a sloppy convergence criterion before starting the production run. 10,41 Our approach allows for rescoring of a large number of predictions with relatively cheap computational methods, while we also take into

14


Page 15 of 24

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


account the flexibility of the protein units before rigid docking, by relaxing the monomer structures to different energy minima. In contrast to other methods, the ZDOCK predictions are reoptimized to proper energy minima using two force fields, Amber and Rosetta. We then compare the lowest energy structures predicted with both force fields, in order to minimize the errors necessarily introduced by the limited sampling and the use of energy minimization with no subsequent exploration of the basin surrounding each minimum. By comparing the low-energy structures, we also remove structures that are falsely predicted to be stable due to the particular force field used. This strategy is different from previous approaches of detecting false positives, and methods based on structure clustering 42 would benefit from rescoring local minima and evaluating MC trajectories with more than one force field. For CCMV, and virus capsid proteins in general, flexibility of parts of the protein is key for their incorporation in all of the different environments in an icosahedral shell. For protein-protein docking, the basic assumption is that the native structure has the largest basin of attraction among all other binding modes. 42,43 Additionally, the flexibility of the type I binding interface for CCMV, also described as a hinge in previous studies, 29 suggests that the free energy landscape is shallow in this region. On the other hand, the experimental hierarchical assembly mechanism points to a significant difference in interface stabilities. Our local docking results with Rosetta support both the increased flexibility and stability of the type I interface, relative to the type II interface. Also, ZDOCK sampled more type I interface structures than type II. Since more than one interface (and hence more than one binding mode) is present in virus capsids, predicting the aggregation of these proteins is a different problem from the usual ligand-receptor binding model. Our blind docking results cannot be compared with other approaches directly, since previous studies were carried out mostly on enzyme-inhibitor and antibody-antigen sets. Interestingly, the most widely used protein-protein docking benchmark set does not contain any homodimers among 230 ligandreceptor pairs of various difficulty, even in its recent version. 44 It remains to be seen how

15



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

our combined approach would perform on a standardized docking benchmark set. Since currently no universal methods exist that would fit all protein-protein docking problems, combined approaches like ours are worth exploring further.

Estimating relative binding affinities with MM-GBSA It is possible to estimate the relative change in binding free energy for the main interface types, compared to that of the most stable type I interface. We selected all type I and type II minima within an iRMS value lower than 3 ˚ A from our ensemble of minima, and estimated their binding free energies with MM-GBSA, following a protocol similar to that in ref. 10 The average association energy was determined as −185.7 ± 10.4 kcal/mol for type I (7 structures), and −108.2 ± 8.9 kcal/mol for type II (5 structures). The MM-GBSA results therefore show that the stability of the type II interface is about 0.6 times that of the type I interface. This is in a good agreement with the calculations done with Rosetta (Figure 1), only the energy scale is different: the type II interface has 0.7 times the stability of the type I interface (-14.71 REU vs. -21.09 REU), if we compare the most stable structures at the bottom of each funnel in Figure 1. Since no type III interfaces were sampled with ZDOCK according to the same iRMSD criteria, we minimized a type III A-B dimer from the full capsid. MM-GBSA results for two runs using different minimization conditions got an association energy of −76.5 ± 6.7 kcal/mol for the type III interface. However, the interface Cα RMSD increased by more than 12 ˚ A during the 6 ns production runs from both starting geometries, with no sign of convergence (Figure S1). This observation, and the lack of a clear binding mode for this interface using Rosetta suggests that the type III interface arises due to the icosahedral symmetry, and does not play a role in protein association. Moreover, the iRMS values during the MD runs are the smallest for all type I structures, followed by the type II structures, and the largest for the type III interface (Figure S1). The MM-GBSA and RosettaDock calculations are therefore consistent with the experimental observation of hierarchical capsid assembly. The Supporting Information 16


Page 16 of 24

Page 17 of 24

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


contains more details about the MM-GBSA protocol and results.

Conclusions Our results suggest that the most important interface for dimer formation is that of type I, and it is likely that this is the type of dimer formed in the initial stages of CCMV capsid assembly. Locally restricted docking with Rosetta shows a more flexible binding mode for the type I interface than that for type II, therefore the basin of attraction of the type I interface is larger, facilitating efficient dimerization. We also found that the combination of ZDOCK predictions with rescoring using both Amber and Rosetta was able to pinpoint the most favoured binding mode for two ss-CCMV capsid protein monomers with a very good accuracy. The second step of the hierarchical self-assembling process, which involves dimerization of type I dimers, has an interface score similar to that of the type II interface, although during this process one type II and one type III interface is formed. The type III interface can be therefore considered irrelevant for CCMV capsid assembly.

Acknowledgement This work was supported by a grant of the Romanian National Authority for Scientific Research and Innovation, CNCS – UEFISCDI, project number PN-II-RU-TE-2014-4-1176.

Supporting Information Available Details of calculations done with Rosetta, MM-GBSA protocol description and results, correlation between reoptimizing low-energy structures with the ff03.r1 and ff14SB force fields, link to structure files for the 6000 ZDOCK predictions minimized with Amber. This material is available free of charge via the Internet at http://pubs.acs.org/.

17



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 24

References (1) Kastritis, P. L.; Bonvin, A. M. J. J. On the Binding Affinity of Macromolecular Interactions: Daring to Ask Why Proteins Interact. J. R. Soc. Interface 2012, 10, 20120835. (2) Basu, S.; Wallner, B. DockQ: a Quality Measure for Protein-Protein Docking Models. PLoS One 2016, 11, e0161879. (3) Méndez, R.; Leplae, R.; De Maria, L.; Wodak, S. J. Assessment of Blind Predictions of Protein–Protein Interactions: Current Status of Docking Methods. Proteins: Struct., Funct., Bioinf. 2003, 52, 51–67. (4) Case, D. A.; Darden, T. A.; Cheatham, T. E.; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Walker, R. C.; Zhang, W.; Merz, K. M.; Roberts, B.; Hayik, S.; Roitberg, A.; Seabra, G.; Swails, J.; Goetz, A. W.; Kolossváry, I.; Wong, K. F.; Paesani, F.; Vanicek, J.; Wolf, R. M.; Liu, J.; Wu, X.; Brozell, S. R.; Steinbrecher, T.; Gohlke, H.; Cai, Q.; Ye, X.; Wang, J.; Hsieh, M. J.; Cui, G.; Roe, D. R.; Mathews, D. H.; Seetin, M. G.; Salomon-Ferrer, R.; Sagui, C.; Babin, V.; Luchko, T.; Gusarov, S.; Kovalenko, A.; Kollman, P. A. AMBER 12; University of California: San Francisco; 2012. (5) Leaver-Fay, A.; Tyka, M.; Lewis, S. M.; Lange, O. F.; Thompson, J.; Jacak, R.; Kaufman, K.; Renfrew, P. D.; Smith, C. A.; Sheffler, W.; Davis, I. W.; Cooper, S.; Treuille, A.; Mandell, D. J.; Richter, F.; Ban, Y.-E. A. E.; Fleishman, S. J.; Corn, J. E.; Kim, D. E.; Lyskov, S.; Berrondo, M.; Mentzer, S.; Popović, Z.; Havranek, J. J.; Karanicolas, J.; Das, R.; Meiler, J.; Kortemme, T.; Gray, J. J.; Kuhlman, B.; Baker, D.; Bradley, P. ROSETTA3: an Object-Oriented Software Suite for the Simulation and Design of Macromolecules. Methods in Enzymol. 2011, 487, 545. (6) Pierce, B. G.; Wiehe, K.; Hwang, H.; Kim, B.-H.; Vreven, T.; Weng, Z. ZDOCK

18


Page 19 of 24

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Server: Interactive Docking Prediction of Protein–Protein Complexes and Symmetric Multimers. Bioinformatics 2014, 30, 1771–1773. (7) Padhorny, D.; Kazennov, A.; Zerbe, B. S.; Porter, K. A.; Xia, B.; Mottarella, S. E.; Kholodov, Y.; Ritchie, D. W.; Vajda, S.; Kozakov, D. Protein–Protein Docking by Fast Generalized Fourier Transforms on 5D Rotational Manifolds. Proc. Natl. Acad. Sci. USA 2016, 113, E4286–E4293. (8) Kynast, P.; Derreumaux, P.; Strodel, B. Evaluation of the Coarse-Grained OPEP Force Field for Protein-Protein Docking. BMC Biophys. 2016, 9, 4. (9) Zhang, X.; Wong, S.; Lightstone, F. Toward Fully Automated High Performance Computing Drug Discovery: a Massively Parallel Virtual Screening Pipeline for Docking and Molecular Mechanics/Generalized Born Surface Area Rescoring to Improve Enrichment. J. Chem. Inf. Model. 2014, 54, 324. (10) Maffucci, I.; Contini, A. Improved Computation of Protein–Protein Relative Binding Energies with the Nwat-MMGBSA Method. J. Chem. Inf. Model. 2016, 56, 1692–1704. (11) Whitesides, G. M.; Grzybowski, B. Self-Assembly at All Scales. Science 2002, 295, 2418–2421. (12) Zlotnick, A.; Aldrich, R.; Johnson, J. M.; Ceres, P.; Young, M. J. Mechanism of Capsid Assembly for an Icosahedral Plant Virus. Virology 2000, 277, 450–456. (13) Xie, L.; Smith, G. R.; Feng, X.; Schwartz, R. Surveying Capsid Assembly Pathways through Simulation-Based Data Fitting. Biophys. J. 2012, 103, 1545–1554. (14) Conway, P.; Tyka, M. D.; DiMaio, F.; Konerding, D. E.; Baker, D. Relaxation of Backbone Bond Geometry Improves Protein Energy Landscape Modeling. Protein Sci. 2014, 23, 47–55.

19



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 24

(15) Torchala, M.; Moal, I. H.; Chaleil, R. A.; Fernandez-Recio, J.; Bates, P. A. SwarmDock: a Server for Flexible Protein–Protein Docking. Bioinformatics 2013, 29, 807–809. (16) Speir, J. A.; Bothner, B.; Qu, C.; Willits, D. A.; Young, M. J.; Johnson, J. E. Enhanced Local Symmetry Interactions Globally Stabilize a Mutant Virus Capsid that Maintains Infectivity and Capsid Dynamics. J. Virol. 2006, 80, 3582–3591. (17) Vriend, G.; Verduin, B.; Hemminga, M. Role of the N-terminal Part of the Coat Protein in the Assembly of Cowpea Chlorotic Mottle Virus: A 500 MHz Proton Nuclear Magnetic Resonance Study and Structural Calculations. J. Mol. Biol. 1986, 191, 453–460. (18) Zhao, X.; Fox, J. M.; Olson, N. H.; Baker, T. S.; Young, M. J. In Vitro Assembly of Cowpea Chlorotic Mottle Virus from Coat Protein Expressed in Escherichia coli and in Vitro-Transcribed Viral cDNA. Virology 1995, 207, 486–494. (19) Kaufmann, K. W.; Lemmon, G. H.; DeLuca, S. L.; Sheehan, J. H.; Meiler, J. Practically Useful: What the Rosetta Protein Modeling Suite Can Do for You. Biochemistry 2010, 49, 2987–2998. (20) Pierce, B. G.; Hourai, Y.; Weng, Z. Accelerating Protein Docking in ZDOCK Using an Advanced 3D Convolution Library. PloS One 2011, 6, e24657. (21) Vreven, T.; Hwang, H.; Weng, Z. Integrating Atom-Based and Residue-Based Scoring Functions for Protein–Protein docking. Protein Sci. 2011, 20, 1576–1586. (22) Duan, Y.; Wu, C.; Chowdhury, S.; Lee, M. C.; Xiong, G.; Zhang, W.; Yang, R.; Cieplak, P.; Luo, R.; Lee, T.; Caldwell, J.; Wang, J.; Kollman, P. A Point-Charge Force Field for Molecular Mechanics Simulations of Proteins Based on CondensedPhase Quantum Mechanical Calculations. J. Comp. Chem. 2003, 24, 1999–2012. (23) Onufriev, A.; Bashford, D.; Case, D. A. Modification of the Generalized Born Model Suitable for Macromolecules. J. Phys. Chem. B 2000, 104, 3712–3720. 20


Page 21 of 24

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(24) Liu, D. C.; Nocedal, J. On Limited Memory BFGS Method for Large Scale Optimization. Math. Prog. 1989, 45, 503. (25) Wales, D. J. GMIN: A Program for Finding Global Minima and Calculating Thermodynamic Properties from Basin-Sampling,. http://www-wales.ch.cam.ac.uk/GMIN, (accessed April 5, 2017). (26) Mantell, R. G.; Pitt, C. E.; Wales, D. J. GPU-Accelerated Exploration of Biomolecular Energy Landscapes. J. Chem. Theory Comput. 2016, 12, 6182–6191. (27) Vertrees, J. InterfaceResidues Pymol Script. https://pymolwiki.org/index.php/ InterfaceResidues, (accessed April 5, 2017). (28) Reddy, V. S.; Natarajan, P.; Okerberg, B.; Li, K.; Damodaran, K.; Morton, R. T.; Brooks, C. L.; Johnson, J. E. Virus Particle Explorer (VIPER), a Website for Virus Capsid Structures and Their Computational Analyses. J. Virol. 2001, 75, 11943–11947. (29) Globisch, C.; Krishnamani, V.; Deserno, M.; Peter, C. Optimization of an Elastic Network Augmented Coarse Grained Model to Study CCMV Capsid Deformation. PLoS One 2013, 8, e60582. (30) Tang, J.; Johnson, J. M.; Dryden, K. A.; Young, M. J.; Zlotnick, A.; Johnson, J. E. The Role of Subunit Hinges and Molecular Switches in the Control of Viral Capsid Polymorphism. J. Struct. Biol. 2006, 154, 59–67. (31) Chaudhury, S.; Berrondo, M.; Weitzner, B. D.; Muthu, P.; Bergman, H.; Gray, J. J. Benchmarking and Analysis of Protein Docking Performance in Rosetta v3.2. PLoS One 2011, 6, e22477. (32) Gray, J. J.; Moughon, S.; Wang, C.; Schueler-Furman, O.; Kuhlman, B.; Rohl, C. A.; Baker, D. Protein–protein Docking with Simultaneous Optimization of Rigid-Body Displacement and Side-Chain Conformations. J. Mol. Biol. 2003, 331, 281–299. 21



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(33) Braun, T.; Orlova, A.; Valeg˚ ard, K.; Lind˚ as, A.-C.; Schröder, G. F.; Egelman, E. H. Archaeal Actin from a Hyperthermophile Forms a Single-Stranded Filament. Proc. Natl. Acad. Sci. USA 2015, 112, 9340–9345. (34) Lyskov, S.; Chou, F.-C.; Conchir, S. .; Der, B. S.; Drew, K.; Kuroda, D.; Xu, J.; Weitzner, B. D.; Renfrew, P. D.; Sripakdeevong, P.; Borgo, B.; Havranek, J. J.; Kuhlman, B.; Kortemme, T.; Bonneau, R.; Gray, J. J.; Das, R. Serverification of Molecular Modeling Applications: The Rosetta Online Server That Includes Everyone (ROSIE). PLoS One 2013, 8, 1–11. (35) Horton, N.; Lewis, M. Calculation of the Free Energy of Association for Protein Complexes. Protein Sci. 1992, 1, 169–181. (36) Reddy, V. S.; Giesing, H. A.; Morton, R. T.; Kumar, A.; Post, C. B.; Brooks, C. L.; Johnson, J. E. Energetics of Quasiequivalence: Computational Analysis of Protein– Protein Interactions in Icosahedral Viruses. Biophys. J. 1998, 74, 546–558. (37) Wales, D. J. Energy Landscapes: Applications to Clusters, Biomolecules and Glasses; Cambridge University Press, 2003. (38) Herbert, A. MaxCluster: A Tool for Protein Structure Comparison and Clustering. http://www.sbg.bio.ic.ac.uk/~ maxcluster, (accessed April 5, 2017). (39) de Hoon, M. J.; Imoto, S.; Nolan, J.; Miyano, S. Open Source Clustering Software. Bioinformatics 2004, 20, 1453–1454. (40) Janin, J. Assessing Predictions of Protein–Protein Interaction: the CAPRI Experiment. Protein Sci. 2005, 14, 278–283. (41) Chen, F.; Liu, H.; Sun, H.; Pan, P.; Li, Y.; Li, D.; Hou, T. Assessing the Performance of the MM/PBSA and MM/GBSA Methods. 6. Capability to Predict Protein–Protein

22


Page 22 of 24

Page 23 of 24

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Binding Free Energies and Re-Rank Binding Poses Generated by Protein–Protein Docking. Phys. Chem. Chem. Phys. 2016, 18, 22129–22139. (42) Kozakov, D.; Schueler-Furman, O.; Vajda, S. Discrimination of Near-Native Structures in Protein–Protein Docking by Testing the Stability of Local Minima. Proteins: Struct., Funct., Bioinf. 2008, 72, 993–1004. (43) Camacho, C. J.; Weng, Z.; Vajda, S.; DeLisi, C. Free energy landscapes of encounter complexes in protein-protein association. Biophys. J. 1999, 76, 1166–1178. (44) Vreven, T.; Moal, I. H.; Vangone, A.; Pierce, B. G.; Kastritis, P. L.; Torchala, M.; Chaleil, R.; Jimnez-Garca, B.; Bates, P. A.; Fernandez-Recio, J.; Bonvin, A. M.; Weng, Z. Updates to the Integrated Protein–Protein Interaction Benchmarks: Docking Benchmark Version 5 and Affinity Benchmark Version 2. J. Mol. Biol. 2015, 427, 3031–3041.

23



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Graphical TOC Entry

24


Page 24 of 24

Predicting the Initial Steps of Salt-Stable Cowpea ... - ACS Publications

Recommend Documents