The Role of Oligomerization and Cooperative Regulation in Protein Function: The Case of Tryptophan Synthase

The oligomerization/co-localization of protein complexes and their cooperative regulation in protein function is a key feature in many biological systems. The synergistic regulation in different subunits often enhances the functional properties of the multi-enzyme complex. The present study used molecular dynamics and Brownian dynamics simulations to study the effects of allostery, oligomerization and intermediate channeling on enhancing the protein function of tryptophan synthase (TRPS). TRPS uses a set of α/β–dimeric units to catalyze the last two steps of L-tryptophan biosynthesis, and the rate is remarkably slower in the isolated monomers. Our work shows that without their binding partner, the isolated monomers are stable and more rigid. The substrates can form fairly stable interactions with the protein in both forms when the protein reaches the final ligand–bound conformations. Our simulations also revealed that the α/β–dimeric unit stabilizes the substrate–protein conformation in the ligand binding process, which lowers the conformation transition barrier and helps the protein conformations shift from an open/inactive form to a closed/active form. Brownian dynamics simulations with a coarse-grained model illustrate how protein conformations affect substrate channeling. The results highlight the complex roles of protein oligomerization and the fine balance between rigidity and dynamics in protein function.


Introduction
The formation of protein oligomeric units often produces increased stability with improved function for the multi-enzyme complexes [1]. The co-localization of protein subunits can shape the active sites, allow allosteric cooperativity, provide an additional level of signaling or regulation, and even permit channeling of intermediates during an enzymatic turnover, which are some of the prime concerns in protein chemistry from the mechanistic point of view [2][3][4][5][6][7][8][9][10]. Such protein dynamics are long recognized to be intimately linked to enzymatic catalysis, but their relationship is exceedingly challenging to delineate [11]. Several experimental and computational studies have probed these fundamental enzymatic processes and their relationships and have provided invaluable insights into the molecular mechanisms [12][13][14][15][16][17]. Hemoglobin is one of the classical and well-studied proteins that exhibit large-scale ligand-induced conformational changes and allosteric cooperativity during the regulation of oxygen transportation. However, for a more complicated and larger system such as tryptophan synthase (TRPS), understanding protein function in relation to the protein dynamics and formation of the multienzyme complex becomes even more challenging.
The current work investigated TRPS, a pyridoxal 59-phosphate (PLP)-dependent abba protein complex that catalyzes the last 2 steps of tryptophan biosynthesis in bacteria, fungi and plants.
Research studies conducted over the past 40 years have revealed interesting structural, dynamic and mechanistic features of this protein. The protein was the first known enzyme to exhibit 2 distinct catalytic activities modulated by allosteric and synergistic interactions and demonstrating an intermolecular substrate channeling process through a 25-Å long tunnel without exposing the intermediate to the environment (see Figure 1(a)). The asubunit of TRPS resembles TIM barrel protein and is composed of 2 functionally important a-loops, L2 (a-residues 53-60) and L6 (a-residues [179][180][181][182][183][184][185][186][187][188][189][190][191][192][193], that surround the a-active site. The significant contributions of these loops in the a-catalysis and a/ b-intersubunit communications have been widely recognized by both experimental and computational work [18][19][20][21][22][23]. Within the superfamily of PLP-dependent enzymes, the b-subunit of TRPS is classified as fold type II (see definition in Text S1) [24] and contains a movable communication domain (COMM domain; bresidues 102-189). The b-H6 of the COMM domain (residues 165-181) preferentially interacts with flexible a-L2 and a-L6 and mediates intersubunit allosteric communication. Both aand bsubunits can adopt open and closed conformations. A fully closed conformation is proposed to be the active state of the protein in terms of catalysis and substrate channeling.
Although the isolated aand b-monomeric units of TRPS can independently catalyze the aand b-reactions, respectively, the rate is very slow [25][26][27]. Steady-state kinetic studies [28] revealed that the rate of the a-reaction in the isolated a-subunit is ,100 times slower than that in the abba tetramer of Escherichia coli TRPS, which has 84% identities and 94% similarities with the Salmonella typhimurium TRPS used in our simulation studies. This observation reflects a strong synergistic effect of subunits on the acatalysis in the multi-enzyme complex. However, the synergistic effects on the b-catalysis are less pronounced. The rate of the breaction in the isolated b-subunit of Zea mays TRPS (ZmTSB1), which shares 96% identity with the bacterial b-subunit of TRPS, is only 1.5 times slower than the oligomeric TRPS of Z. mays [29].
While studying the stability of TRPS, Yutani and co-workers found that the isolated aand b-subunits of Pyrococcus furiosus TRPS, which share 35% and 59% sequence identity with the aand b-subunits of the S. typhimurium TRPS, respectively, are highly stable [30][31]. The study concluded that entropic effects are the major factors contributing to the stability. Similar results have been observed for Thermus thermophilus, a hyperthermophile with 30% and 55% identical amino acid sequences to the corresponding aand b-subunits of the S. typhimurium TRPS, which indicate the importance of entropic effects in stability of the monomeric subunits [32]. Other kinetic studies investigated the homologs of the S. typhimurium a-subunit, such as BX1 (33% identical to the S. typhimurium a-subunit) and indole-3-glycerol phosphate lyase (IGL) from Z. mays. Both enzymes can efficiently catalyze the a-reaction without the other protein partner, but BX1 and IGL are about 1400 and 1150 times, respectively, more efficient than the isolated a-subunit of the E. coli TRPS [33]. The faster reaction rate for BX1 may be due to a highly stable Glu134 (structurally and functionally equivalent to the a-Glu49 of TRPS). Unlike the flexible a-Glu49 of TRPS, Glu134 of BX1 is rigid and preferably stays in the active conformation [34]. This finding suggests that efficient catalysis may require a fine balance between stability and flexibility of enzymes, although the detailed molecular aspects of such linkages are not clear.
In this study, we addressed fundamental questions of protein chemistry, including 1) the importance of oligomerization of protein subunits, 2) understanding subunit cooperativity and correlative motions, 3) the linkage between allostery and cooperativity with protein function, and 4) protein conformational changes in substrate channeling. We performed several sets of explicit water molecular dynamics (MD) simulations of a/bdimeric and isolated aand b-monomeric units of the S. typhimurium TRPS with and without ligands. Notably, the isolated aand b-monomeric units are folded proteins and are stable in solution experimentally, but their catalysis rates are reduced [34]. The trajectories were analyzed, and intra-and inter-subunit correlated motions were illustrated. The ligand-protein interaction energies, entropic effects, and H-bond networks were also studied. Brownian dynamics simulations with a coarse-grained model were performed on selected protein conformations from the MD simulations to study substrate channeling.

Materials and Methods
Construction of the ligand-free (LF) open conformation a/b-dimeric unit Since the protein data bank only contains a-subunit with a closed a-L6 loop, we performed a 15-ns MD simulation with a generalized Born (GB) implicit solvent model to obtain an open a-L6 loop conformation [35]. This method has been already employed for studying the HIV-1 protease flaps to successfully demonstrate the open and closed states of this protein [36]. The initial structural coordinates for the a-subunit were obtained from the Salmonella typhimurium TRPS (PDB entry 2J9X); the a-site ligand was manually removed [37]. The coordinates of three missing residues (Ala190, Leu191, and Pro192) in the a-L6 loop were taken from a completely closed S. typhimurium TRPS (PDB entry 3CEP) [38]. After a subsequent minimization, equilibration and MD simulations with the GB model in the Amber package [39], several open conformations of the ligand-free a-subunit were collected on the basis of the distance between a-Thr183 (a-L6) and a-Asp60 (a-L2). The open conformations of the ligandfree a-subunit were combined with a ligand-free open b-subunit (PDB entry 1QOQ) to construct several ligand-free TRPS with open aand open b-subunit [40]. The modeled a/b-dimeric units were minimized and equilibrated in explicit water. The systems were then subjected to a minimum of 13-18 ns of explicit MD simulations and important distances were subsequently analyzed. The most stable ligand-free a/b-dimeric unit in terms of smooth distance fluctuations was then selected for a 60-ns MD simulation by use of the NAMD 2.6 program [41].

Construction of the ligand-bound (open conformation) and ligand-bound-reference (closed conformation) dimeric units
A ligand-bound complex was constructed by placing both aand b-site ligands in the binding sites. IGP was docked into the asite of the ligand-free complex obtained from the procedure described in the previous section (the detail parameters of proteinligand docking are given in Text S1). Since the side-chains of the a-site produced considerable changes during the free protein simulation (in particular the a-Phe212), molecular docking programs could not reproduce the crystal structure conformation of IGP. Therefore, the substrate was manually placed into the binding site, and the distances of catalytically important residues a-Asp60 and a-Glu49 with IGP were maintained, as suggested by experiments. The b-site ligand, aminoacrylate, was docked into the b-subunit of the ligand-free a/b-complex by use of the Autodock4 package [42]. The choice of IGP and aminoacrylate as ligands for aand b-sites, respectively, ensures the closed conformation of the a/b-complex. The system containing aand b-site ligands is termed the ligand-bound (LB) complex. After subsequent minimization and equilibration, a 100-ns MD trajectory was collected to observe the possible ligand-induced conformational changes in the complex. Since the simulation may require a very long time (probably a couple hundred ns long) to

Author Summary
Conformational changes of enzymes are often related to regulating and creating an optimal environment for efficient chemistry. An increasing number of evidences also indicate that oligomerization/co-localization of proteins contributes to the efficiency of metabolic pathways. Although static structures have been available for many multi-enzyme complexes, their efficiency is also governed by the synergistic regulation between the multi-units. Our study applies molecular dynamics and Brownian dynamics simulations to the model system, the tryptophan synthase complex. The multi-enzyme complex is a bienzyme nanomachine and its catalytic activity is intimately related to the allosteric signaling and the metabolite transfer between its aand b-subunits connected by a 25-Å long channel. Our studies suggest that the binding partner is crucial for the ligand binding processes. Although the isolated monomers are stable in the ligand-free state and can form stable interaction if the substrate is in the final bound conformation, it has higher energy barrier when ligand binds to the active site. We also show that the channel does not always exist, but it may be blocked before the enzyme reaches its final bound conformation. The results highlight the importance of forming protein complexes and the cooperative changes during different states. exhibit the switching of the aand b-subunits from open/semiopen (LB state) to the completely closed states, we also run a reference simulation with a completely closed conformation. Therefore, another TRPS system, with IGP and aminoacrylate in the aand b-site, respectively, was prepared by using the initial coordinates from a crystal structure (PDB entry 3CEP). This is our reference structure with completely closed aand b-subunits, which we termed the ligand-bound-reference (LBR) complex. We created a 50-ns MD simulation after subsequent minimization and equilibration processes.

Construction of the isolated monomeric units of TRPS
The isolated monomeric units for all three states (LF, LB and LBR) were simply prepared by splitting the a/b-dimeric units into their subsequent aand b-monomers, so that the initial geometries of isolated monomeric units were exactly the same as their corresponding subunits in the dimeric unit for comparison.

Computer simulation protocol
For the molecular dynamics simulations, the ff03 amber force field and general amber force field (GAFF) were applied to both a/b-dimeric and isolated aand b-monomeric units (LF, LB and LBR TRPS complexes) [43][44]. An antechamber was used to create the topology and coordinate files for the ligands [45]. The protonation states for histidines, aspartates and glutamates were assigned by the MCCE program [46]. The TRPS dimeric units contain one aand one b-subunit, whereas isolated monomers contain only one of each subunit.
Although no substrates bound to the LF dimeric and isolated monomeric units, a PLP molecule was kept as a co-factor in the bactive site. The systems were electronically neutralized by the addition of 14 Na + ions for the a/b-dimeric units and 6 and 8 Na + ions for the isolated aand b-monomeric units. The LB TRPS represents a transition stage of the ligand binding process and was constructed by manually docking a substrate into a free subunit (see reference [47] for details). The system includes 3-indole-D-glycerol-39-phosphate (IGP) in the a-active site and aminoacrylate in the bactive site; systems were subsequently neutralized by the addition of 13, 5 and 8 Na + ions for the a/b-dimeric and isolated aand bmonomeric units, respectively. Both LF and LB complexes have one Na + ion placed close to the b-active site, as suggested by experiments. The LBR complex refers to a completely closed state of TRPS comprised of IGP and aminoacrylate in the band bactive sites of the complex, respectively. The carbonyl group of aminoacrylate was unprotonated, and six crystal waters were kept in the b-site. The Cs + ion located close to the b-active site in the crystal structure (PDB entry 3CEP) was replaced with the Na + ion; 13 more Na + ions were added to neutralize the a/b-dimeric unit; and 5 and 8 Na + ions were used to neutralize the isolated aand bmonomers, respectively.
All 9 complexes were solvated by use of a 12 Å TIP3P water box with the xleap program in the Amber10 package [39]. Each dimeric unit has about 86000 atoms, whereas isolated monomers have #48000 atoms. The initial energy minimization for water molecules involved the sander program in Amber10. The NAMD 2.6 program was then used for further minimization, equilibration and production runs. Before equilibration, the systems were gradually heated from 250 to 300 K for 30 ps. The resulting trajectories were collected every 1 ps. The total trajectory lengths for the a/b-dimeric units were 60, 100 and 50 ns for LF, LB and LBR states, respectively. For the isolated a-monomeric units, the trajectories were 50 ns long for both the LF and LBR states, and 150 ns for the LB state. The production runs for the isolated bmonomeric units were 56, 126 and 45 ns for the LF, LB and LBR states, respectively. The NPT ensemble was applied, and periodic boundary conditions were used throughout the MD simulations. A temperature of 298 K was maintained by use of a Langevin thermostat with a damping constant of 2 ps 21 , and the hybrid Nose-Hoover Langevin piston method was used to control pressure at 1 atm. The SHAKE algorithm was used to constrain the length of all bonds involving hydrogens; therefore, the time step was set to 2 fs. The non-bonded interactions were truncated at a distance of 14 Å with a switching beginning at 12 Å . The particle mesh Ewald method was used to treat long-range electrostatic interactions beyond the cut-off limit. The VMD program [48] was used for visualization and graphical representation. PyPAT script was used to analyze the H-bond network and MutInf [49] for the correlated motions in simulated trajectories. RMSF and entropy were calculated by Bio3D [50] and T-Analyst [51], respectively.
The Brownian dynamics simulation algorithm, together with a coarse-grained model (CGBD), was used to study the motions of the indole molecule in the channel formed by the aand b-subunits. The CGBD simulation method has been well described [52][53]. In our simulation, the amino acids are represented by one bead placed at the Ca of each residue [54]. Most residues were assigned an effective radius from an existing publication [55]. For residues in the active sites and along the channel, the bead radius was measured by the distance between the Ca and side-chain based on a crystal structure (PDB entry 3CEP). For indole, each ring is represented by one bead, and an effective radius was based on the size of the pyrrole and benzene ring of 1.6 Å and 1.9 Å , respectively.
The protein is held rigid, and the motion of each bead of indole is simulated with use of the BD algorithm of Ermak and McCammon [56] and Shen et al. [57]. Although the slower protein fluctuations might have a role during indole channeling, the coupling between protein conformational changes and indole motion was not taken into account in this study [58][59]. Multiple protein conformations were chosen for the CGBD simulations. The diffusion coefficient used in the algorithm to move a bead was computed by the Stokes-Einstein equation, and the viscosity of water was set to 1 cp (T = 293 K). In our coarse-grained model, the beads of indole are linked by a virtual bond, and Coulombic and van der Waals interactions were applied for intermolecular interactions [54][55][56]. A Lennard-Jones type functional form was used for van der Waals interactions, U vdw = 0.5[((r i +r j )/ (r ij )) 8 21.5((r i +r j )/(r ij )) 6 ], where r i and r j are the effective radii of beads i and j, respectively. The Coulombic interaction was approximated by U elec = q i q j /e ij r ij , and a distance-dependent dielectric coefficient (e ij = 4r ij ) was used to avoid unrealistic in vacuo Coulombic interactions [60][61].
Conformations for the simulations are snapshots taken from 0, 6, 12, 24, 30, 40 and 50-ns MD simulations in the LBR state; 2, 12, 24, 48-ns MD simulations in the LF and LB states. All the snapshots were superimposed on the crystal structure (PDB entry 3CEP) and the coarse-grained indole molecule was placed in the same position in the a-active site shown in the crystal structure. For each snapshot, 500 different random number seeds were used to study motions of indole as it approached the b-active site. The simulations used a 50-fs time step and were run for 2-4 ms. A simulation was terminated if indole reached the b-site or escaped farther than 40 Å of the a-active site. If indole cannot reach the b-site within 4 ms, then we consider that the channel is blocked. If indole diffuses farther than 40 Å of the a-active site, then we consider that indole escapes, since it is unlikely that indole diffuses back to the active sites. We computed a distance between one bead of indole and the center of mass of the b-state in a given protein conformation to determine whether the indole reacted. If the distance was closer than 5 Å , then the indole reacted at the bactive site.

Interaction energy and entropy calculation
The MM-PBSA approach was used to compute the ligandprotein interaction energies. The total energy E tot (r) can be divided into two terms: potential energy term, U(r), and solvation energy term, W(r), both functions of the coordinate r. The molecular mechanical energies were computed in a single MD step in the Sander module using a cutoff value of 40 Å for the nonbonded interactions. The solvation energy can be further decomposed into a Poisson-Boltzmann term, W PB , for electrostatic solvation free energy [62], and a cavity/surface area term, W np , for nonpolar solvation free energy [63][64]. For the electrostatic component of the solvation energy, the dielectric constant of the interior protein (solute) was set to 1, whereas an implicit solvent dielectric constant of 80 was defined for the solvent region. The nonpolar solvation free energy was approximated with the commonly used solvent-accessible surface area (SASA) model. The SASA was estimated with a 1.6-Å solvent-probe radius as implemented in Sander. Amber10 was used to compute all energy terms for each snapshot saved during the MD simulations, with waters removed [39]. The change in mean energy on molecular interactions can be split as follows: where DU c represents the changes in valance energy (bond, angle, dihedral and improper dihedral energies), DU vdw represents van der Waals interactions, DU ele represents Coulombic interactions, and DW PB and DW np represent polar and nonpolar solvation free energy, respectively. Each individual interaction energy term is calculated according to the following equations: where D,E Protein-ligand . is the ligand-protein interaction energy. Note that the valance energy term is cancelled during the calculations. The configurational entropy S consists of conformational and vibrational parts, which describe the number and width of occupied energy wells, respectively [65][66][67], computed from each dihedral angle. The configurational entropy is calculated by the Gibbs entropy formula [68]: where p(x) is the probability distribution of dihedral x and R is the gas constant. T-analyst was used to compute the Gibbs entropy, and only the internal dihedral degree of freedom of rotatable dihedrals is considered in the entropy calculations. The absolute temperature T was set to 298 K in this study. The change in configurational entropy of dihedrals of interest between a bound and free state can be obtained by TDS config. = TS bound 2TS free .

Results/Discussion
TRPS is one of the best-characterized examples of an oligomeric enzyme with stringent allosteric regulation of the catalytic reaction. The enzyme has been proposed to cycle between a low-activity open conformation in the ligand-free (LF) state and a high-activity closed conformation in the ligandbound-reference (LBR) state. The allosteric interactions are significantly influenced by the presence of aand b-site ligands. Experiments suggest that destabilizing the a/b-interface or separating the aand b-subunits loses allosteric communication, thus resulting in impaired catalysis, particularly at the a-site [22].
The 9 simulations starting from the a/b-dimeric unit or the isolated monomers with different states i.e. ligand-free (LF; IGPfree and/or aminoacrylate-free but PLP), ligand-bound (LB; IGPbound and/or aminoacrylate-bound to the semi-open conformation proteins), and ligand-bound-reference (LBR; IGP-bound and/or aminoacrylate-bound to the closed conformation proteins) allow us to investigate the cooperativity between subunits and protein allostery induced by ligand binding. Moreover, we used Brownian dynamics simulations to study the coupling between the conformational changes and substrate channeling processes.

Allosteric communications in the free and bound dimeric complex
Effective local or allosteric protein communication is a key to protein function. In most macromolecules, these communications are usually governed by non-bonded inter/intra-molecular interactions, such as van der Waals and electrostatic attractions and hydrophobic effects. Among these attraction forces, changes in hydrogen bond networks and surface areas are useful quantitative measurements for protein communication. Figure 1(b) demonstrated that interactions at the a/b-interface in TRPS combine hydrophobic interactions [69], and salt bridges and H-bonds. Experimental mutational studies for some of these interacting residues show that the salt bridges and H-bonds regulate allosteric and synergistic motions in the protein complex. A quantitative comparison of some H-bond networks, across the subunits and within the subunits, at the a/b-interface of LF and LBR dimeric units is shown in Figure 2(a, b). The analysis reveals a stronger communication at the a/b-interface of the LBR dimeric unit than at that of the LF dimeric unit. This finding suggests that binding of ligands in the aand b-active sites of TRPS enhances the subunit communications, which are necessary to synchronize the catalysis taken in both aand b-active sites located 25 Å apart from each other.
Correlated motions in proteins are ubiquitous and often related to protein functions. Assessing such correlations is therefore crucial for understanding protein function. Although we observed more inter-subunit interactions in the LBR state, the correlations are more pronounced in the LF state. The complex is also more flexible in the LF state, and the motions are not random but are in concert. Figure 3 shows a comparative correlation of regions important for subunit communication, such as a-L2, a-L6, b-H6 of COMM domain and residues at the a/b-interface of the TRPS dimeric unit in the LF and LBR states obtained by the use of the MutInf package [49]. With a few exceptions in the b-subunit, in general, the correlation is weaker at/near the dimeric interface in the LBR state than the LF state; loops a-L2 (red rectangular box) and a-L6 loops (blue rectangular box) show significant correlation in the LF state. The correlation map suggests that the a-subunit (a-L2, a-L6 and the interfacial residues) and b-H6 of the COMM domain (pink rectangular box) has weak correlation in the LBR state. A possible reason for a weaker correlation is that stronger inter-subunit interactions rigidify those regions (Tables 1 & 2) upon binding of the ligands, resulting in smaller magnitudes of correlative motions. We suggest that in the LF state, the concerted motions may help and guide the loops and the COMM domain to close when substrates bind to the active sites.

Structural flexibility of isolated monomers versus a/bdimeric unit
Fully closed protein conformations are believed to appear only when both aand b-site ligands are present in the a/b-dimeric complex; they provide the optimized geometry necessary for enzyme catalysis. The closed conformations can optimize substrate-protein interactions to stabilize the substrate in the active site. To quantify the stability of substrates binding to the a/ b-dimeric unit versus the isolated aor b-monomer, we performed end-point energy calculation, also known as MM-PBSA calculations. Although more rigorous free energy calculation methods, such as umbrella sampling or metadynamics, need to be applied to get detailed free energy profile, it may need excessively large computational power to fully sample the energy landscape for a system with this big size [70][71]. A simple thermodynamic cycle and single-trajectory post-processing allow for efficiently computing the various contributions and differences in ligand binding to the dimeric and isolated monomeric units.
Because the catalytic rates are greatly reduced in the substrateisolated monomeric complexes, we anticipated that both ligands might show weaker intermolecular attraction in the monomers. Unexpectedly, both aand b-substrates in the substrate-isolated monomeric complexes showed fairly strong intermolecular attractions in the LBR TRPS state than in other states, which suggests that the monomers are nearly as stable as the dimeric unit. However, substrates in the LB monomer have higher intermolecular energies and are unstable, and the conformations of ligand-monomer complexes deviate from their dimer conforma-  tions, especially in the a-subunit. Table 3 gives a comparison of ligand-protein interaction energies in the LB and LBR monomeric and dimeric complexes. The values of DE total for the isolated LBR monomers (a = 27.6 and b = 280.5 Kcal/mol) are similar to those in the LBR dimeric unit (a = 210.9 and b = 276.7 Kcal/ mol), which suggests that the ligand-protein intermolecular attractions do not have significant differences between the isolated LBR monomers and the dimer. The changes in the electrostatic (D,U ele +W PB .) and the non-polar solvation (DW np ) energy terms upon dissociation of the dimeric unit into the monomeric units are insignificant in the LBR states. In the LB states, the transition states during ligand binding processes, substrates interact weakly with the protein in both aand b-monomers and the a/bdimeric unit. Interestingly, the interactions are much weakened in the isolated a-monomer, which indicate that without forming an a/b-dimeric unit, ligand binding substantially disturbs the stability of the protein. Overall, both aand b-substrates are less stable in the LB state than are ligands in the LBR state. The LB state is in association processes, whereas ligands are binding to TRPS. These findings suggest that the a/b-dimeric unit helps both aand b-site ligands bind in the active sites and bring the proteins to the closed conformations through a systematically advanced allosteric communication across the a/b-interface. Absence of interface communication (i.e., the isolated monomers) detains the transition of open conformations to closed conformations and results in the deceleration of catalysis in monomer complexes.
Of interest is knowing whether the isolated aand bmonomeric units are more disordered than a/b-dimeric units, which may be less favorable for ligand binding. Therefore, we calculated the root mean square fluctuations (RMSFs) of Ca atoms and torsional entropy for each residue for the first ,50 ns long trajectory of the LF and LBR states. Figure 4 shows a comparison of the RMSF values of the isolated a-monomeric unit and the a/ b-dimeric units for both LF and LBR states, which match well with the trends of the fluctuations in the B-factors of the TRPS crystal structures. The RMSF plot clearly indicates that most of the regions in the isolated a-monomeric unit are more rigid, as compared to the a-subunit of the a/b-dimeric unit in the LF state, while an opposite trend can be observed for the LBR state. The effect of the ionic strength on the dynamics of the hydrophobic surface for the isolated a-monomeric unit in the LBR state seems negligible. The RMSF plot obtained from a 50 ns long explicit water MD simulation with 100 mM NaCl concentration is compared with those of the isolated a-monomeric unit and the a-subunit of the a/b-dimeric unit, and is given in Figure 1 in Text S1. For the b-monomeric units, in general, the difference in the RMSF values is insignificant. To quantitatively account for these flexibilities, torsional entropy was computed for the isolated monomeric and dimeric units in different states. The entropy computed for the peptide bond L angles was similar with all simulations, so we focused on other more flexible dihedral angles. The total entropic contributions for the backbone (K and J) and sidechains (SCs) indicated that in the LF state, both isolated aand b-monomeric units are surprisingly more rigid than the dimeric unit (see Tables 1 & 2). Khare et al. [72] have observed a similar behavior in the wild-type Cu, Zn superoxide dismutase (SOD1) enzyme, where some residues are more rigid in the monomeric SOD1 as compared to dimer and are coherent with the NMR data. In TRPS, the difference is particularly significant in the sidechain rotation. Regions involved in ligand binding and closing the binding sites, such as aand b-active sites, a-L6, a-L2, and b-COMM, show a pronounced decrease in sidechains motions of the isolated monomeric units, thus contributing to their rigidity. The hydrophobic binding interface between the subunits provides alternative contact points that allow sidechains of residues in the dimeric unit to adopt different binding conformations (data not shown). In addition, the correlated motions through non-direct sidechains contacts also increase protein flexibility, so such correlated motions vanish in the monomer.
Upon ligand binding, the protein flexibility was reduced largely in the dimeric unit; however, surprisingly, no significant entropic penalty was found in the isolated LBR monomers (TDS = ,1.6 and ,1.7 kcal/mol for the aand b-subunits, respectively). When the substrate binds in the a-active site, the dihedrals entropy of the a-subunit loses 57.6 kcal/mol in the dimeric unit (Table 1). Because ligand IGP has intra-molecular interactions and is not very flexible in its free state, the entropy loss from reducing the flexibility of a few rotatable bonds of IGP is not significant (,2 kcal/mol). The difference is comparatively less sizable in the b-subunit; binding the b-ligand to the active site produces a protein dihedral entropy loss of 32.5 kcal/mol in the dimeric unit ( Table 2). Interestingly and unexpectedly, without the partners, the isolated monomers are more rigid in the LF state. Although binding a chemical ligand to a protein may always result in losing  the configuration entropy of the chemical compound, binding a protein ligand to a protein partner may have more complex behavior, such as gaining some flexibility in the TRPS system. A detailed study may be required to fully characterize and understand this behavior. As the LF monomer is more rigid, when ligand IGP binds to the active site, the flexibility changes between the LF monomer and the LBR monomer are also less substantial than those in the dimeric unit. In the LBR state, comparing the total entropy calculations shows that both aand b-monomers are more flexible than the dimeric unit.
In the isolated monomeric form, after ligand binding (the LBR state), the subunit has more freedom to change its conformation slightly to minimize the entropic penalty associated with gain of enthalpy in ligand-protein binding. For example, Glu49, Asp60, Gln65 and Asp130, which interact with the a-site ligand or communicate with the b-subunit, are able to form H-bonds with  different atoms. The carboxylate or amide groups of these residues flip along with the sidechain (See Figure 2 in Text S1), which preserves the flexibility but forms multiple sets of H-bonds to gain reasonable substrate-protein interactions. Similarly, the LBR state shows a frequent flipping of carboxylate group coupled with the sidechain rotation in Glu350 and Glu172 of the b-active site in the isolated b-monomer. Figure 5 displays the percentage of H-bond networks for LBR monomeric and dimeric units for residues at the a/b-interface and active sites. Details regarding average distances and angles of H-bond are given in the Table 1 in Text S1. We found that in TRPS, generally, the loss of inter-subunit H-bond at the a/b-interface in the isolated monomeric unit is partly compensated by the formation of new H-bond networks within the subunit (see Figure 5a), as was also reported for other monomeric proteins [73]. We observed that the total number of H-bonds in the LBR monomeric states increased by at least 3-4 times as compared to the dimeric units (data not shown). Therefore, the isolated monomeric units are not less stable than the dimeric units energetically. In contrast, the dimeric unit has less room to adopt different protein conformations, which results in larger entropy loss in the LBR state. However, some residues (blue circles in Figure 6(c-d)) at the b-interface and in the b-active site show strong correlations with other interfacial residues and residues present in the b-active site of the LBR b-dimeric complex. These correlations are almost diminished in the isolated b-monomeric complex.
Overall, the effects of ligand binding and oligomerization on the 2 subunits are considerably different. The reasons may be that i) the b-subunit is larger and more rigid than the a-subunit, ii) the b-active site is buried within the subunit, right beneath the COMM domain and located relatively far from the interface as compared with the a-active site, and iii) the motion of the COMM domain in the b-active site is small as compared with the motion of loops in the a-subunit. The correlations within the b-subunit are minor as well ( Figure 6). For example, residues involved in the communication with the a-subunit, 165 to 181 in the b-H6 of the COMM domain (pink rectangular box), are correlated weakly with the b-interface (green underline) and the b-active site residues (cyan underline) in both the monomeric and dimeric units.

The power of two: Role of forming the dimeric complex
An increasing number of studies show that co-localization of proteins contributes to the efficiency of cellular signaling events and metabolic pathways [74][75]. TRPS is one of the model systems, and the dimeric unit is the minimal function structure. To mimic nature's synergy, one recent strategy is to engineer proteins that consider their spatial organization [76]. However, for enzymes such as TRPS, which are involved in regulation and synchronization in producing intermediate and final products, simply assembling multiple proteins in close proximity may not be enough. The dimeric unit forms a channel for efficient intermediate transportation, but the aand b-subunits also use the inter-subunit interactions to assist in conformational transitions and synchronize the reactions in both active sites.
Our studies suggest that without a protein partner, both of the isolated aand b-monomers form a stable and fully closed conformation when ligands are both bound in the active sites, which is the active form of the enzyme. However, the monomers, in particular the isolated a-monomer, may require an extended time to transit from an open/inactive form to a closed/active form. Forming the dimeric unit does not rigidify TRPS to form a pre-organized conformation for ligand access and to reduce entropic loss upon ligand binding. However, instead, it stabilizes the protein when the protein conformation is perturbed by the substrates during the binding processes. As a result, the dimeric unit has a smoother active-inactive transition. Notably, for both the isolated monomers and dimeric unit, the proteins sample both open and partially closed conformations, but the open (inactive) form is favored while the ligand is unbound in the LF state.
Presumably, because the hydrophobic interface provides alternative sidechain contacts and inter-subunit interactions, the dimeric unit is more flexible than the isolated monomer. Although the more flexible LF state in the dimeric TRPS results in larger configuration entropy loss upon ligand binding, we suggest that it also contributes to ligand recruitment. While a substrate is loosely bound to the binding site, the active-inactive transition rates increase, as was recently suggested by Zhou [77]. The binding sites are moving toward the fully closed conformation, and the binding mechanism gradually shifts from population shift (conformation selection) to induced fit [78][79][80]. However, as revealed by our simulations, the more unstable monomeric conformations in the ligand binding processes introduce a larger transition barrier; thus the transition rates can be decreased significantly. The dimeric unit uses the inter-subunit interactions to make the conformational transition easier.
In the LB state, the DE tot of the isolated a-monomer is ,9 kcal/mol larger than that of the a/b-dimeric unit, while the DE tot for the b-monomeric and the dimeric unit lies within the standard error (see Table 3). The value suggests that the transition rate may be decreased by several orders of magnitude in the isolated a-monomer but reduced only a little in the isolated b-monomer. The results are in good agreement with experiments showing that the catalytic rate is ,100 times slower in the isolated a-monomer but only 1.5 times slower in the isolated b-monomer as compared with the abba tetramer [28][29]. The calculation further supports our conjecture that one major role of oligomerization in TRPS is to help the ligand binding processes.
In the LBR state, the isolated monomers show frequent flipping of the carboxylate group in key catalytical residues, such as a-Glu49, but the flipping rarely occurs in the dimeric unit. Multiple sets of H-bonds are established by the flipping of a carboxylate or an amide group and sidechain rotations, so the ligand-protein interactions are not weakened. However, the fluctuations can decrease the catalytic rate in the isolated monomeric units. Our work suggests that for residues directly involved in the catalysis, rigid sidechains are preferred for optimized protein function. A similar point has been concluded for the homomeric BX1 protein, whereby the protein has a rigid Glu134, the residue having the same role as a-Glu49, to enhance the catalytic rates [34]. Coarse-grained Brownian dynamics simulation for intermediate channeling One of the unique features of the TRPS dimeric unit is substrate channeling. Conformational changes may affect the availability of the channel, and a fully closed conformation is necessary to avoid intermediate escape [81]. Protecting indole from diffusing away from TRPS is crucial for producing the final product, tryptophan, because the intermediate is relatively unstable in solution.
Considering the significance of substrate channeling and the challenge of studying the process experimentally, we carried out CGBD simulation to explore the indole channeling processes (See Figure 7).
The transportation of indole in the LBR state is smooth and rapid. Almost all, 99.6%, of indole can reach the b-active site within 4 ms, and the average travel time is 39 ns. In contrast, on the basis of 4 different LF protein conformations taken from the atomistic simulations, only ,50% of indole can reach the b-active site in the LF state. Note that we manually placed an indole to the a-active site in the LF state to simulate indole diffusion when TRPS is in open conformation. The travel time of indole towards the b-active site in the LF state is similar to that in the LBR state, but about a half of indole escapes the a-active site from the open a-loop6 (Figure 7a). In the LB state, where the protein is undergoing transition from an open to a fully closed form, indole does not always flow smoothly into the b-active site. A simulation from a snapshot taken from a 24ns MD simulation showed that no indole could reach the b-active site, because of the channel blockage (Figure 7b). Other simulations from snapshots taken from 2-to 48-ns MD simulations revealed a leak in the a/b-interface resulting in only 73% of indole successfully arriving at the b-active site.
Our work suggests that the channel also has a dynamic characteristic, which is substantially influenced by the conformational changes at the active sites. An efficient substrate channeling with a maximal success rate is only possible when both subunits are in fully closed conformation, which is in good agreement with the experiments. Since both indole and the channel are mainly non-polar, no major attraction forces that steer indole to diffuse from the a-site to the b-site were observed. Instead, indole spends longer time in positions that have larger cavities, as the molecule can freely diffuse to all direction before reaching the b-active site. The detailed channeling profile has been explored with the CGBD, and the population of indole staying in the channel formed by one conformation in the LBR state is shown in Figure 8b. The peaks correspond to large space appeared in the channel. Because our model provides fairly large space in the aactive site, indole usually needs to diffuse around the site before finding the right direction to move forward. Moreover, our trajectories show that indole may diffuse back and forth a couple of times in the channel before finally reaching the b-active site, which may be one reason that the diffusion time is an order of magnitude slower than indole diffusion in water. Our coarsegrained model keeps the protein rigid, so it cannot represent correlation between intermediate diffusion and protein conformational changes. However, as indole is a small and neutral molecule, it is unlikely to have prominent intermediate-protein correlations to accelerate the intermediate diffusion. For systems where protein motions strongly correlate with ligand channeling, a fully flexible protein system with the use of multi-bead coarse-grain models may need to be applied to more accurately capture the role of protein motions [82][83].
The significance of protein oligomerization in nature is widely recognized. TRPS is a good model system revealing the crucial role of oligomerization in assuring successful ligand binding and enhancing the rates of chemical catalysis. This study showed that the oligomerization of the aand b-subunits not only provides a direct channel for efficient intermediate transportation but also permits allosteric cooperativity via inter-subunit communications to assist with conformational transitions necessary to synchronize the reactions in both aand b-active sites.

Supporting Information
Text S1 Supporting information for ''The Role of Oligomeri-