Allosteric regulation of β-reaction stage I in tryptophan synthase upon the α-ligand binding

Tryptophan synthase (TRPS) catalyzes the last two steps of tryptophane (Trp) biosynthesis and is a bifunctional enzyme consisting of the α- and β-subunits. Each subunit has its own ligand and the binding and releasing of the ligands greatly affect the enzyme reactions at the α- and β-active sites.^1,2 To understand the complete catalytic reaction mechanism of TRPS, it is necessary to elucidate the inter-subunit communications in different ligand-bound states. In other words, how the binding affinity and the reaction activity at an active site are affected by those at another site is a key question to be answered. Such mechanisms known as the allosteric communication and regulation are of critical importance in many enzymes and proteins.^3–6 Allostery in TRPS has been investigated both functionally and structurally.^1,2,7–29 Recently, TRPS has attracted much attention as a target for drugs aimed at inhibiting the tryptophan synthesis in the pathogenic bacteria, Mycobacterium tuberculosis (M.tb), in order to treat tuberculosis (TB).^30–32 Effective allosteric inhibitors of TRPS have been reported to kill M.tb.^33,34 An atomistic description of the allosteric regulation is essential for the design of efficient inhibitors.

TRPS from Salmonella typhimurium has been extensively studied for decades.^19–21 The α- and β-subunits in TRPS have distinct ligands at the active sites, which are separated by ∼35 Å. A flexible loop 6 (αL6) in the α-subunit and a communication (COMM) domain in the β-subunit shown in Fig. 1(a) not only serve as gates for the α- and β-active sites but also mediate the interaction between subunits.^22,23 The complete catalytic cycle of Trp synthesis by TRPS is illustrated in Fig. S1 of the supplementary material. At the α-active site, a C–C bond in the α-ligand, 3-indole-D-glycerol-3′-phosphate (IGP), is enzymatically cleaved to indole and glyceraldehyde-3-phosphate (G3P). Indole is transported to the β-active site through an internal tunnel, and the β-active site catalyzes the reaction of the indole and L-serine (L-Ser) with the pyridoxal phosphate (PLP) co-factor to produce an L-Trp. The β-reaction consists of two stages, stages I and II (see Figs. S1 and S2). Initially, PLP covalently binds to βK87 [denoted E(Ain)] and reacts with the incoming L-Ser to form an external aldimine [denoted E(Aex₁)] via a gem diamine intermediate [E(GD₁)]. Thereafter, three further processes occur in order, as shown in Fig. 1(b): (1) E(Aex₁) is converted to the quinonoid intermediate [denoted E(Q₁*)] by a proton transfer from the ligand to the βK87, (2) the dihedral angle of the hydroxy group of E(Q₁*) rotates to form E(Q₁); (3) E(Q₁) is converted to the α-aminoacrylate intermediate [denoted E(A-A)] by a proton transfer from the βK87 to the ligand and the dissociation of a water molecule. Structural and kinetic studies with site-directed mutagenesis have reported that the conformational changes and the ligand binding in each subunit strongly affect the enzyme activity.^7–18 Spectroscopic studies with Stopped-Flow Kinetic UV/Vis reported that the binding of IGP (or its analog) at the α-active site increases the reaction activity of the β-reaction stage I by about ∼3 to 10-fold.^7,24,25 Inter-subunit communication between the α- and β-active sites plays an important role in the enzymatic functions of TRPS.

Despite many biochemical experiments showing evidence of the allosteric communication in TRPS, many important questions remain to be answered. For instance, the mechanism by which IGP binding at the α-active site can enhance the reaction activity of the β-reaction stage I at the distal β-active site is still unknown. Hereafter, we use the notation of A|B, where A and B indicate the states of the α- and β-active sites, respectively. The x-ray crystal structures were solved for the wild-type (WT) TRPS in the Apo|Serine pyridoxal phosphate Schiff base (PLS) [Protein Data Bank Identifier (PDBID): 1KFJ¹³] and F9F|PLS (PDBID: 2CLL²⁴) and a mutant (βK87T) in the Apo| E(Ain). L-Ser (PDBID: 1UBS³⁵), where F9F is an analog of IGP. Although these structures provide detailed information about the ligand-bound states, it is difficult to answer the above question by simply comparing the static structures.

Molecular dynamics (MD) simulations have been used to examine dynamic aspects of proteins and enzymes, such as folding/unfolding, conformational changes, and ligand binding. Many MD simulation packages are now available.^36–39 We have also developed an MD program, GENESIS (https://www.r-ccs.riken.jp/labs/cbrt/).^40,41 Currently, MD simulations on specialized supercomputers, Anton/Anton2/Anton3⁴² and computational algorithms allow hundreds of micro-seconds, or even milli-second simulations. Enhanced sampling methods enable the simulation of binding events in enzymes with accurate free-energy landscapes of biomolecules.^43–45 Replica-exchange molecular dynamics (REMD),^46–48 generalized replica exchange with solute tempering (gREST),^49,50 Gaussian accelerated MD (GaMD),⁵¹ and metadynamics⁵² are representative algorithms for predicting the thermodynamic properties of proteins or other biomolecules. However, classical MD simulations using molecular mechanics force fields^53,54 cannot describe the chemical reactions that involve bond breaking and formation. The hybrid quantum mechanics/molecular mechanics (QM/MM) method proposed by Warshal, Levitt, and Karplus more than 40 years ago, could overcome these problems.^55–57 In GENESIS, QM/MM calculations are available by combining with highly parallelized QM software, QSimulate-QM (https://qsimulate.com/), as well as other QM software, Gaussian,⁵⁸ Q-Chem,⁵⁹ TeraChem,⁶⁰ and DFTB+.⁶¹ 

Both classical MD simulations and QM/MM calculations for TRPS have been reported in the literature.^62–71 Maria-Solano et al. computed the potential of mean force (PMF) along the open-to-closed conformational changes of the COMM domain in the α- and β-complex and a stand-alone β-subunit with mutations to find the enzyme pathways which are related to the stand-alone activity of TRPS using shortest path map (SPM) tools.⁶⁷ They also succeeded in identifying which specific amino acid substitution should be introduced for increasing the stand-alone activity of TRPS along the enzyme pathways obtained by SPM tools using ancestral sequence reconstruction.⁶⁸ Recently, we examined inter-subunit communication with MD simulations and umbrella sampling (US).⁷⁰ The simulation revealed that the IGP binding at the α-subunit induces a wide opening of the COMM domain due to a hydrogen bond switching around the interface region of the α- and β-subunits, making it easier for L-Ser to enter into the β-active site (see Fig. S1). Teixeira et al. calculated the energy profiles and conformations of the reactant/transition state/product of the whole reactions in the L-Trp synthesis by performing internal reaction coordinate (IRC) calculations and vibrational frequency calculations.⁷¹ However, the mechanisms underlying the increase in reaction activity in the β-reaction stage I after the IGP-binding at the α-active sites remain unclear. In the aforementioned study, the energy profile was calculated using the crystal structure in the IGP|E(Aex₁) state. Thus, the energy profiles in the β-reaction stage I before the IGP-binding remain unknown. Here, we studied three reaction steps in the β-reaction stage I, E(Aex₁) → E(Q₁*) → E(Q₁) → E(A-A) shown in Fig. 1(b) (see Fig. S2 for the complete β-reaction). The minimum-energy path (MEP) was obtained for the Apo|E(Aex₁) structure using the QM/MM string method.⁷² Then, the potential of mean force (PMF) along the pathways was computed using the QM/MM umbrella sampling (US) method.⁷³ The MEP search and US simulations were carried out using density functional theory (DFT) with B3LYP-D3/cc-pVDZ and aug-cc-pVDZ levels of theory, respectively, for the QM calculations. Our simulations suggest that hydrogen bond switching from βD305–E(Q₁*) to βD305-βR141 is essential in reaction step 2 of the β-reaction stage I, suggesting a mechanism underlying the increase in β-reaction activity upon IGP-binding at the α-active sites. Finally, we discuss the relationship between computational results with known TRPS crystal structures and extensive MD simulation sampling.

The QM/MM method treats the active site, where the chemical reaction takes place, using quantum chemical (QM) calculations, and the surrounding protein and bulk environment using the MM force field. The QM calculation solves the electronic Schrödinger equation given in atomic units as

where i, a, and m are the indices of the electrons, nucleus, and MM atoms, respectively. Z_a and q_m are the charge of nucleus and MM atoms, respectively, and r_XY denotes the distance between X and Y. By noting that the electronic energy $Ee$ is a function of the coordinates of the nucleus and MM atoms (R_a and R_m, respectively), the QM energy is written as

Then, the total QM/MM energy is given by

where $VLJQM/MM$ is the Lennard-Jones interaction between the QM and MM atoms, and V^MM is the force field. We refer to previous work⁷⁴ for further details on the theory. Various functions based on the QM/MM energy are available in GENESIS, such as minimization, MD and REMD simulations, vibrational analyses, and reaction path search.

The QM calculations were carried out using external QM programs. We have implemented a general interface,⁷⁵ in which GENESIS creates the input file, runs the QM program, and reads the necessary information from the output file at runtime. Various QM programs are available for this scheme. Recently, an interface was developed for QSimulate-QM, which combines the two programs through dynamic libraries.⁷⁶ The scheme achieves a high-performance QM/MM-MD simulation because it passes the information on memory without the file I/O, and also because the Message Passing Interface (MPI) processes are used seamlessly in each program. GENESIS/QSimulate-QM can perform a few nanoseconds per day with the density functional tight binding (DFTB) method, and tens of picoseconds per day using the density functional theory (DFT) method.

Umbrella sampling (US)⁷³ method is one of the importance sampling methods that use a harmonic biasing potential,

where $ξq$ is the collective variable (CV), k is the force constant, and ξ₀ is the midpoint value of the potential that are user–adjustable parameters. The umbrella potential V(q) is added to the potential energy E₀, which is obtained by the Molecular Dynamics (MD) simulation, and the total potential energy is defined by

The umbrella potential allows the system to explore the region around $ξq=ξ0$⁠, which is difficult to reach by the conventional MD simulation to overcome the high-energy barrier around the ξ₀. In the US simulation, M independent simulations are performed, which have different umbrella potentials, $Vmm=1,2,…M$ with anchor positions, $ξmm=1,2,…M$ and force constants $kmm=1,2,…M$⁠. The m dependent umbrella potential, $Vmq$ is defined as

In practice, the replica-dependent k_m values are often set to be identical. The total potential energy is given by

The multistate Bennett acceptance ratio (MBAR) method⁷⁷ computes the physical quantities of the unbiased system, which is based on the potential energy E₀. After M independent MD simulations are carried out, MBAR solves the following equation:

where N_i is the number of snapshots in the i-th window, and q_in is a coordinate from the n-th snapshot and i-th window. The dimensionless Helmholtz free energy $f̂i$ of the original unbiased system can be obtained by solving Eq. (9) self-consistently. The PMF, $Fξ$⁠, which is the Helmholtz free energy as a function of the CV of the unbiased system, is given by

where $pξ$ is the probability distribution and Z is the partition function.

An all-atom model of the Apo|E(Aex₁) bound state was prepared based on the x-ray crystal structure (PDBID: 1KFJ¹³). The ligand PLS was manually replaced with the E(Aex₁) in the β-active site and the conformation is Open/Open which means that αL6 is not obtained in the x-ray crystal structure and the salt-bridge between βR141 and βD305 is not formed. The atomic structures of missing residues were modeled by GalaxyFill⁷⁸ CHARMM-GUI⁷⁹ and a proton donor in the β-reaction, βK87, was deprotonated. The information of all-atom models is summarized in Table I. Topology files needed for MD simulations were generated by the antechamber and tleap module of AMBER18.³⁶ The force field was set to FF14SB⁸⁰ for TRPS and Generalized AMBER Force Field 2 (GAFF2)^81,82 for the ligands. Restrained electrostatic potential (RESP) charge of the ligand was obtained at the B3LYP/6-31+G(d, p) level using GAUSSIAN16.⁵⁸ TRPS was solvated in a water box of size 120 × 120 × 120 ų. A TIP3P water model was used as the force field. The system included 150 mM of K⁺ and Cl⁻ ions as shown in Fig. S3(a). All molecular graphic images were generated using Visual Molecular Dynamics (VMD)⁸³ and PyMOL.⁸⁴ 

The system was first energy minimized, and then equilibrated using MD simulations. The equilibration steps are presented in Table S1. Classical MD simulations were performed in the NVT and NPT ensembles for 1.2 ns in the periodic boundary condition (PBC). The temperature and pressure were controlled at 300 K and 1 atm, respectively, using the Bussi thermostat.⁸⁵ The time step was set to 2.5 fs for the simulations. The reference system propagator algorithm (RESPA)⁸⁶ was used as the integrator. Long-range electrostatic interactions were calculated using the smooth Particle-Mesh Ewald (PME) method^87,88 with a cut-off distance of 8.0 Å. All bonds involving hydrogen atoms were constrained using the SHAKE/RATTLE⁸⁹ and SETTLE⁹⁰ algorithms. GENESIS 1.7.1 program^91,92 was used in all MD simulations.

After the equilibration in PBC, the system was truncated in two water spheres with a 44.0 Å radius centered on the Cγ atom of αD60 and C1 of the β-ligand as illustrated in Fig. S3. An additional equilibration was performed in the NVT ensemble for 1.1 ns in non-PBC. The equilibration steps are summarized in Table S2. The temperature was controlled at 300 K using the Bussi thermostat.⁸⁵ The time step was set to 2.0 fs for the simulations with the velocity Verlet (VVER) integrator. The switching and cut-off distance was set to 16.0 and 18.0 Å, respectively.

The reaction pathways and the energy profiles of the three steps of the β-reaction stage I shown in Fig. 1(b) were obtained by the string method.^93,94 The atoms within 9.0 Å of E(Aex₁)/E(Q₁)/E(Q₁*) were relaxed, whereas others were fixed. The QM region was set as E(Aex₁)/E(Q₁)/E(Q₁*), and the side chain of the βK87 consisted of 49 atoms. QM calculations were performed at the B3LYP-D3/cc-pVDZ level.^95–98 32, 16, and 32 images were used for step 1, 2, and 3, respectively. The MEP search region was set to be the same as the QM region.

US-MD simulations were performed along the MEP. The main interatomic distances that contribute to the reaction are considered to be the CVs. The CVs used for each step are illustrated in Fig. 2 by colored lines. 32, 16, and 32 windows were set in steps 1, 2, and 3, respectively. The midpoint values of each window (ξ₀) are listed in Tables S4–S6. The force constant of the umbrella potential was set to 200.0 kcal mol⁻¹Å⁻² in steps 1 and 3, and 100.0 kcal mol⁻¹Å⁻² in step 2. QM/MM-MD simulations were performed for 6 ps per window at the B3LYP-D3/aug-cc-pVDZ level of theory in the NVT ensemble. The proton that undergoes proton transfer reaction was unconstrained, and thus the system was propagated with a time step of 0.5 fs using the VVER integrator. The temperature was maintained at 300 K by using a Bussi thermostat. The structure was sampled every ten steps, and the trajectory data of the last 5 ps were analyzed to obtain the PMF. The GENESIS analysis tools were used to perform path CV, MBAR and PMF analyses.

The energy profiles of each step in the β-reaction stage I of the Open/Open state obtained by the MEP search are shown in Fig. 3(a) and Table II, and the conformations of the reactant, transition state (TS), and product are shown in Fig. 4. The energy profile obtained by the MEP search is the average of three calculations using different initial coordinates obtained by equilibration.

In the first step, proton H1 bound to the donor of the proton transfer, carbon C1 (d_H1−C1 = 1.10 Å), of E(Aex₁) in the reactant (I in Fig. 4) is transferred to the acceptor, nitrogen NZ, of βK87 (d_H1−NZ = 1.03 Å), yielding intermediate E(Q₁*) (III in Fig. 4). In this process, βK87 is protonated, and hydrogen HZ1 of βK87 forms a hydrogen bond with the oxygen O7 of the phosphate group. In this process, the carbon, C1, hybridization state of the E(Aex1) changes sp² from sp³. In the TS (II in Fig. 4), the proton, H1 is located in the middle of the donor, C1, and the acceptor, NZ. The barrier of the reaction is calculated to be 26.7 kcal/mol. This step is the rate determining step of the β-reaction stage I. This is mirrored in the variation of the atomic distances r₁, r₂, and r₄ along the MEP, as shown in Fig. 2(c) (left). Note that the hydrogen bond between oxygen O1 of E(Aex₁) and oxygen O of βG303 remained intact.

In the second step, we expected that the dihedral angle defined by atoms N1-C1-C2-O1 would rotate 180°, which places the hydroxy group close to the acid-base catalytic group, βK87, and facilitates subsequent proton transfer from βK87 to E(Q₁) in the third step [see Fig. 1(b)]. However, the result of the MEP search shows that the dihedral angle rotates only 150° and that the hydrogen bond donor, oxygen O1, of the hydroxy group, is located far from the acceptor, nitrogen NZ, and proton H1 of βK87 (d_H1−O1 = 4.05 Å) (V in Fig. 4). In Fig. 2(c) (middle), the distance between O1 and OD1 of βD305, represented by r₆, is kept around 2.6 Å and the distance between O1 and NZ, represented by r₂, is larger than 4.5 Å. This result suggests that a stable hydrogen bond is formed between βD305 and the hydroxy group of the E(Aex₁), and it stabilize E(Aex₁) for an inactive state. The hydrogen bond inhibits the rotation of the hydroxy group. Consequently, the intermediate E(Q₁) has an unfavorable conformation for the proton transfer in the next step.

The third step consists of two sub-steps. In the first sub-step, the hydroxy group of the intermediate E(Q₁) is shifted close to βK87 by breaking the hydrogen bond between βD305 and the hydroxy group, forming a new hydrogen bond between the hydroxy group and βK87 (d_H1−O1 = 1.24 Å) (V* in Fig. 4). In the second sub-step, the proton (H1) of βK87 is transferred to the hydroxy group, thereby forming a water molecule consisting of H4-O1-H1 (VII in Fig. 4) with the bond cleavage of the substrate C–O bond. The released water forms a hydrogen bond with the nitrogen NZ of βK87 (d_H1−NZ = 1.81 Å). r₁, r₂, and r₃ in Fig. 2 (right) show that the hydroxy group initially approached to βK87, followed by a proton transfer from nitrogen NZ to oxygen O1. Note that the hydrogen bond between βK87 and the phosphate group represented by r₅ is kept in step 3.

The hydrogen bond between the oxygen OD1 of βD305 and hydrogen H4 of the hydroxy group is crucial because it inhibits the smooth dihedral angle rotation of the hydroxy group and the second proton transfer from βK87 to the hydroxy group. This hydrogen bond was found in the x-ray crystal structures of wild-type (WT) TRPS not only in the Apo|PLS (PDBID: 1KFJ¹³) but also in F6F|PLS (PDBID: 2CLM²⁴), F9F|PLS (PDBID: 2CLL²⁴), and F19F|PLS (PDBID: 2CLO²⁴), as shown in Fig. S4. Furthermore, the same hydrogen bond was present in the TRPS mutant (βK87T) of TRPS in the Apo|PLS (PDBID: 1UBS³⁵) and mutant (αT183V) in the Apo|PLS (PDBID: 1KFE¹³). The hydrogen bond did not break during step 2, as shown in Fig. 4. These results indicate that the hydrogen bond is stable in E(Aex₁), E(Q₁*), and E(Q₁).

The US-MD simulations at the B3LYP-D3/aug-cc-pVDZ level yield the PMF along the MEP in the β-reaction stage I, as shown in Fig. 3(b) and Table III. The PMF is 1/2–2/3 smaller than the potential energy (PE) changes. The dynamics of TRPS accounted for by the MD simulations, as well as the difference in the basis set (cc-pVDZ and aug-cc-pVDZ), cause the difference between PMF and PE. The relative energies of the β-reaction stage I between reactant E(Aex₁) and the product E(A-A) is endothermic by 12.0 kcal/mol in the PE, whereas it is exothermic by 3.5 kcal/mol in the PMF. The first step of the β-reaction stage I, which is a rate-determining step, is an endothermic reaction, and the reaction energy by the US-MD simulations is 14.0 kcal/mol and the activation barrier is 21.0 kcal/mol. The second step is an endothermic reaction, for which the reaction energy by the US-MD simulations is 4.3 kcal/mol and the activation barrier is 2.6 kcal/mol. Teixeira et al. reported that the second step of the β-reaction stage I in the IGP| E(Q₁*) is an exothermic reaction.⁷¹ The hydrogen bond βD305-E(Q₁*) in the second step may account for this discrepancy. The third step is an exothermic reaction, and the reaction energy is −21.8 kcal/mol. The activation energy from the intermediate state after the first sub-step to the reactant [E(A–A, VII)] is −11.5 kcal/mol.

In the x-ray crystal structure (PDBID: 2J9X²⁴) of G3P|P1T, where P1T has a similar conformation to E(A-A), βD305 forms a hydrogen bond with βR141 instead of P1T. Therefore, we carried out an MEP search and a US-MD simulation with an initial structure of the Open/Closed state, in which βD305 forms a salt bridge with βR141, to clarify the role of the hydrogen bond between βD305 and E(Q₁*) in the second step of the β-reaction stage I. The MEP search and the US-MD simulation procedures used were the same as previously described. Figures 5(a) and 5(b) show the energy profile of the three steps in the β-reaction stage I in the presence of the salt bridge, βD305-βR141 as shown in Table II and Table III. Figures 6(a) and 6(b) show the activation and relative energies of steps 1-3. The results show that the rate-determining step is step 1, which is the same as the calculation with βD305-E(Aex₁) and the reaction at the second step is exothermic with the reaction energy of −6.2 kcal/mol when βD305-βR141 is formed, whereas the reaction is endothermic with the reaction energy of 4.4 kcal/mol when βD305-E(Q₁*) is formed.

In the second step, the dihedral angle defined by atoms N1-C1-C2-O1 was smoothly rotated from E(Q₁*) [III in Fig. 5(c) to E(Q₁] [V in Fig. 5(c)] via the TS [IV in Fig. 5(c)]. E(Q₁) is stabilized by the hydrogen bond between hydrogen H1 of βK87 and oxygen O1 of E(Q1) (d_H1−O1 = 1.16 Å). Therefore, the reaction in the second step is exothermic.

In the third step, the conformation of the intermediate E(Q₁) of step 2 is favorable for the proton transfer reaction from the proton donor, nitrogen NZ of βK87 to the acceptor, oxygen O1 of E(Q₁). The hydroxy group was already close to the acid-base catalytic group, βK87 after the dihedral angle rotation in the second step, as shown in Fig. 7. This conformation smoothly connects to the third step without switching from the open to the closed state by moving the hydroxy group toward βK87. Therefore, in the third step, there were no sub-steps as compared with the simulations with the hydrogen bond between βD305 and E(Q₁). These results show that the conformational change from the open to the closed state and hydrogen bond switching from βD305-E(Aex₁) to βD305-βR141 is an important factor in regulating the β-reaction stage I.

In this study, QM/MM calculations of the β-reaction stage I in TRPS were performed in the Open/Open and the Open/Closed state. MEP searches and US-MD simulations at the B3LYP-D3/cc-pVDZ and aug-cc-pVDZ levels, respectively, in the Apo|E(Aex₁) show the reaction pathways, energy profiles, and detailed conformations of the reactant/TS/product in the β-reaction stage I. Teixeira et al. reported the energy profiles and conformations of the reactant/TS/product in the whole reaction scheme of TRPS using the IRC calculations and the vibrational frequency calculations at the B3LYP/6-31++g(2df, pd) level.⁷¹ Our PMF from the US-MD simulations is in good agreement with their results, except for the second step of the β-reaction stage I. Our result shows an endothermic reaction, whereas their result is exothermic. This difference appears to result from the initial coordinates of the TRPS in the two calculations. Teixeira et al. used the x-ray crystal structure of the Closed/Closed state of the F9F|7MN (PDBID: 3PR2⁹⁹), which is similar to that of IGP|E(Q₂). It contains a salt bridge, βD305-βR141. Meanwhile, the initial coordinates of our calculations were based on the x-ray crystal structure in the Apo|PLS states (PDBID:1KFJ¹³), which forms a different state, Open/Open, with hydrogen bond, βD305-PLS.

Based on the small local conformational changes and large differences in the energy profile, we consider that this might be important to answer the question in the allosteric regulation of TRPS: how IGP binding at the α-active site can enhance the reaction activity of the β-reaction stage I at the distal β-active site. There are many x-ray crystal structures of wild-type (WT) and mutants of TRPS, whose α-active site is either apo or in one of the IGP-analog bound states. Figure S4 shows the WT for Apo|PLS (PDBID: 1KFJ¹³), F6F|PLS (PDBID: 2CLM²⁴), F9F|PLS (PDBID: 2CLL²⁴), F19F|PLS (PDBID: 2CLO²⁴), αD60N for IPL|PLS (PDBID: 1BEU¹⁰⁰), αT183V for Apo|PLS (PDBID: IKFE¹³), βK87T for Apo|PLS (PDBID: 1UBS³⁵), G3P|PLS (PDBID: 2TSY³⁵), and IPL|PLS (PDBID: 2TRS³⁵). Interestingly, the hydrogen bond between βD305 and E(Aex₁) is present in six out of the nine x-ray crystal structures listed above, in both the apo (1KFJ,¹³ 1KFE,¹³ and 1UBS³⁵) and in the IGP-analog bound states (2CLM,²⁴ 2CLL,²⁴ and 2CLO²⁴). By contrast, the salt bridge between βD305 and βR141 was observed only in the IGP-analog bound states (2TSY³⁵ and 2TRS³⁵). This structural evidence suggests that the free-energy barrier for the allosteric switching from the open to the closed state in the β-subunit involving the hydrogen bond switching is marginal.

To understand the allosteric regulation mechanisms for the β-reaction stage I in TRPS in more detail, we carried out a number of classical MD simulations for 50 ns based on the x-ray crystal structure in the Apo|E(Aex₁) state based on PDBID: 1KFJ,¹³ which forms a hydrogen bond between βD305 and the ligand. In all simulations, the hydrogen bond was not stable and was switched to that between βD305 and βR141 within tens of nanoseconds, as shown in Fig. S5. The salt bridge, βD305–βR141 has been considered as a measure to examine whether the COMM domain in the β-subunit is open or closed.^12,101–103 However, in our recent MD/US study,⁷⁰ the salt bridge did not show a good correlation with the conformational changes of the COMM domain for opening/closing. Note that in our simulations, the open and close conformations in the β-subunit are defined based on the root mean square deviation (RMSD) of the COMM domain, which was measured by aligning the β-subunit except for the COMM domain in reference to an x-ray crystal structure 2J9X. Then, the COMM domain is referred to as open conformation when the RMSD is larger than 2.0 Å. Instead, the hydrogen bonding networks at the interface of the two subunits and between the β-subunit and β-ligand were shown to be more essential determinants of the TRPS structure. Based on these results, we can consider two hypotheses on the allosteric regulation of the β-reaction stage I: (i) although the crystal structure used in the study is categorized as an Open|Open state, the β-subunit structure (or the orientation of the COMM domain) can be more wide open in solution. Note that our previous MD simulations, in which the system was modeled based on the x-ray crystal structures 1K8X (Apo|PLP), 1WBJ (G3P|PLP), and 2J9X (G3P|P1T) found that the COMM domain of Apo|E(Ain), IGP|E(Ain), and IGP|E(A–A) is in the open, widely open, and closed forms, respectively, in solution.⁷⁰ Therefore, the solution structure of Apo|E(Aex₁) requires a careful examination using a similar strategy. If Apo|E(Aex₁) showed a more opened structure, the distance between βD305 and βR141 would be longer, suggesting an increased possibility of forming a hydrogen bond between βD305 and the β-ligand. (ii) Another possibility is that the solution structure in the Apo|E(Aex₁) state shows a closed conformation due to L-Ser at the β-active site. Even in this case, the salt bridge interaction between βD305 and βR141 can fluctuate, as shown in the Apo|E(Ain) and IGP|E(A-A) states in our previous study. In this case, longer MD simulations or those with enhanced conformational sampling are necessary to examine the marginal stability between the two hydrogen bonds in the presence and absence of the α-ligand at the binding site. In either case, the possibility of conformational dynamics in the solution structures of TRPS with different ligand-bound states should not be neglected. X-ray crystal structures give an atomistic description of the active sites in proteins at a moment in time. Meanwhile, it is important to consider dynamic aspects using MD and QM/MM simulations as we did in the previous and current studies.

In summary, we learned from the QM/MM calculations that the allosteric regulation in TRPS in which reactivity at the β-active sites upon the binding of a ligand at the distal α-active site is increased can be explained via a local interaction change, such as the allosteric switching from the open to closed state involving the hydrogen bond switching from βD305-ligand to βD305-βR141. Existing crystal structures of TRPS suggest that the former hydrogen bond is marginally stable in Apo|PLS as an inactive state, whereas the latter is found only in the presence of a ligand at the α-active site. Computational and structural information suggest that inter-subunit communication is regulated dynamically and can change a few local but essential interactions in the vicinity of the distal reaction active site. Such dynamic regulations of local and global structures can be applicable to many other allosteric enzymes, which should be examined by structural and computational studies.

See the supplementary material for the scheme of the Trp synthesis of TRPS and the computational details of the equilibration and US-MD. Movies of the β-reaction stage I in the Apo|E(Aex₁) of TRPS with the hydrogen bond in the βD305-ligand or βD305-βR141 have been provided.

This research was partially supported by RIKEN Pioneering Research Projects (Dynamic Structural Biology/Glycolipidologue Initiative) (to Y.S.), RIKEN Incentive Research Projects (to S.I.), Program for Promoting Research on the Supercomputer Fugaku [Biomolecular dynamics in a living cell (Grant No. JPMXP1020200101)/MD-driven Precision Medicine (Grant No. JPMXP1020200201)], MEXT/KAKENHI Grant Nos. JP19H05645, JP21H05249 (to Y.S.), JP20H02701, and JP22H04761 (to K.Y.). We used the computer system HOKUSAI (Project ID: Q22535), provided by the RIKEN Information System Division, and Oakbridge-CX, provided by the University of Tokyo (hp220110) through the HPCI System Research Project.

The authors have no conflicts to disclose.

Shingo Ito: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (equal); Visualization (lead); Writing – original draft (lead); Writing – review & editing (supporting). Kiyoshi Yagi: Conceptualization (equal); Investigation (supporting); Methodology (lead); Project administration (equal); Software (lead); Supervision (equal); Validation (lead); Writing – review & editing (equal). Yuji Sugita: Conceptualization (equal); Funding acquisition (lead); Project administration (lead); Resources (lead); Supervision (lead); Writing – review & editing (lead).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

A|B

A and B indicate the ligands/states of the α- and β-active sites

COMM

Communication

E(A–A)

α-aminoacrylate intermediate

E(Aex₁)

External aldimine

E(Ain)

Internal aldimine

E(Q₁)

Quinonoid intermediate after the β-reaction step 2

E(Q₁*)

Quinonoid intermediate before the β-reaction step 2

F9F

N-(4′-trifluoromethoxybenzenesulfonyl)-2-aminoethyl phosphate

G3P

Glyceraldehyde-3-phosphate

IGP

3-indole-D-glycerol-3′-phosphate

L-Ser

L-Serine

M.tb

Mycobacterium tuberculosis

MD

Molecular dynamics

MEP

Minimum-energy path

PLP

Pyridoxal phosphate

PLS

Serine pyridoxal phosphate Schiff base

TB

Tuberculosis

Trp

Tryptophane

TRPS

Tryptophan synthase

US

Umbrella sampling

αL6

Flexible loop 6 in the α-active site