Full Enyzme Complex Simulation: Interactions in Human Pyruvate

Subscriber access provided by READING UNIV

Article

Full Enyzme Complex Simulation: Interactions in Human Pyruvate Dehydrogenase Complex Samira Hezaveh, An-Ping Zeng, and Uwe Jandt J. Chem. Inf. Model., Just Accepted Manuscript • DOI: 10.1021/acs.jcim.7b00557 • Publication Date (Web): 03 Jan 2018 Downloaded from http://pubs.acs.org on January 4, 2018

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 35 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

Full Enzyme Complex Simulation: Interactions in Human Pyruvate Dehydrogenase Complex Samira Hezaveh1, An-Ping Zeng1, and Uwe Jandt1* Institute of Bioprocess and Biosystem Engineering, Hamburg University of Technology, Denickestr. 15, 21071 Hamburg, Germany

1

Institute of Bioprocess and Biosystem Engineering

Hamburg University of Technology Denickestr. 15, 21071 Hamburg, Germany

*Corresponding author: Dr. Uwe Jandt Fax: +49 (0)40 42878 2909, Tel: +49 (0)40 42878 2847

Authors information: Samira Hezaveh: [email protected] An-Ping Zeng: [email protected] Uwe Jandt: [email protected]

1 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 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 2 of 35

Abstract The pyruvate dehydrogenase complex (PDC) is a large macromolecular machine consisting of dozens of interacting enzymes that are connected and regulated by highly flexible domains, also called swinging arms. The overall structure and function of these domains and how they organize the complex function is rarely elucidated yet. This lack of structural and dynamic understanding is frequently observed in multi domain enzymatic complexes. Here we present the first full and dynamic structural model of full human PDC (hPDC), including the linking arms binding to the surrounding E1 and E3 enzymes via their binding domains with variable stoichiometries. All the linking domains were modelled at atomistic and coarse-grained level, and the latter was parameterized to reproduce the same properties of those from atomistic model. Radius of gyration of the wild type full complex and functional trimeric subunits were in agreement with available experimental data. Furthermore, the E1 and E3 population effect was studied on the overall structure of the full complex. The results indicated that, by decreasing the number of E1s, the flexibility of the now non-occupied arms increases. Furthermore, their flexibility depends on the presence of other E1 and E3 in the proximity, even if they are associated to other arms. As one consequence, the radius of gyration reduces with decreasing number of E1. This effect provides also an indication on an optimal configuration of E1 and E3, based on the assumption that a certain stability of the enymatic cloud is necessary to avoid free metabolic diffusion of intermediates (metabolic channeling). Our approach and results open a window for future enzyme engineering in a more effective way by evaluating the effect of different linker arms lengths, flexibility and combination of mutations on the activity of other complex enzymes that involve flexible domains, including for example processive enzymes.


Page 3 of 35 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


Introduction Structure determination and model build-up for multienzyme complexes are challenging due to the dynamic and flexible interactions between the different subunits that introduce significant complexity. One of these complex enzymes with the main function in coupled metabolic channelling is the human pyruvate dehydrogenase complex (hPDC) which, in its basic form, is assembled from four subunit components: E1: pyruvate decarboxylase, E2: dihydrolipoamide acetyltransferase, E3: dihydrolipoamide dehydrogenase and E3-BP (E3 binding protein, formerly also referred to as protein X).1 E2 and E3BP form the inner huge 60-meric core of the hPDC with a dodecahedric structural shape2 and with two distinct populations of trimeric building blocks suggested recently from two variants of a substitutional model.3,4,5 Our previous results from coarse-grained simulation on the same matter showed a higher stability for the core model 48E2+12E3BP.

6

Based on this model,

two types of trimers are assumed to form, 3E2 and 2E2+1E3BP. E1 and E3 enzymes bind to the E2 and E3BP via binding domains (BD), respectively. The E2/E3BP core bridges E1 and E3 via its lipoyl domains known as swinging arms by visiting the active sites of E1, E2 and E3, respectively.7,8 With its unique properties, including self-assembly, modularity and metabolic channelling mechanism, this type of enzyme complexes can be a high value target for engineering for novel multi enzymatic reaction cascades.9,10,11 In this study, using the previously developed model of hPDC core,12,6 we built the whole hPDC complex including the linking arms binding to the E1 and E3 enzymes via their binding domains at atomistic to coarse-grained (CG) level. This allows for the interaction and dynamics study of the three enzymes including the role and dynamic of the arms which will further help to design different enzymes with variable arms lengths and flexibilities and also to create the foundation for the synthesis of artificial complexes in the future. So far, to the best of our knowledge, there are neither simulation studies available on the enzymatic 3 ACS Paragon Plus Environment


Page 4 of 35

interactions within the hPDC or any other PDC enzyme nor the arm dynamics and features, mainly because the structure based modelling of multi enzymatic systems have so far been restricted by their huge size.10 There are only few experimental studies done on other type of PDC to study the effect of linker’s length and sequence on PDC activity.13 It was previously shown by simulation data that the C-terminal truncated core exhibited an obvious disintegration of the complex 6, which was determined by an increase in radius of gyration. These results are in agreement with experimental data in which the double truncated versions of inner core in the hPDC showed better activity compared to other truncation variants. 23 In this work, we studied enzymatic interactions within the hPDC by Molecular Dynamics (MD) simulations at atomistic and coarse-grained (CG) level based on MARTINI model.14 The results showed a clear dependency of radius of gyration to the number of peripheral enzymes which decreases strongly by reducing the number of E1s and is hardly dependent on the number of E3. Also the stability of individual arms depends on the presence of E1 in the vincinity, even if they are not bound to the arm of interest. Overall this leads to an configuration of approx. 25-30 E1 and 8-12 E3 copies where the assumed “metabolic channeling” effect is optimally preserved by reducing the free diffusion path lengths from the catalytic core through the surrounding “cloud” of E1 and E3 enzymes. The investigations continue in this direction in order to understand how the number of E1 might affect the complex activity. Since similar flexible structures are also generally essential for the function of other types of multi enzymatic complexes, e.g. processive enzymes and cellulases,15 therefore investigations are also ongoing in this direction. This paper is organized as follows: The details of the modelling and force field used at coarsegrained level are given in the Material and Methods section. In the Results and Discussion section, we present the monomeric and trimeric unit simulations of the complex and 4 ACS Paragon Plus Environment

Page 5 of 35 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


parameterization of the arms are based on comparison of atomistic and coarse-grained simulation data. Then the simulation of full hPDC and variable E1 and E3 stoichiometries is shown. At this level, different proportion of E1s and E3s and their effect on the complex structure and radius of gyration will be analyzes and discussed.

Methods Atomistic Model The full models of monomers and trimers of hPDC were built up based on the previously developed models of the E2 and E3BP sub-units12 and then the arms including the binding domains were added. Since their structure is not known, the were positioned in a random extended configuration, with the only geometric condition to allow for positioning of all 60 arms plus 48 times E1 and 12 times E3 in a fully equipped 60-mer core. The E1 and E3 were then positioned at their related binding domains (non covalent binding). For the E3 bound to its binding domain there was already a PDB (1EBD) available, but not for E1. Therefore we used the autodock function in Chimera16 to detect the possible binding zone. All structure and simulation data are available upon request. Coarse-grained Model The coarse-grained (CG) models used for MD simulations of proteins are based on the MARTINI force field.14 The model of hPDC core subunit has been built based on the initial atomistic and coarse-grained models.12,6 An elastic network was applied to MARTINI model in order to maintain the stability of the protein secondary structure.17 Simulation setup



Page 6 of 35

The simulations and analysis were performed using GROMACS software package version 5.0.2.18,19 The Visual Molecular Dynamics (VMD) program was used for the visualizations.20 In accordance to the corresponding experimental works 23, all simulations were performed at neutral pH. For coarse-grained simulations using MARTINI single point water (W), a cut-off of 12 Å was applied for Lennard-Jones (LJ) and Coloumbic interactions. For the system setups with polarisable water (PW) electrostatics were treated via Particle mesh Ewald (PME) technique. The LJ potential was smoothly shifted to zero between 0.9 and 1.2 nm, and the Coulomb potential was smoothly shifted to zero between 0.0 and 1.2 nm. The temperature and pressure were maintained to the reference values (for the pressure, P0=1 bar) using the vrescale and Parrinello-Rahman with the coupling time constant τT =1.0 ps for temperature and τP=12 ps for the pressure. A time step of 5 fs was used for equilibration and was increased to 20 fs for the simulations. Equlibration was conducted until the density in the simulation box remained stable; in case of trimer simulation the equilibration time was 100 ps, while for full PDC, 50 ps was sufficient. All simulations were positioned in triclinic (trimers) or cubic (full PDC) boxes, with periodic boundary conditions applied. A minimum distance of proteins to their periodic image of at least 2.3 nm (trimers, atomistic) and 8 nm (full PDC, CG) was maintained throughout all simulations; for full PDC, this led to cubic box side lengths of ~81.7 nm. Arguably due to the very filamentous structure of the simulated enzyme complexes, single point MARTINI water beads frequently froze during equilibration when using standard simulation setups, i.e. formed crystal-like structures. This effect was completely avoided by introducing a minor random shuffling of the initial of water bead positions by ±0.5Å, followed by another energy minimization step before equilibration (50 or 100 ps) and finally simulation. All simulations were performed at 310 K. All the errors on the calculated properties have been evaluated using the block averaging method.21 Trimers 6 ACS Paragon Plus Environment

Page 7 of 35 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


After preliminary set-up of monomeric systems, trimers were constructed. We set up the two types of trimer systems, once for 3E2 and another once for 2E2+1E3BP, with and without binding of the E1 and E3. The systems were solvated both in W and PW with ~200000 water beads, respectively, in a box of ~30 nm/side. After energy minimization the MD runs were performed at NPT ensemble for ~100-150 ns. The systems were neutralized adding Clcounterions. Same setups were done also at atomistic level using GROMOS united atom (UA) and OPLS all atom (AA) force fields as a reference for the arm parameterization. Whole hPDC complex The full hPDC complex was simulated for wild-type (WT) and C-terminal truncated core versions with different numbers of E1 (ranging from 0 to 48) and E3 (ranging from 0 to 12). The systems were positioned in a cubic box of ~82 nm/side. Equilibration was performed for 50 ps for each run. Simulations were performed for at least 24 ns (total 106 simulations); due to high computational costs only few were extended to 100 ns (0xE1+0xE3; and 0xE1, 16xE1, 28xE1, 48xE1 with each 12xE3).

Results and Discussion Stability of additional subunits Simulation stability of catalytic domains of the core subunits (E2 and BP), for example based on root mean square fluctuation (RMSF), has been evaluated and discussed before12. For the connection of E3-BD to E3, a comparison of RMSF of simulated backbone atoms to RMSF calculated from the most similar available crystal structure B-factor data (PDB code 1EBD from G. stearothermophilus) has been performed to assess agreement with respect to crystal data22, as far as possible. The result is depicted in Supplementary Figure S1 in supplementary information (SI). It shows that, within the subsequently used time frame from 10 to 24 ns, the 7 ACS Paragon Plus Environment


Page 8 of 35

RMSF is well comparable to the crystallographic data. The average RMSF of atomistic simulation of (0.12 ± 0.047) nm is 0.93 times the calculated average RMSF of 1EBD. Prominent increases of simulated RMSF (especially around Asp-133, Ser-262, Gly-344) correlate to flexible domains exposed to the solvent. Likewise, the RMSF of E1, averaged over all backbone atoms within the same time range, is (0.12 ± 0.056) nm. Trimers without E1/E3 The hPDC 48E2+12E3BP model consists of two types of trimers: 3E2 (8 copies) and 2E2+1E3BP (12 copies). We built these two types of trimers based on the previous atomistic model6 and integrated the arms to the catalytic domains of E2 and E3BP. 3E2 trimer, without E1, is shown in Figure 1. E2 has a double lipoyl domain while E3BP has only one. The two types of arms each contain one binding domain (BD); a E1 specific BD at the E2 arm (visualized in Figure 1), and a E3 specific BD at the E3BP arm (not shown).

Figure 1. Visualization of a 3E2 trimer from (a) side view and (b) top view with linker arms. E1 and E3 are not attached. Monomers are identical and were colored differently for better clarity.


Page 9 of 35 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


In order to validate and assess the stability of the E2 and E3BP trimers without E1 and E3 attached to the BDs, we analysed their radius of gyration (Rg) as shown in Figure 2 and in comparison to the atomistic model. The CG model turned out to be very flexible, or unstable, in comparison to the atomistic ones, especially when not applying the elastic network (ElNeDyn). Therefore, in order to reproduce better the properties of atomistic simulations, we optimized the constraint parameter for backbone angles at CG level based on comparing the Rg of the total trimers including the linker arms with their atomistic counterparts. The finally chosen constraint parameter for backbone angles was 100 kJ/mol instead of the standard value of 40 kJ/mol. The modified CG model follows the atomistic trend much better than the standard one. The atomistic Rg value of trimer was about 5.75 nm at 100 ns which was very close to the result of 5.97 nm from dynamic light scattering (DLS) data extracted from the hydrodynamic radius of 7.699 nm.23 The results of other 3E2 trimers (not shown here) were similar to these 2E2+1E3BP ones. To check the effect of different electrostatic treatments in the system, the same simulations were also performed in PW. However no remarkable difference was observed. We also briefly compared the same atomistic model with simulations using OPLS force-field rather than GROMOS to see whether presence of arms can show any drastic difference in the dynamic of the system. Both force-fields gave similar results (Figure S2), therefore, further investigations focused on GROMOS force field only.



Page 10 of 35

Figure 2. Radius of gyration of 2E2+1E3BP trimer with lipoyl arms but no E1 and E3 attached; atomistic versus CG. The ElNeDyn constraint parameters for backbone angles are 40 kJ/mol (“CG ElNeDyn”) and 100 kJ/mol (“CG ElNeDyn modified”). W represents the standard MARTINI water and PW depicts polarizable water.

Trimers with E1 and E3 The same simulations as before were repeated, but this time including E1 and E3 connecting to E2 and E3BP binding domains (BD), respectively (Figure 3).


Page 11 of 35 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


Figure 3. Visualization of 2E2+E3BP trimer. E1s and E3 are bound to their binding domains (BD).

Results of radius of gyration at CG level compared to those from atomistic ones are shown in Figure 4. The best match belongs to the CG ElNeDyn with modified parameters when using PW. To briefly investigate potential influence of high ionic concentrations, we also performed simulation in 50 mM ionic concentration by equally adding Na+ and Cl- ions and the results are still very similar to those from PW with no increased ionic concentration. This is while the system in the MARTINI standard water shows a drastic difference and deviate strongly from the atomistic trend. Therefore for the further simulation of the complete complex we chose CG ElNeDyn in PW and in PW with increased ionic strength.



Page 12 of 35

Figure 4. Radius of gyration of 2E2+1E3BP trimer with lipoyl arms, and E1 and E3 attached; atomistic versus CG. The constraint parameter for backbone angles is 100 kJ/mol.

Previously published experimental data showed that trimeric individuals of hPDC after breakdown of the catalytic core to its subunits did not show any activity while the agglomerates of the same subunits showed highest activity.23 It was speculated that this may be related to the overall positioning of the E1 and E3 enzymes with respect to the catalytic core subunits, i.e. unphysiological positioning might lead to unefficient metabolic channeling or metabolic intermediate transfer via the flexible arms. To investigate this putative effect, the average distance between the peripheral enzymes of E1 and E3 were evaluated for individual trimers of the 3E2s and 2E2+1E3BP, as shown in Figure 5. In overall E1s showed a tendency to get apart from each other in both types of trimers. The average minimum distance between them was calculated as 11 to 15 nm at the end of simulation. This is while E1 and E3 distance reduced and in average it turned out to be about 8 to 11 nm. From our simulation data this can be due to an inefficient access of lipoyl domains to the active sites as the distance between 12 ACS Paragon Plus Environment

Page 13 of 35 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


them increased. Further it is likely that free diffusion of metabolites away from the catalytic core and through the E1 and E3 cloud can lead to a loss of volatile intermediates and unefficient metabolic channeling. In contrast, in a full core or agglomerates of trimeric subunits, the symmetric geometry and presence of side enzymes is supposed to prevent the occurrence of this, leading to active enzyme complexes.

Figure 5. Minimum distance vs. time for two types of trimers (a) 2E2+E3BP and (b) 3E2. In 3E2 trimers the sets are referring to pair distances of E1s.

Whole hPDC complex WT hPDC



Page 14 of 35

The whole hPDC structure was built after parameterization of the linker arms and binding of E1 and E3 to their respective binding domains of E2 and E3. The system was initially built for the full possible numbers of 48 E1s and 12 E3s; variants with reduced numbers of both enzymes were derived from that. Figure 6 shows the stepwise modeling visualization of the system starting from the full core of E2+E3BP, then adding the linker arms and binding domains and at the end attaching of the other peripheral E1 and E3 enzymes.

Figure 6. Full hPDC modelling starting from full core of 48E2+12E3BP with E3BPs in red and E2s in green (a) then adding arms of E2 in orange and E3BP in gray (b), E3s in purple (c) and eventually E1s in yellow (d).

For practical reasons and also because of the fact that isolated complexes after full equilibration are presumably not realistic, because they exhibit strong agglomeration as


Page 15 of 35 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


well,23 we instead focused on the dynamics of the complexes and the relative deviation between them when different numbers of E1 and E3 are involved. The radius of gyration for the complex variant including arms, but without E1 and E3, can be roughly compared with experimental data, which showed an average radius of hydration RH=26.5 nm for the single species of the full core23. For globular proteins, this would correspond to a radius of gyration of approx. Rg~20.5 nm. However, for such rotationally symmetric, but extremely filamentous proteins, it can be assumed that the ratio RH to Rg is larger, i.e. experimental Rg is smaller. The simulation of the variant without E1 and E3 attached reached Rg =20.5 nm at a time of of approx. t=22-24 ns, so we considered this time frame of the simulations to be most realistic in comparison with the single species population under experimental conditions. The arms in the simulations without E1/E3 however collapse further, presumably completely on the core surface, which is considered not realistic. This is most probably because no interactions between different PDC (cores) could be simulated due the restrictions in simulation system size. It is interesting that the system was able to accommodate all possible number of E1s and E3s. This was not expected, since from experimental data, E3 has been reported to be present in 6 to 12 copies while number of E1s detected is about only 30.24,25 Since the number of E1s seems to be important in the structure and possibly the activity of the system, we prepared series of simulations with different E1 and E3 combinations. Each configuration was built on the initial, fully occupied setup (48xE1, 12xE3), then E1 and/or E3 were randomly removed until the desired total number of E1 and/or E3 was reached. Therefore, the spatial distribution of the remaining E1 and E3 as well as the position of unoccupied places between was not biased by any constraint. Potential effects arising from any non-homogenous distribution are expected to be averaged by the relatively large number of total conducted independent simulations (106 total). 15 ACS Paragon Plus Environment


Page 16 of 35

The radius of gyration was calculated for the same at the same time point of t=24 ns, which is depicted in Figure 7; the Rg trends vs. time are shown in Figure S3. The results for these short simulation times provides no reliable information on equilibrated state but allows for relative comparison of the dynamics of the different variants. The radius of gyration depends mostly on different E1 numbers rather than E3s. The most compact system occurs when we reduced E1s to 0. Visually, this reduction seems to be the result of arms bending because of more space available on the surface of the core. The presence of peripheral enzymes in general is important for the dynamic and behavior of arms on the surface of the core. To which extend the presence proximate E1 and E3 influence not only the stability of the whole population of arms but also the arms individually is investigated in the following.

Figure 7. Radius of gyration of WT full hPDC for different number of E1 and E3 up to 24 ns.

In order to assess the influence of the presence of E1 and/or E3 in the proximity of the flexible arms on the individual dynamics of the latter, even without necessarily being attached to them, we calculated a collapse ratio of each individual arm. This was straightforward 16 ACS Paragon Plus Environment

Page 17 of 35 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


defined as the relative Rg (R’), computed by Rg at t=24 ns (the maximum time achieved by all simulations covering different stoichiometries) divided by the corresponding Rg at t=10 ns (the time point where all separate simulations start). This R’ was plotted versus the smallest distance of the same arm to the closest E1 and/or E3 (D0) at the beginning of the simulation (Figure 8 a). This yields a total of 1248 data points, exhibiting at least three clearly distinguishable populations, which are further analyzed in Fig 8 b: Population I), with D0 6nm: These show the least stability (R’[D0>6nm]=0.80±0.086, n=173). Due to the large number of samples, all of these populations can be distinguished with very high significance of p

Full Enyzme Complex Simulation: Interactions in Human Pyruvate

Recommend Documents