Human Pyruvate Dehydrogenase Complex E2 and ... - ACS Publications

Apr 22, 2016 - Subunits: New Models and Insights from Molecular Dynamics ... Institute of Bioprocess and Biosystem Engineering, Hamburg University of ...
4 downloads 0 Views 9MB Size
Article pubs.acs.org/JPCB

Human Pyruvate Dehydrogenase Complex E2 and E3BP Core Subunits: New Models and Insights from Molecular Dynamics Simulations Samira Hezaveh, An-Ping Zeng, and Uwe Jandt* Institute of Bioprocess and Biosystem Engineering, Hamburg University of Technology, Denickestrasse 15, 21071 Hamburg, Germany S Supporting Information *

ABSTRACT: Targeted manipulation and exploitation of beneficial properties of multienzyme complexes, especially for the design of novel and efficiently structured enzymatic reaction cascades, require a solid model understanding of mechanistic principles governing the structure and functionality of the complexes. This type of system-level and quantitative knowledge has been very scarce thus far. We utilize the human pyruvate dehydrogenase complex (hPDC) as a versatile template to conduct corresponding studies. Here we present new homology models of the core subunits of the hPDC, namely E2 and E3BP, as the first time effort to elucidate the assembly of hPDC core based on molecular dynamic simulation. New models of E2 and E3BP were generated and validated at atomistic level for different properties of the proteins. The results of the wild type dimer simulations showed a strong hydrophobic interaction between the C-terminal and the hydrophobic pocket which is the main driving force in the intertrimer binding and the core self-assembly. On the contrary, the C-terminal truncated versions exhibited a drastic loss of hydrophobic interaction leading to a dimeric separation. This study represents a significant step toward a model-based understanding of structure and function of large multienzyme systems like PDC for developing highly efficient biocatalyst or bioreaction cascades.



INTRODUCTION The pyruvate dehydrogenase complex (PDC) is a large macromolecular machine that shows unique and fascinating properties like self-assembly, high structural stability, coupled metabolic channelling within multistep reactions, and regeneration of cofactors. In brief, it mainly converts pyruvate into acetyl-CoA and thereby links the anaerobic glycolytic energy pathway to the aerobic TCA cycle.1 Mammalian and specifically human PDC (hPDC) is assembled from several copies of four subunit components: pyruvate decarboxylase (E1), dihydrolipoamide acetyltransferase (E2), dihydrolipoamide dehydrogenase (E3), and E3-binding protein (E3BP).2 However, PDCs from other species such as yeast or microbes lack the E3BP subunit, and prokaryotic PDC can show cubic instead of dodecahedric core structure. For example, Bacillus stearothermophilus (BST) inner core (PDB ID 1B5S)3 is formed by 20 identical building units, i.e., E2 trimers, while in hPDC both E2s and E3BPs form at least two distinct populations of trimeric building blocks in not yet resolved stoichiometry.4−6 Both E2 and E3BP are composed of lipoyl domains (LDs) or swinging arms, subunit binding domains (SBD), linker domains, and inner domains.7 The inner domains form the core, and the E1 heterotetramers and E3 homodimers bind noncovalently to the subunit binding domains of E2 and E3BP, respectively. The E2/E3BP core structure is unique with a © XXXX American Chemical Society

special assembly, and it plays an important role in the PDC activity as it bridges E1 and E3 via its swinging arms where it can visit the active sites of E1, E2, and E3, respectively. This type of enzymatic mechanism is usually referred as “substrate channelling” from which the reactions cascades benefit in several ways, such as less diffusion loss, less inhibition by intermediates, protection of intermediates, and recycling of cofactors. These properties, combined with the ability to selfassemble and the modular structure, make the E2/E3BP and swinging arm-based channeling mechanism a high value target for the engineering of novel multienzymatic reaction cascades.8−10 So far the available experimental studies provide neither the molecular structure of hPDC core nor its assembly mechanism nor the dynamics. Fortunately, computational material science can provide powerful tools to explore these phenomena with high accuracy. Computer simulations can complement the experimental work to study the large complex at a molecular level. To the best of our knowledge, so far no computational studies have been conducted on the hPDC or other PDCs mechanism and no molecular dynamics (MD) simulation Received: March 15, 2016 Revised: April 20, 2016

A

DOI: 10.1021/acs.jpcb.6b02698 J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B



studies were done on the core itself. This research work aimed to study hPDC core subunits with respect to their key role in enzyme functioning and substrate channeling. To begin we applied the Molecular Dynamics simulation at the atomistic level. Although the atomistic level is not preferable due to its limitations with respect to time and system size to study the protein complexes such as the PDC,9 it will be later used as a base model for building higher-scale models of the system, i.e., coarse-grained and meso-scale which allow exploring these systems on scales of different orders of magnitude in length and time. This will further enable us for the targeted manipulation and also to create the foundation for the synthesis of artificial complexes in the future. For this purpose we made homology models of hPDC core subunits based on a recent template11 and simulated based on the Amber03 force field. Using atomistic models, different dynamic and structural properties of the E2 and E3BP monomers were studied. PDC core in general has a unique assembly in which trimeric units can bind noncovalently by the so-called “sock-and-ball” joints or “anchoring effect,” i.e., a few key residues from a hydrophobic pocket interact with the C-terminal α-helix from a neighboring trimer.12 The anchoring occurs when the Cterminal of one monomer binds noncovalently in the hydrophobic pocket of the adjacent one. This assembled structure has been shown to form at pH 7.4 and is quite stable, and even altering the pH as low as 5.0 cannot affect or interrupt its stability.13 Peng et al.13 designed a pH-sensitive E2 with application in drug delivery by introducing histidine residues to the intra- and intertrimers key residues. From their results only the intertrimer mutants triggered the repulsive interactions and solubility of the core at pH 5.0 while the intratrimer mutants still maintained their structure. Their observation demonstrated that the closely coupled amino acids interactions at the intertrimer interfaces are determinant in the formation of the full core.13−15 Therefore, the assembly can only be disturbed if the intertrimer interaction gets hindered. One of the main reasons to study the core disassembly is to understand the complexity of PDC core and why nature set up such a huge system and whether the mini cores consisting of only trimers are still active. So far some experimental investigations have been done on truncation of some key residues from C-terminal in other species mainly with 20 identical trimers of E2s.16,17 However, to the best of our knowledge there has been no such work done on hPDC especially because of higher complexity of the core due to the presence of E3BP. Moreover, no detailed MD simulation based analysis on any PDC core from other species has been done so far. In this direction we performed dimer simulations of these subunits in order to understand the anchoring effect as one of the unique features in the PDC core assembly. This paper is organized as follows. The details of the homology modeling, the force field, and simulation setups are given in the Models and Methods section. In the Results section, we present for the first time the atomistic model and the simulation results of E2 and E3BP monomers for which different calculated properties of dynamic and structural properties of the proteins are discussed. Further, the anchoring effect will be shown by simulation of the dimeric interactions of E2 and E3BP for the wild type (WT) and truncated (t) versions. Finally, in the Discussion and Conclusion sections overall properties and stability of the models of hPDC subunits as well as dimeric interactions will be discussed and summarized.

Article

MODELS AND METHODS

Homology Modeling. Despite several available experimental data on the structure of the PDC subunits in other species, still the structure of the hPDC core is not known. To the best of our knowledge there are no crystal structures of E2 and E3BP in hPDC available. The difficulty in the crystallization of E2-E3BP core is due to the long flexible linking arms and regions on subunits and also the core structure itself.18 However, there is one E2 monomeric structure model based on the cryo-electron microscopy (cryo-EM) technique19 built based on two bacteria models: B. Stearothermophilus (BST) and Azotobacter vinelandii (AV) with PDB IDs of 1B5S3 and 1EAE,20 respectively, and it is available under the PDB ID of 3B8K.19 In this paper we refer to this model as E2-EM. We investigated the accuracy of this model by running 40 ns MD simulation with Amber03 force field. The model did not show good structural stability, which can be related to the high percentage of loops in this model. This is while the crystallographic structures of E2 monomers from other species with higher resolutions (e.g., BST and AV) always contain more stable regions (small B-factor values) especially close to the active site. Therefore, it is more likely that the human model also has a similar structure rather than a simple and very flexible loop region as in the E2-EM model for which the root-meansquare distance (rmsd) of backbone atoms turned out to be very high around 0.57 ± 0.01 nm even after excluding the N terminal flexible loop. On the other hand, as mentioned for this model only a monomeric structure has been reported. However, for our subsequent study to build a higher order hPDC core structure we needed a reliable template for which trimeric structure is also available. Hence we preferred to build a new model based on a recent trimeric template with a higher resolution rather than a monomeric one built on the previous lower resolution templates. Building trimers from the latter would need further optimizations and approximations, especially at the interfaces, which we considered not to be a good option. Therefore, for our study we built new homology models for both E2 (inner domain: 239 residues) and E3BP (inner domain: 229 residues) based on a recent template from Escherichia coli (E. coli) with PDB ID 4N7211 using the SWISS model.21−24 The Automated Mode was applied in which the amino acid sequence is used as an input data. In this method the suitable templates are identified based on Blast and HHblits and ranked according to their estimated quality.25,26 We chose the E. coli template because of its higher resolution, 2.25 Å, in comparison to the E2-EM model templates from BST and AV with 4.4 and 2.6 Å resolutions, respectively. The final EM effective resolution for tE2 (catalytic domain of E2) was reported to be around 8.8 Å.19 Table 1 shows the SWISS analysis on different templates used for modeling of human E2. The sequence alignments of all species are reported in Figure S1 in the Supporting Information. Simulation Setup. All MD simulations and analysis were performed using GROMACS software package version 5.0.2.27,28 with the Amber03 force field.29 The Visual Molecular Dynamics (VMD) program was used for the visualizations.30 The available protein structures and coordinates were downloaded from the Protein Data Bank (PDB). In the simulations, all-bonds were constraint using the LINCS algorithm.31 The integration time step was set to 2 fs. The temperature was maintained to the reference values using the V-rescale with B

DOI: 10.1021/acs.jpcb.6b02698 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

system, but in the case of E3BP WT the system was already neutral. E2 and E3BP Dimer. Dimers of E2 and E3BP were simulated for WT as well as C-terminal truncated versions of E2 and E3BP for 8 (t8) and 7 (t7) amino acids, respectively. Two monomers (WT or truncated ones) were positioned and orientated in a box according to the biological assembly34−36 of BST 60mer in order to study the interactions of the same in an assembled core; i.e., they are facing each other from their Cterminals. Monomers were positioned about 2 nm from each other measured as the distance between the C-terminal of one to the hydrophobic pocket of the other one. The system was then solvated with ∼70 000 TIP3 water molecules in a box of ∼13 nm/side. Distance of proteins to the boundary was ∼1 nm. The MD runs were performed at NPT ensemble for 40 ns at 310 K. Systems were neutralized by adding two Cl − counterions to the E2 WT and the truncated systems. In the case of E3BP systems, no counterion was required.

Table 1. Different Templates Used for Modelling of Human E2 template

E. coli (4N72)a

BST (1B5S)b

AV (1EAA)b

E2-EM (3B8K)

identity coverage technique resolution (Å) type of core

31.0% 96% X-ray 2.25 24-mer

32.6% 92% X-ray 4.4 60-mer

29.1% 96% X-ray 2.6 24-mer

100% 100% Cryo-EM 8.8 60-mer

a b

Template and its PDB ID used in this work to model E2 and E3BP. Templates and their PDB ID used to model E2-EM.

coupling time constant τ = 0.1 ps. The pressure was kept maintained at 1 bar using Parrinello−Rahman barostat with coupling constant p = 2 ps. The isothermal compressibility value used for the pressure coupling algorithm was 4.5 × 10−5 bar−1 for water. The long-range interactions were applied using the Particle Mesh Ewald (PME) method with a real space cutoff of 1.0 nm, a Fourier mesh spacing of 0.16 nm, and fourth-order interpolation. The Lennard-Jones interactions were calculated using a cutoff of 1.0 nm. The systems were first energy minimized using the steepest descent algorithm32 for 50 000 steps to eliminate any overlap or clash between the atoms. After energy minimization, all atoms were assigned velocities from the Maxwell−Boltzmann distribution at 310 K. All systems were equilibrated in two phases. The first phase was conducted under an NVT ensemble for at least 100 ps with position restraints on the heavy atoms of the protein leading to relaxation of the solvent molecules around it. After the system temperature was stabilized, the pressure was maintained by conducting another 2 ns simulation under NPT ensemble. After the systems were well equilibrated, position restraints were removed and production MD was performed. All the errors on the calculated properties have been evaluated using the block averaging method.33 E2 ad E3BP Single Monomer. Single monomers were simulated for wild types (WT) of E2 and E3BP, respectively. Each monomer was positioned in a box of ∼9 nm/side filled with ∼23 000 TIP3 water molecules. Distance of protein to the boundary was ∼1 nm. The MD runs were performed at NPT ensemble at 310 K, initially for 40 ns with an extension to another 200 ns. Counterions were added in order to neutralize the systems; i.e., one Cl− counterion was added to the E2 WT



RESULTS AND DISCUSSION Homology Model. New homology models were built for E2 and E3BP. Figure 1 shows the superimposed E2 and E3BP homology models based on E. coli in comparison to the previous E2-EM. E2 has an active site shown for both models in Figure 1 with key residues of His212-Arg213-Val214-Val215Asp216. Single Monomer. The structural stability of our models were examined by analyzing root mean square deviation (rmsd), root mean square fluctuation (rmsf), radius of gyration (Rg), and secondary structure analysis with respect to initial structure during the MD simulations of single monomers of E2 and E3BP. Root Mean Square Deviation (rmsd) and Fluctuation (rmsf). In order to check the overall stability of the structure, total rmsd was calculated on the backbone atoms and by excluding the N-terminal flexible loops as shown in Figure 2. The length of the N-terminal flexible loop was defined as 10 and 20 residues for E2 and E3BP, respectively. Rmsd analysis with N-terminal exclusion from the initial 40 ns simulation showed a good stability for E2 and E3BP for which the values turned out to be 0.19 ± 0.01 and 0.24 ± 0.05 nm, respectively. In order to confirm the stability of our models, the simulations were extended to 240 ns. Figure 2a, b shows rmsd as a function of time for both proteins. The average rmsd for E2 and E3BP

Figure 1. E2 (green) and E3BP (red) new homology models in comparison to E2-EM model (yellow). C and N-terminals are tagged along with the lower loop (L), upper loop (U) in E2-EM, and upper β-hairpin (UB) in E2 and E3BP. C

DOI: 10.1021/acs.jpcb.6b02698 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 4. E2 trimer for hPDC (top view). Active site (orange) is wrapped among several regions from the same monomer (green) and the interface. β-hairpin (green) is the UB region in Figure 1. For clarity only one monomer is colored and N terminal is tagged.

Figure 2. Total rmsd calculated for (a) E2 and (b) E3BP (red curves) and the N-terminal excluded one (black curves) in comparison to (c) those from 40 ns simulation of E2-EM.

Figure 5. RDF plots of OActive site-OWater for E2 and E2-EM monomers.

In order to get more details on the fluctuations of different residues, the rmsf analysis was performed. Rmsf values for the backbone atoms were initially calculated up to 40 ns as shown in Figure 3a, to be comparable to the results of E2-EM. The results of extended simulations up to 240 ns for E2 and E3BP are shown in Figure 3b. As is shown, the stabilities of the new E2 and E3BP models are highly preserved during the course of the simulation even up to 240 ns, while in the case of E2-EM several regions show high flexibility. The first region with high rmsf value for all models is the N-terminal, which is very flexible especially for E3BP. However, in the case of E3BP even the helix region right after its N-terminal has high rmsf value that can explain the higher rmsd of E3BP than that of E2 (Figure 2b). These regions are especially among those with very high crystallographic B-factor values in their template itself. Another important region exhibiting strong differences between E2-EM and the new models is located at the residue numbers 87−98. This region is basically a β-hairpin in E2 and E3BP (Figure 1, UB) the same as in E. coli, and is one of the regions located very close (about 1 nm away) to the active site of E2. As shown in Figure 4, this region is among those which sandwich the E2 active site and normally interacts with the Nterminal part of the adjacent monomer in a trimer by forming hydrogen bonds. But as is observed, the E2-EM model is very flexible in this region due to the existence of a loop with a very high rmsf value of about 0.76 nm, shown as the U region in

Figure 3. Rmsf of E2, E3BP for (a) from 10 to 40 ns in comparison to the E2-EM and (b) from 40 to 240 ns analysis.

after 40 ns (up to 240 ns) was calculated as 0.24 ± 0.02 and 0.41 ± 0.03 nm, respectively. These results show that the structures are still within a good range of stability especially for E2. However, for E3BP there is an increase in rmsd value around 40 ns, and after that it reaches a plateau. On the other hand, the rmsd value for the E2-EM model turned out to be 0.57 ± 0.01 nm. The instability of E2-EM is mainly related to the higher number of loops (apart from the N-terminal flexible loops) in comparison to the newly modeled E2. The existence of more loops in the catalytic domain can be inferred by comparing the total and N-terminal excluded rmsd values in Figure 2c. As shown in the figure the trend and value of the Nterminal excluded curve is not very different from its total rmsd. This is an indication of instability in catalytic domain itself. D

DOI: 10.1021/acs.jpcb.6b02698 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 6. Conformation of first cluster (red) for E2 in comparison to the initial energy minimized structure (yellow) at two different side views (a) and (b).

which is unfavorable for the intertrimer anchoring effect. This was observed from our simulations and will be discussed later in the results of the dimer simulations. For all proteins there is a big rmsf peak around residue numbers from 180 to 190. This area is related to the lower loop (L) region (Figure 1) with a longer flexible loop connected to the β-sheets for E3BP (10 residues) than E2 (4 residues). From the X-ray structure this region has also high crystallographic Bfactor values. The higher flexibility of this region for E3BP in comparison to E2 is clearer from the higher and wider peak after 40 ns until 240 ns (Figure 3b). This region for E2-EM is mainly a loop with no β-sheets that leads to a higher rmsf value (∼0.64 nm) than in E2 (∼0.52 nm) within the 40 ns of simulations. The last region of interest is the active site, which is located in the residues from 212 to 216 (Figure 4), which is clearly less stable in E2-EM (rmsf ∼ 0.30 nm) than in E2 new model (rmsd ∼ 0.16 nm). As discussed, this can be related to the very flexible upper loop in E2-EM. Moreover, further investigations

Table 2. Cluster Analysis Details for E3BP cluster no.

no.

%

1 2 3 4

9714 5945 1548 964

51 31 8 5

Figure 1, rather than a stable β-hairpin for proper wrapping of the active site that can also affect the active site stability. From the figure the position of the loop is completely open and tilted backward in contrast to the β-hairpin in E2, which is right above the active site. The next important region is around the residue numbers 120−130 where the upper part of the hydrophobic pocket is located. This area is very stable for E2 (rmsf 0.22 nm) and E3BP (rmsf 0.10 nm). Obviously the higher flexibility and major unfolding (will be shown later) of this region for E2-EM (rmsf about 0.39 nm) makes the hydrophobic pocket less stable

Figure 7. Conformation of the first cluster (red) and second cluster (green) for E3BP in comparison to the initial energy minimized structure (yellow) at two different side views (a) and (b). E

DOI: 10.1021/acs.jpcb.6b02698 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B Table 3. Average Percentage of Secondary Structure Elements %

simulation time (ns)

structure

coil

β-sheet

β-bridge

bend

turn

α-helix

5 helix

3 helix

E2-EM E2 E3BP

40 240 240

53 68 63

31 21 23

15 23 21

0 2 1

12 7 10

10 10 10

28 33 31

4 0 0

0 4 4

they mainly interact via hydrophobic interactions. These investigations suggest that the active site stability is highly dependent on the presence of these intraprotein interactions, and in their absence an alteration in the active site stability can occur. Figure 5 shows the radial distribution function (RDF) for OActive site−OWater. From the plot, E2-EM active site is exposed to more water density especially in 1 nm distance that is shown by higher RDF values. This can prove the importance of the β-hairpin in the active site proper wrapping that leads to its better stability. It should be noted that the conformations in a single monomer may differ from its trimer due to some interfacial interactions which exist in the trimeric form, e.g., the Nterminal and the lower flexible loop that are directly located in the interface of the trimer. As discussed and was shown in Figure 3, these regions are very flexible in all models even in our new models (E2 and E3BP) in their monomeric form but can be stabilized via trimerization. On the other hand, there are some important regions far from the interface, e.g., the regions that form the hydrophobic pocket, and their stability is important in the self-assembly of the core. These areas are in fact in the outer regions of the trimers and in the water interface and are less likely to be affected by the trimerization. From our results in the E2-EM model, we observed low stability of these regions (Figure 3a) that can lead to lack of anchoring stability. This issue will be discussed later in the dimer simulation results. Cluster Analysis. The cluster analysis was used to detect conformational diversity of the protein structure. For this purpose we used the GROMOS clustering algorithm37 on the backbone atoms with a chosen cutoff of 0.15 nm; i.e., similar conformations with an rmsd less than 0.15 nm are added in a cluster, and finally the conformation with the smallest rmsd from other members of the cluster is considered as the representative structure of that cluster. The calculations were done after 50 ns. The total numbers of clusters are 7 and 13 for

Figure 8. Rg as a function of time for E3BP and E2.

Figure 9. Top view of E2-EM WT dimer at (a) 0 ns and (b) 40 ns. The big arrows show the position of the α-helix, which is an upper part of the hydrophobic pocket. N and C-terminals are tagged as N and C, respectively.

in this direction showed intraprotein interactions in the vicinity of E2 active site, which are missing in E2-EM. In E2, the active site is buried under the upper β-hairpin (about 1 nm away), and

Figure 10. E2 WT dimer simulation, (a) and (b) at 0 ns top and side view, respectively. (c) and (d) showing side view at 40 ns. F

DOI: 10.1021/acs.jpcb.6b02698 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 11. E2-t8 dimer simulation, (a) top view at 0 ns. (b), (c), and (d) showing side views of same orientation fixed on right monomer at 0, 5−10, and 20−40 ns, respectively. The left monomer gradually moves up from position 1 to 3. In (d) the lower loop of the top monomer directed itself into the page.

Figure 12. E3BP-t7 dimer simulation, (a) top view at 0 ns. (b), (c), and (d) showing side views of same orientation fixed on right monomer at 0, 5, and 40 ns, respectively. The left monomer gradually moves around the other one and rotates from position 1 to 3.

between the first two clusters in E3BP is due to the higher flexibility of the lower loop region which leads to a partial unfolding of the β-sheet in the case of second cluster as shown in the figure. Secondary Structure (SS) Analysis. DSSP (Define Secondary Structure of Proteins) method38 was used to derive the secondary structure of the two subunits. The secondary structure remains fairly conserved during the course of the simulations for E2 and E3BP as shown in contrast to the E2EM (Figure S3 in the Supporting Information) where the number of structures, i.e., α-helix, β-sheet, bridges, and turns, decreases while the number of coils increases. Table 3 gives the percentage of secondary structure for E2-EM in comparison to

E2 and E3BP, respectively. The total number of structures for matrix was 19 001 for both. In the case of E2, the first cluster accounts for 98% itself. Figure 6 shows the representative conformation of the first cluster for E2 in comparison to the initial structure. The cluster numbers as a function of time are reported in Figure S2 in the Supporting Information. In the case of E3BP the first two clusters account for 51% and 31%. Table 2 summarizes the cluster information for the same. Figure 7 shows the superimposed representative conformations of the first two clusters for E3BP. E2 showed the least conformational diversity during the simulation. However, there is a higher diversity in the case of E3BP. The major difference G

DOI: 10.1021/acs.jpcb.6b02698 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

are fairly stable especially for E2. In the case of E3BP there is a partial unfolding of α-helix into turn which occurs partially in the N-terminal and the lower loop area (L), but the rest is very stable. Further to identify the tertiary structure stability of the proteins, the average minimum distance matrix was calculated and plotted. The results showed a loss of tertiary structure stability for E2-EM within 40 ns while E2 and E3BP are stable up to 240 ns. The plots are reported in the Supporting Information in Figure S5 in comparison to the ones from E2EM in Figure S6. Radius of Gyration (Rg). The Rg values as a function of time are presented for both subunits in Figure 8. The Nterminal flexible loop was excluded from the calculations. As shown in the figure, the Rg for E2 is very constant and reaches a plateau after ∼20 ns at the average value of 1.82 ± 0.006 nm. In the case of E3BP, however, up to 180 ns the average Rg value is 1.81 ± 0.003 nm but in the last 70 ns of the simulation, the Rg of the protein decreased by ∼5% to an average value of 1.75 ± 0.001 nm. This reduction as observed can be due to a minor bending of the partially unfolded α-helix located right after the N-terminal. Dimer Simulation. In order to study the anchoring effect at the molecular level and study the effect of the C-terminal on this interaction, we performed simulation of E2-EM, E2, and E3BP dimers for the WT and for the aforementioned Cterminals truncated systems. As mentioned before due to the unfolding of the upper part in the hydrophobic pocket (Figure 3a), the anchoring effect in E2-EM can be retarded and weakened, and at the same time unfoldings in other parts of the protein was observed. Figure 9 shows the E2-EM dimer at 0 ns and after 40 ns. The Cterminals finally bend outward and unfold while the two monomers are lost the original orientation. This can be seen by a qualitative comparison in the angle difference between the big black arrows on the left monomers in Figure 9a, b. The geometry of the anchoring is completely lost, and even the αhelix of the hydrophobic pocket in the right monomer is unfolded as shown in Figure 9b, right monomer. The simulation results for E2 dimers showed that the initial orientation (Figure 10a, b) of the monomers with respect to each other is quite stable and the C-terminal anchoring in the hydrophobic pocket is well maintained up to 40 ns, as shown in Figure 10c, d. Here no major unfolding was detected as in E2EM. Similar results were observed for E3BP dimers (Figure S7 in the Supporting Information). The elimination of the C-terminal that leads to the loss of the hydrophobic interaction showed a dimer separation for E2-t8 in Figure 11. Monomers initial orientation as in Figure 11a, b was not conserved, and the monomers started to tilt and reorient themselves and lost the original orientation due to lack of the anchoring effect (Figure 11c, d). Although they gradually got separated, they are still attracted to each other from α2 and α3 from 5 to 20 ns. These results are compatible with the experimental results of Thermoplasma acidophilum bacteria by Marrott et al.16 They studied the effect of truncation on the core assembly and activity of the enzyme after minimizing the size of the system. After truncation of five residues from the Cterminal, their crystallographic data showed the presence of hexamers in which the trimers where attracted from same regions. After 20 ns and up to 40 ns the monomers were only attracted from their upper β-hairpins (UB region in Figure 1) with a very small contact area.

Figure 13. Average distance between Glu127 to Leu231 in E2 (a) and Gln115 to Leu222 in E3BP (b) for truncated and WT versions from both sides (1 and 2).

Figure 14. Number of contacts for hydrophobic atoms in a dimeric interaction for E2 WT and E2-t8 (a) and E3BP WT and E3BP-t7 (b).

E2 and E3BP. The number of structure especially for α-helix, βsheet in E2-3B8K (28% helical and 15% β-sheet) is much lower than that of E2 (33% helical and 23% β-sheet) and E3BP (31% helical and 21% β-sheet). Instead the number of coils in E2− 3B8k (31%) is higher than that in E2 (21%) and E3BP (23%), and even it increases during only 40 ns simulation time, which is an indicator of an instable catalytic domain. To be more specific, the evolution of secondary structure per residue along the trajectory was also reported as a function of time (Figure S4 in the Supporting Information). The structures H

DOI: 10.1021/acs.jpcb.6b02698 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B Table 4. Number of Contacts (within 0.6 nm Distance) and Minimum Distance for Hydrophobic Interactions E2

BP

WT

a

minimum distance (nm) number of contacts

0.33 ± 0.003 57.92 ± 1.96

minimum distance (nm) number of contacts

0.33 ± 0.001 48.94 ± 1.07

t8 Up to 5 ns 0.39 ± 0.01 5.35 ± 0.58 After 10 ns 0.39 ± 0.03b 10.71 ± 1.93b

WT

t7

0.36 ± 0.002 15.03 ± 3.84

0.92 ± 0.05 0.00a

0.36 ± 0.003 18.61 ± 0.40

0.65 ± 0.09b 1.05 ± 0.31b

Less than 0.02. bNot related to intertrimer interface.

Figure 15. Hydrophobic surface of residues (red and green) in E2 dimer for side (a) and top views (b) and in E3BP-t7dimer for side (c) and top views (d) at 0 ns. The blue color shows other residues in a transparent mode.

Hydrophobic Interactions. Since understanding the protein interaction is required to elucidate the anchoring effect, the hydrophobic interactions in a dimer were studied. These interactions were quantified for E2 and E3BP as shown in Figure 14 in which the number of hydrophobic interactions were measured before and after elimination of the C-terminal. These interactions were evaluated by calculation of the number of contacts between the hydrophobic heavy atoms (nonhydrogen) in the interface of the two monomers and for the distances less than 0.6 nm. Here it should be noted that the number of contacts for the truncated versions in the intertrimer interface should be considered only up to 5 ns when the two monomers more or less face each other from the C-terminals as in WT (intertrimer interface in core) after which they reorient and lose the interface interaction and calculations are not comparable with WT results. Within 5 ns of simulation time, the contact number is remarkably smaller for E2-t8 and especially for E3BP-t7 in comparison to their WTs. The details for the number of contacts and minimum distance for the hydrophobic interactions are shown in Table 4. Calculations were done up to 5 ns and then after 10 ns. The number of hydrophobic atoms in E2 and E3BP WTs are 480 and 402, respectively. To be more specific, E2 has five hydrophobic residues (234Pro, 235Ile, 237Met, 238Leu, 239Leu), while E3BP has four residues (225Pro, 226Ile, 228Leu, 229Ala) in their C-terminal. For this reason the number of contacts in E2 WT is higher than that in the E3BP

In the case of E3BP-t7 the monomers initial orientation (Figure 12a, b) was also lost, and the dimer separation occurred by losing the connection from one C-terminal and tilting and orienting itself backward while still connected from side parts as shown in Figure 12c, d. It showed a kind of rotation around the right monomer by separating from one C-terminal. Figure 13a shows the average distance between the hydrophobic pocket residue (Glu127) and the end terminal residue of E2-t8 (Leu231) in comparison to the same in their WT version. The distance between the targeted residues was very well maintained (1.97 ± 0.02 and 1.91 ± 0.04 nm for both sides) in the WT, which indicates an effective and strong anchoring. However, in the case of E2-t8 there was a rapid increase in the distance (for the last 20 ns it reached 3.55 ± 0.06 and 3.63 ± 0.11 nm and did not reach a plateau). This indicates a separation from both hydrophobic pockets. Similar results were observed for E3BP WT and E3BP-t7 by calculating the average distance between hydrophobic pocket residue (Gln115) and the end terminal residue of E3BP-t7 (Leu222). From Figure 13b, the WT versions are strongly bound from both sides (with averaged distance of 1.62 ± 0.01 and 1.74 ± 0.03 nm), but upon truncation of the C-terminals the distance started to increase for one side (reached average value of 3.40 ± 0.02 nm) while it decreased for the other side (reached average value of 0.94 ± 0.08 nm). This was visually shown before in Figure 12 where separation occurred while both monomers completely lost their original orientation with respect to each other that indicated the instability of the E3BPt7 dimer. I

DOI: 10.1021/acs.jpcb.6b02698 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

much higher stability (rmsf = 0.16 nm) than that of E2-EM (rmsf = 0.30 nm). Further the protein interactions and the anchoring effect were studied in dimers of E2 and E3BP where the monomers were oriented as in a 60-meric core to mimic the same interactions. The anchoring effect was observed by simulation of wild-type (WT) and truncated E2 (E2-t8) and E3BP (E3BP-t7) versions. In both WT subunit dimers, the original orientation was well preserved and the C-terminals were tightly bound to the hydrophobic pockets. On the contrary, upon elimination of the C-terminals, within a short simulation time the two monomers lost their mutual orientation in both cases and detached from their binding pockets. The detachment was shown by the distance increase between the C-terminals and the hydrophobic pockets in E2-t8 and E3BP-t7 in comparison to those in WT. Further, a dramatic decrease in the number of hydrophobic interactions in the truncated versions proved a dominant hydrophobic force in the dimers interaction and in the core assembly. This phenomenon is shown for the first time from our simulations at the molecular level for hPDC subunits. The results correlate well with earlier experimental investigations on other PDC species with 20 identical E2 containing trimers; for human PDC however, no experimental proof was reported so far. In our group experimental work is ongoing in this direction, and our preliminary results indicate the disassembly of hPDC for several double truncations, e.g. , E2-t8 and E3BP-t7, for the first time (data not shown). Our atomistic models will serve as a starting point for building higher-scale models of the hPDC core and complex system, i.e., at coarse-grained and meso scale, which will allow exploring them on a scale of different order of magnitude in size and time for the design of multienzyme systems. This will further enable targeted manipulation of such complex enzyme system and thus create the foundation for the synthesis of artificial complexes and enzyme cascade in the future.

Figure 16. Minimum distance for hydrophobic atoms in a dimeric interaction for E2 (a) and E3BP (b). Black curves are fitted to the Gaussian distribution.



WT version. At the same time the minimum distance is only slightly smaller in E2 than in E3BP. For the truncated versions the minimum distance for E3BPt7 is bigger than E2-t8 due to the presence of residues on the upper regions of the dimer interface which are absent in E3BPt7. Figure 15 shows the residues hydrophobic surface for E2-t8 and E3BP-t7. The minimum distance distributions of residues with hydrophobic interactions are shown in Figure 16. As seen before, the minimum distance between the two hydrophobic groups increased especially in the case of E3BP-t7 (Figure 15) and shows a wider Gaussian distribution in comparison to the E2-t8. This shows less hydrophobic interactions and less affinity of E3BP-t7 monomers than E2-t8.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.6b02698. Sequence alignments, cluster numbers, secondary structure changes, and minimum distance matrix analysis (PDF)



AUTHOR INFORMATION

Corresponding Author

*Fax: +49 (0)40 42878 2909. Tel.: +49 (0)40 42878 2847. Email: [email protected].



Notes

CONCLUSION AND OUTLOOK In this study we have presented new homology models of E2 and E3BP monomers of human pyruvate dehydrogenase complex (hPDC) core. We showed that the previous E2-EM model is not stable and not accurate enough. To validate our new models, extensive MD simulations were conducted on E2 and E3BP. For a more quantitative and mechanistic assessment several dynamic and structural properties such as rmsd, rmsf, Rg, and secondary structure analysis were carried out. Overall, both models, especially the E2 model, showed very good structural stability and consistency within the extended simulation time in comparison to the low-resolution cryo-EM model (E2-EM). The active site of the new E2 model showed a

The authors declare no competing financial interest.



ACKNOWLEDGMENTS We acknowledge funding by the German Federal Ministry of Education and Research (BMBF, grant 031A128).



REFERENCES

(1) Patel, M. S.; Patel, U.; Sen, A.; Sethia, G. C.; Hens, C.; Dana, S. K.; Feudel, U.; Showalter, K.; Ngonghala, C. N.; Amritkar, R. E. Experimental observation of extreme multistability in an electronic system of two coupled Rossler oscillators. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 2014, 89 (2), 022918.

J

DOI: 10.1021/acs.jpcb.6b02698 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B (2) Brautigam, C. A.; Wynn, R. M.; Chuang, J. L.; Chuang, D. T. Subunit and catalytic component stoichiometries of an in vitro reconstituted human pyruvate dehydrogenase complex. J. Biol. Chem. 2009, 284 (19), 13086−13098. (3) Izard, T.; Aevarsson, A.; Allen, M. D.; Westphal, A. H.; Perham, R. N.; de Kok, A.; Hol, W. G. J. Principles of quasi-equivalence and Euclidean geometry govern the assembly of cubic and dodecahedral cores of pyruvate dehydrogenase complexes. Proc. Natl. Acad. Sci. U. S. A. 1999, 96 (4), 1240−1245. (4) Wallis, N. G.; Allen, M. D.; Broadhurst, R. W.; Lessard, I. A. D.; Perham, R. N. Recognition of a surface loop of the lipoyl domain underlies substrate channelling in the pyruvate dehydrogenase multienzyme complex. J. Mol. Biol. 1996, 263 (3), 463−474. (5) Vijayakrishnan, S.; Callow, P.; Nutley, M. A.; McGow, D. P.; Gilbert, D.; Kropholler, P.; Cooper, A.; Byron, O.; Lindsay, J. G. Variation in the organization and subunit composition of the mammalian pyruvate dehydrogenase complex E2/E3BP core assembly. Biochem. J. 2011, 437, 565−574. (6) Brautigam, C. A.; Wynn, R. M.; Chuang, J. L.; Machius, M.; Tomchick, D. R.; Chuang, D. T. Structural insight into interactions between dihydrolipoamide dehydrogenase (E3) and E3 binding protein of human pyruvate dehydrogenase complex. Structure 2006, 14 (3), 611−621. (7) Perham, R. N. Swinging arms and swinging domains in multifunctional enzymes: Catalytic machines for multistep reactions. Annu. Rev. Biochem. 2000, 69, 961−1004. (8) Zhang, Y. H. P. Substrate channeling and enzyme complexes for biotechnological applications. Biotechnol. Adv. 2011, 29 (6), 715−725. (9) Jandt, U.; You, C.; Zhang, Y. H. P.; Zeng, A. P. Compartmentalization and metabolic channeling for multienzymatic biosynthesis: Practical strategies and modeling approaches. Adv. Biochem. Eng./Biotechnol. 2013, 137, 41−65. (10) Fu, J. L.; Yang, Y. R.; Johnson-Buck, A.; Liu, M. H.; Liu, Y.; Walter, N. G.; Woodbury, N. W.; Yan, H. Multi-enzyme complexes on DNA scaffolds capable of substrate channelling with an artificial swinging arm. Nat. Nanotechnol. 2014, 9 (7), 531−536. (11) Wang, J. J.; Nemeria, N. S.; Chandrasekhar, K.; Kumaran, S.; Arjunan, P.; Reynolds, S.; Calero, G.; Brukh, R.; Kakalis, L.; Furey, W.; et al. Structure and function of the catalytic domain of the dihydrolipoyl acetyltransferase component in Escherichia coli pyruvate dehydrogenase complex. J. Biol. Chem. 2014, 289 (22), 15215−15230. (12) Mattevi, A.; de Kok, A.; Perham, R. N. The pyruvate dehydrogenase multienzyme complex. Curr. Opin. Struct. Biol. 1992, 2 (6), 877−887. (13) Peng, T.; Lim, S. Trimer-based design of pH-responsive protein cage results in soluble disassembled structures. Biomacromolecules 2011, 12 (9), 3131−3138. (14) Lim, S.; Peng, T.; Sana, B. Protein cages as theranostic agent carriers. Springer 2013, 39, 321−324. (15) Peng, T.; Lee, H.; Lim, S. Isolating a trimer intermediate in the self-assembly of E2 protein cage. Biomacromolecules 2012, 13 (3), 699−705. (16) Marrott, N. L.; Marshall, J. J. T.; Svergun, D. I.; Crennell, S. J.; Hough, D. W.; van den Elsen, J. M. H.; Danson, M. J. Why are the 2oxoacid dehydrogenase complexes so large? Generation of an active trimeric complex. Biochem. J. 2014, 463, 405−412. (17) Koike, K.; Suematsu, T.; Ehara, M. Cloning, overexpression, and mutagenesis of cDNA encoding dihydrolipoamide succinyltransferase component of the porcine 2-oxoglutarate dehydrogenase complex. Eur. J. Biochem. 2000, 267 (10), 3005−3016. (18) Vijayakrishnan, S.; Kelly, S. M.; Gilbert, R. J. C.; Callow, P.; Bhella, D.; Forsyth, T.; Lindsay, J. G.; Byron, O. Solution structure and characterisation of the human pyruvate dehydrogenase complex core assembly. J. Mol. Biol. 2010, 399 (1), 71−93. (19) Yu, X. K.; Hiromasa, Y.; Tsen, H.; Stoops, J. K.; Roche, T. E.; Zhou, Z. H. Structures of the human pyruvate dehydrogenase complex cores: A highly conserved catalytic center with flexible N-terminal domains. Structure 2008, 16 (1), 104−114.

(20) Mattevi, A.; Obmolova, G.; Kalk, K. H.; Teplyakov, A.; Hol, W. G. J. Crystallographic analysis of substrate binding and catalysis in dihydrolipoyl transacetylase (E2p). Biochemistry 1993, 32 (15), 3887− 3901. (21) Arnold, K.; Bordoli, L.; Kopp, J.; Schwede, T. The SWISSMODEL workspace: A web-based environment for protein structure homology modelling. Bioinformatics 2006, 22 (2), 195−201. (22) Biasini, M.; Bienert, S.; Waterhouse, A.; Arnold, K.; Studer, G.; Schmidt, T.; Kiefer, F.; Cassarino, T. G.; Bertoni, M.; Bordoli, L.; Schwede, T.; et al. SWISS-MODEL: Modelling protein tertiary and quaternary structure using evolutionary information. Nucleic Acids Res. 2014, 42 (1), 252−258. (23) Bordoli, L.; Kiefer, F.; Arnold, K.; Benkert, P.; Battey, J.; Schwede, T. Protein structure homology modeling using SWISSMODEL workspace. Nat. Protoc. 2009, 4 (1), 1−13. (24) Guex, N.; Peitsch, M. C.; Schwede, T. Automated comparative protein structure modeling with SWISS-MODEL and SwissPdbViewer: A historical perspective. Electrophoresis 2009, 30, S162− S173. (25) Altschul, S. F.; Madden, T. L.; Schaffer, A. A.; Zhang, J. H.; Zhang, Z.; Miller, W.; Lipman, D. J. Gapped BLAST and PSI-BLAST: A new generation of protein database search programs. Nucleic Acids Res. 1997, 25 (17), 3389−3402. (26) Remmert, M.; Biegert, A.; Hauser, A.; Soding, J. HHblits: Lightning-fast iterative protein sequence searching by HMM-HMM alignment. Nat. Methods 2011, 9 (2), 173−175. (27) Van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. GROMACS: Fast, flexible, and free. J. Comput. Chem. 2005, 26 (16), 1701−1718. (28) Lindahl, E.; Hess, B.; van der Spoel, D. GROMACS 3.0: A package for molecular simulation and trajectory analysis. J. Mol. Model. 2001, 7 (8), 306−317. (29) Duan, Y.; Wu, C.; Chowdhury, S.; Lee, M. C.; Xiong, G. M.; Zhang, W.; Yang, R.; Cieplak, P.; Luo, R.; Lee, T.; et al. Point-charge force field for molecular mechanics simulations of proteins based on condensed-phase quantum mechanical calculations. J. Comput. Chem. 2003, 24 (16), 1999−2012. (30) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual molecular dynamics. J. Mol. Graphics 1996, 14 (1), 33−38. (31) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS: A linear constraint solver for molecular simulations. J. Comput. Chem. 1997, 18 (12), 1463−1472. (32) Cauchy, A. Méthode générale pour la résolution des systèmes d’équations simultanées. C. R. Math. 1847, 25, 536−538. (33) Flyvbjerg, H.; Petersen, H. G. Error-estimates on averages of correlated data. J. Chem. Phys. 1989, 91 (1), 461−466. (34) Lawson, C. L.; Dutta, S.; Westbrook, J. D.; Henrick, K.; Berman, H. M. Representation of viruses in the remediated PDB archive. Acta Crystallogr., Sect. D: Biol. Crystallogr. 2008, 64, 874−882. (35) Krissinel, E.; Henrick, K. Inference of macromolecular assemblies from crystalline state. J. Mol. Biol. 2007, 372 (3), 774−797. (36) Henrick, K.; Thornton, J. M. PQS: A protein quaternary structure file server. Trends Biochem. Sci. 1998, 23 (9), 358−361. (37) Daura, X.; Gademann, K.; Jaun, B.; Seebach, D.; van Gunsteren, W. F.; Mark, A. E. Peptide folding: When simulation meets experiment. Angew. Chem., Int. Ed. 1999, 38 (1−2), 236−240. (38) Kabsch, W.; Sander, C. Dictionary of protein secondary structure−pattern-recognition of hydrogen-bonded and geometrical features. Biopolymers 1983, 22 (12), 2577−2637.

K

DOI: 10.1021/acs.jpcb.6b02698 J. Phys. Chem. B XXXX, XXX, XXX−XXX