Insight into Structural Characteristics of Protein-Substrate Interaction in Pimaricin Thioesterase

As a polyene antibiotic of great pharmaceutical significance, pimaricin has been extensively studied to enhance its productivity and effectiveness. In our previous studies, pre-reaction state (PRS) has been validated as one of the significant conformational categories before macrocyclization, and is critical to mutual recognition and catalytic preparation in thioesterase (TE)-catalyzed systems. In our study, molecular dynamics (MD) simulations were conducted on pimaricin TE-polyketide complex and PRS, as well as pre-organization state (POS), a molecular conformation possessing a pivotal intra-molecular hydrogen bond, were detected. Conformational transition between POS and PRS was observed in one of the simulations, and POS was calculated to be energetically more stable than PRS by 4.58 kcal/mol. The structural characteristics of PRS and POS-based hydrogen-bonding, and hydrophobic interactions were uncovered, and additional simulations were carried out to rationalize the functions of several key residues (Q29, M210, and R186). Binding energies, obtained from MM/PBSA calculations, were further decomposed to residues, in order to reveal their roles in product release. Our study advanced a comprehensive understanding of pimaricin TE-catalyzed macrocyclization from the perspectives of conformational change, protein-polyketide recognition, and product release, and provided potential residues for rational modification of pimaricin TE.

As a potent polyene antibiotic produced by Streptomyces natalensis, pimaricin (i.e. natamycin) primarily functions in the treatment of fungal infections caused by Candida, Fusarium, Penicillium, and Aspergillus organisms [10]. It is also known as an additive in food industry [11]. Pimaricin was approved by the Food and Drug Administration (FDA) as a drug for fungal keratitis in 1978 [12]. Ergosterol constitutes a major component in fungal and trypanosomatids plasma membrane, while absent in animal cells [13]. Pimaricin serves to bind specifically to ergosterol, downregulate vesicle trafficking, suppress membrane protein transport, and interfere with endocytosis, as well as exocytosis without permeabilizing the membrane [14][15][16]. Its strong performance in clinical trials makes pimaricin appealing to researchers, and its biosynthetic pathway modification and drug design have become new science hotspots [17].
Pimaricin is synthesized by type I polyketide synthases (PKSs), which consists of several covalently-connected multi-domain "modules." Each module contains a set pattern of domains, with acyltransferase (AT) adding acyl-CoA building blocks, acyl carrier protein (ACP) carrying the polyketide between modules, and ketosynthase (KS) accepting the growing chain from ACP [18]. An extra combination of domains, such as ketoreductase (KR), dehydratase (DH), methyltransferase (MT) were responsible for the production of distinctive macrolactones [19][20][21]. Situated in the last module, the thioesterase (TE) domain off-loads the ACP-tethered polyketide from PKS via macrocyclization or hydrolysis.
Consistent with serine hydrolase, a two-step mechanism has been proposed for TE-mediated catalysis of macrocyclic polyketides [22]. The first step is acylation of a universally conserved serine residue in the catalytic triad, generating an acyl-enzyme intermediate and stabilized for a considerable time [23]. The second step takes place with an intra-molecular hydroxyl group on polyketide which initiates a nucleophilic attack and leads to cyclization, or hydrolysis of the acylenzyme intermediate with no efficient intra-molecular nucleophile present.
In our previous work concerning 6-deoxyerythronolide B synthase (DEBS)-TE [24], a hydrogen bond emerged between the polyketide chain terminal hydroxyl group O11 and carbonyl oxygen O (Figure 1), as accompanied by the swing of C11 ethyl of substrate. This structure has been reported in Trauger's work in 2001, where it was referred to as the "pre-organization state" (POS). According to Trauger [25], the "product-like" conformation might contribute to pre-organization of the substrate for cyclization. The conformation with a hydrogen bond, forming between the O11 and Nε of His259 in the catalytic triad, was defined as an active state in our study. This conformation maintained for ~100 ns in our simulations with considerable steadiness. However, the distance of O11-C1 for nucleophilic attack was larger than 6 Å in active state. Finally, an advantageous conformer (the prereaction state, PRS) was found [24] after ~270 ns MD simulation, which possessed both hydrogen bond O11-NεH259 and an appropriate distance between O11 and Cl to facilitate nucleophilic attack. Critical to mutual recognition and catalytic preparation between TE and substrate, the PRS seemed decisive in the occurrence of macrocyclization. To understand the molecular basis of pimaricin-TE (pima-TE) catalyzed macrocyclization, MD simulations were employed on enzyme-substrate complex. POS, active state, and PRS were discovered during 5 × 50 ns molecular dynamics (MD) simulations, and the conformational transition between POS and PRS was explored. The structural characteristics of POS and PRS were uncovered by conducting analyses of hydrogen-bonding and hydrophobic interactions. Additional simulations on several mutants (including Q29A, M210G, R186F, R186Y and S138C) were carried out to validate the functions of several key residues in substrate recognition and product release. At last, the binding energies of enzyme-product complex were obtained through MM/PBSA calculations, and with critical residues highlighted. We also provided an explanation on the departure of product from the active site. To understand the molecular basis of pimaricin-TE (pima-TE) catalyzed macrocyclization, MD simulations were employed on enzyme-substrate complex. POS, active state, and PRS were discovered during 5 × 50 ns molecular dynamics (MD) simulations, and the conformational transition between POS and PRS was explored. The structural characteristics of POS and PRS were uncovered by conducting analyses of hydrogen-bonding and hydrophobic interactions. Additional simulations on several mutants (including Q29A, M210G, R186F, R186Y and S138C) were carried out to validate the functions of several key residues in substrate recognition and product release. At last, the binding energies of enzyme-product complex were obtained through MM/PBSA calculations, and with critical residues highlighted. We also provided an explanation on the departure of product from the active site.

Key Structural Conformations in MD Simulations
Intrigued by the recognition mechanism of pima-TE, 5 × 50 ns MD simulations were performed on constructed complex. Root-mean-square deviation (RMSD) analysis revealed that all five trajectories attained equilibrium ( Figure S1). The root-mean-squared fluctuation (RMSF) values highlighted six parts on pima-TE. Firstly, the lid region was violently jacked up by the erected polyketide ( Figure 2). As a polyene molecule with 26-atom skeleton, pimaricin's accommodation would require a larger space, compared with pikromycin, a 14-membered ring. It was conceivable that the relaxation of the substrate would incur close contact with the lid. Next, as components of the channel, αL2, as well as loop l1 & l2, presented evident structural dynamism and various coiling states, while αL3 exhibited negligible fluctuation. Helix αL2 was proposed to wield a larger influence on protein-intermediate recognition than αL3, and l1 and l2 were responsible for the exit and entrance size. At last, RMSF indicated that loop l3 adopted larger fluctuations than loop l1 and l2, and the b-factor calculation [26,27] disclosed an inherent mobility of loop l3.

Key Structural Conformations in MD Simulations
Intrigued by the recognition mechanism of pima-TE, 5×50 ns MD simulations were performed on constructed complex. Root-mean-square deviation (RMSD) analysis revealed that all five trajectories attained equilibrium ( Figure S1). The root-mean-squared fluctuation (RMSF) values highlighted six parts on pima-TE. Firstly, the lid region was violently jacked up by the erected polyketide ( Figure 2). As a polyene molecule with 26-atom skeleton, pimaricin's accommodation would require a larger space, compared with pikromycin, a 14-membered ring. It was conceivable that the relaxation of the substrate would incur close contact with the lid. Next, as components of the channel, αL2, as well as loop l1 & l2, presented evident structural dynamism and various coiling states, while αL3 exhibited negligible fluctuation. Helix αL2 was proposed to wield a larger influence on protein-intermediate recognition than αL3, and l1 and l2 were responsible for the exit and entrance size. At last, RMSF indicated that loop l3 adopted larger fluctuations than loop l1 and l2, and the bfactor calculation [26,27] disclosed an inherent mobility of loop l3. Next, conformational variations at active site in each trajectory were carefully studied. The entire 250 ns trajectory was divided into three categories, based on distance measurement. With a hydrogen bond coming into being between terminal hydroxyl O7 and NεH261 (distance (O7-NεH261) ≤ 3.0Å), the intermediate was regarded ready to be de-protonated by H261, namely, an active state. The active state Next, conformational variations at active site in each trajectory were carefully studied. The entire 250 ns trajectory was divided into three categories, based on distance measurement. With a hydrogen bond coming into being between terminal hydroxyl O 7 and Nε H261 (distance (O 7 -Nε H261 ) ≤ 3.0 Å), the intermediate was regarded ready to be de-protonated by H261, namely, an active state. The active state was observed in all five trajectories (8.7, 3.1, 17.1, 79.5, and 23.4%, respectively), with the highest proportion in md4 ( Figure 3). Moreover, the terminal O 7 was proposed to be conducive for nucleophilic attack onto carboxyl C 1 with distance (O 7 -C 1 ) ≤ 4.5 Å. The PRS was defined as both criteria were met, and was present in md4 for 4700 frames (18.8%, Figure S2); in other trajectories, PRS appeared with a significantly lower frequency, testifying to its unsteadiness as a transient state.  Figure 3). Moreover, the terminal O7 was proposed to be conducive for nucleophilic attack onto carboxyl C1 with distance (O7-C1) ≤ 4.5Å. The PRS was defined as both criteria were met, and was present in md4 for 4700 frames (18.8%, Figure S2); in other trajectories, PRS appeared with a significantly lower frequency, testifying to its unsteadiness as a transient state. Distinguished from PRS and active state that ultimately lead to macrocyclization, inactive conformations are susceptible to hydrolysis. Notably, among inactive conformations, the POS, which is characterized by a hydrogen bond between O7 and carbonyl oxygen O1 of substrate, was also observed within md1 for 11896 frames (47.6%, Figure S2), whereas it was nearly absent in others (3, 2, 20 and 0 frames in md2-5).
The conformational flip took place in 18-22 ns of md1 trajectory, with conformation altering progressively from PRS (−180°) to POS (−60°). As presented in Figure 4, terminal hydroxyl O7 firstly swung up and disassociated from H261, followed by Cβ-Cγ twisting clockwise and terminal methyl oriented towards the entrance (I→II). Further, the intermediate swelled to diminish distance O7-O1. After quick adjustment, POS came into being and maintained for rest of the trajectory (II→III). Distinguished from PRS and active state that ultimately lead to macrocyclization, inactive conformations are susceptible to hydrolysis. Notably, among inactive conformations, the POS, which is characterized by a hydrogen bond between O 7 and carbonyl oxygen O 1 of substrate, was also observed within md1 for 11896 frames (47.6%, Figure S2), whereas it was nearly absent in others (3, 2, 20 and 0 frames in md2-5).

Conformational Transition between POS and PRS
Next, the transformation between POS and PRS was studied using dihedral angle C α -C β -C γ -O 7 as an indicator of polyketide terminal rotation. In PRS, bond C α -C β ran anti-parallel against C γ -O 7 (−180 • ), but in POS, the dihedral angle was altered to an acute angle fluctuating between (−30 • , −70 • ).
The conformational flip took place in 18-22 ns of md1 trajectory, with conformation altering progressively from PRS (−180 • ) to POS (−60 • ). As presented in Figure 4, terminal hydroxyl O 7 firstly swung up and disassociated from H261, followed by C β -C γ twisting clockwise and terminal methyl An energy calculation was also conducted to investigate the structural stability of aforementioned conformations. As expected, POS harbored a lower energy than PRS by 4.58 kcal/mol, indicating the steadiness of the O7-O1 intra-molecular hydrogen bond. On the other hand, the active state was calculated to be 0.18 kcal/mol less stable than PRS. The slight difference prompted that conformational transition between active state and PRS would easily achieve through Cβ-Cγ bond rotation.
In conclusion, a conformational transformation between POS and PRS was accomplished through dihedral flip and conformation adjustment, and the energies on POS, active state, and PRS were computed to understand the reaction process.

Hydrophilic and Hydrophobic Interactions in Pima-TE System
Based on MD simulations, hydrogen bonding and hydrophobic interactions between pima-TE and substrate were studied. As exhibited in Figure 5, in PRS, loop l1 (residue 170-177) played a crucial part in fastening the substrate. The atoms O2, O3, O4 or O5 were anchored by H172 (13.35%), T177 (15.20%) and Q174 (4.55%) without fixed pattern. Residue Q29, stretching downward from the lid region, served as a crane to lift up O6 and gave rise to an erected molecule (39.26%). The main chain of M210 fixed O4 at the center of the molecule (28.56%), while its side chain laid parallel to the hydrophobic area of conjugated polyene (99.92%). The hydrophobic segments of T73, L183 and Y180 were also oriented towards the polyketide chain and worked jointly to conserve a water-free subenvironment with frequencies of 94.88, 89.27 and 78.50%, respectively.
An energy calculation was also conducted to investigate the structural stability of aforementioned conformations. As expected, POS harbored a lower energy than PRS by 4.58 kcal/mol, indicating the steadiness of the O 7 -O 1 intra-molecular hydrogen bond. On the other hand, the active state was calculated to be 0.18 kcal/mol less stable than PRS. The slight difference prompted that conformational transition between active state and PRS would easily achieve through C β -C γ bond rotation.
In conclusion, a conformational transformation between POS and PRS was accomplished through dihedral flip and conformation adjustment, and the energies on POS, active state, and PRS were computed to understand the reaction process.

Hydrophilic and Hydrophobic Interactions in Pima-TE System
Based on MD simulations, hydrogen bonding and hydrophobic interactions between pima-TE and substrate were studied. As exhibited in Figure 5, in PRS, loop l1 (residue 170-177) played a crucial part in fastening the substrate. The atoms O 2 , O 3 , O 4 or O 5 were anchored by H172 (13.35%), T177 (15.20%) and Q174 (4.55%) without fixed pattern. Residue Q29, stretching downward from the lid region, served as a crane to lift up O 6 and gave rise to an erected molecule (39.26%). The main chain of M210 fixed O 4 at the center of the molecule (28.56%), while its side chain laid parallel to the hydrophobic area of conjugated polyene (99.92%). The hydrophobic segments of T73, L183 and Y180 were also oriented towards the polyketide chain and worked jointly to conserve a water-free sub-environment with frequencies of 94.88, 89.27 and 78.50%, respectively.
Under the circumstance of POS, Q29 (19.38%) and H172 (10.71%) assisted the intermediate erection, except for this time the molecule slighted twisted to stabilize the intra-molecular hydrogen bond, which enabled S33 from the lid region to contribute in the hydrogen bonding network (8.90%). Identical to PRS, M210 again assumed the role of both a hydrophilic stake (56.98%) and a hydrophobic driving force (99.28%). Besides Y180 (36.17%) and L183 (57.06%) as in PRS, L213 also participated in hydrophobicity maintenance (25.43%). Under the circumstance of POS, Q29 (19.38%) and H172 (10.71%) assisted the intermediate erection, except for this time the molecule slighted twisted to stabilize the intra-molecular hydrogen bond, which enabled S33 from the lid region to contribute in the hydrogen bonding network (8.90%). Identical to PRS, M210 again assumed the role of both a hydrophilic stake (56.98%) and a hydrophobic driving force (99.28%). Besides Y180 (36.17%) and L183 (57.06%) as in PRS, L213 also participated in hydrophobicity maintenance (25.43%).
In a word, binding modes of pimaricin polyketide with TE shared considerable similarity between PRS and POS, with Q29, M210 and residues on loop l1 interacting with the chain via hydrogen bonding, and M210, Y180, and L183 contributing to hydrophobic network.

Mutation 1-Q29A
According to the analyses of wild type simulations, Q29, located at lid region of pima-TE, could mediate the distance O7-NεH261 within a favorable range through bonding with O6 of substrate. When mutated to Ala, with Q29's side chain shortened and hydrogen bond abolished, it was speculated that the substrate would fall off from its original position. Here, the distance between O6 and Cα of Q/A29 (designated as O6-CQ/A29) was utilized to depict the substrate's spatial displacement. As shown in Figure 6, the distance O6-CQ/A29 fluctuated acutely in Q29A trajectory md2 and md3, with the substrate either overlength (2I, 3I), hydrogen bonding to other residues (3II), or drifting aimlessly (2II, 3III). The erratic change also decreased PRS formation by a large margin (4.37% vs. 0.54%, Table S1). On the other hand, POS was observed in Q29A md1 with a frequency of 80.84% (1I, 1II).  [28]. Backbone of the polyketide substrate was colored in yellow, and residues providing hydrogen bonding and hydrophobic interactions in slate_blue and brown.
In a word, binding modes of pimaricin polyketide with TE shared considerable similarity between PRS and POS, with Q29, M210 and residues on loop l1 interacting with the chain via hydrogen bonding, and M210, Y180, and L183 contributing to hydrophobic network.

Mutation 1-Q29A
According to the analyses of wild type simulations, Q29, located at lid region of pima-TE, could mediate the distance O 7 -Nε H261 within a favorable range through bonding with O 6 of substrate. When mutated to Ala, with Q29's side chain shortened and hydrogen bond abolished, it was speculated that the substrate would fall off from its original position. Here, the distance between O 6 and Cα of Q/A29 (designated as O 6 -C Q/A29 ) was utilized to depict the substrate's spatial displacement. As shown in Figure 6, the distance O 6 -C Q/A29 fluctuated acutely in Q29A trajectory md2 and md3, with the substrate either overlength (2 I , 3 I ), hydrogen bonding to other residues (3 II ), or drifting aimlessly (2 II , 3 III ). The erratic change also decreased PRS formation by a large margin (4.37% vs. 0.54%, Table S1). On the other hand, POS was observed in Q29A md1 with a frequency of 80.84% (1 I , 1 II ).

Mutation 2-M210G
To validate M210's function in hydrophobic interaction network, M210 was mutated into Gly. The substrate backbone's distance RMSD (dRMSD) was calculated with the first frame as a reference. As seen in Figure 7, the larger dRMSD signified variation in substrate conformation, and its irregularity suggested volatility. Furthermore, new patterns of hydrogen bonding were observed in mutant (Figure 7): In md1, the polyketide chain leaned towards αL3 and interacted with N214 (23.08%); in md2 and md3, the substrate slightly rotated and bonded with D179 on αL2 (76.35% and 87.10%). Having lost M210 as a hydrophobic barrier, the polyketide chain would adjust its position, and M206 from the neighboring cycle of αL3 exhibited hydrophobicity. Owing to the altered interaction network, it was hard for the substrate to attain PRS as in wild type complex (4.37% vs. 0.72%).
Specifically, the aforementioned conformational change of substrate also produced a similar effect on loop l1 seeking hydrogen bonding, and gave rise to a shrinking channel exit, while a bigger exit would be favored in product release. Taken together, M210 was crucial in maintaining the polyketide chain in between αL2 and αL3, which was conduced to protein-substrate interaction and an advantageous channel exit shape.

Mutation 2-M210G
To validate M210's function in hydrophobic interaction network, M210 was mutated into Gly. The substrate backbone's distance RMSD (dRMSD) was calculated with the first frame as a reference. As seen in Figure 7, the larger dRMSD signified variation in substrate conformation, and its irregularity suggested volatility. Furthermore, new patterns of hydrogen bonding were observed in mutant ( Figure 7): In md1, the polyketide chain leaned towards αL3 and interacted with N214 (23.08%); in md2 and md3, the substrate slightly rotated and bonded with D179 on αL2 (76.35% and 87.10%). Having lost M210 as a hydrophobic barrier, the polyketide chain would adjust its position, and M206 from the neighboring cycle of αL3 exhibited hydrophobicity. Owing to the altered interaction network, it was hard for the substrate to attain PRS as in wild type complex (4.37% vs. 0.72%).
Specifically, the aforementioned conformational change of substrate also produced a similar effect on loop l1 seeking hydrogen bonding, and gave rise to a shrinking channel exit, while a bigger exit would be favored in product release. Taken together, M210 was crucial in maintaining the polyketide chain in between αL2 and αL3, which was conduced to protein-substrate interaction and an advantageous channel exit shape.

Mutation 3-R186F & R186Y
As a residue containing multiple hydrophilic groups, R186 bonded with O7 for a rather high probability in wild type complex simulations. As seen in Figure 8a, in md1-3, high frequency of O7-NR186 bonding could account for the scarce existence of O7-NεH261 interaction. To promote PRS formation, R186 was firstly mutated to Phe.

Mutation 3-R186F & R186Y
As a residue containing multiple hydrophilic groups, R186 bonded with O 7 for a rather high probability in wild type complex simulations. As seen in Figure 8a, in md1-3, high frequency of O 7 -N R186 bonding could account for the scarce existence of O 7 -Nε H261 interaction. To promote PRS formation, R186 was firstly mutated to Phe.
To our disappointment, the frequency of PRS formation did not improve (4.37% vs. 1.64%). It was determined that E80 was coupled with R186 and R266, to pose a spatial barrier at the entrance and prevent the admission of other substrate, while functionally maintaining the closure and hydrophobicity of substrate pocket. Nevertheless, when R186 was mutated to Phe, a crack appeared (Figure 8b), and frequency of E80-R266 interaction was lowered as well (Table S2). Worse still, lacking the tying force, the distance Cα F186 -Cα E80 also increased, implying a larger entrance (Figure 8c).
Based on our findings, pima-TE was re-modified into R186Y mutant. This time, we endowed a hydroxyl to the side chain of mutated residue to bond with E80, while the remainder stayed hydrophobic. We were more than pleased to find a significant rise in PRS ratio (4.37% vs. 18.14%), with Y186-E80 bonding partly restored (Table S3). Of particular note, a close-to-reaction PRS conformation appeared in md3 and maintained for over 10 ns, with the terminal methyl oriented towards the entrance and forcing O 7 closer to C 1 (Figure 8d). We thus regard R186Y as a promising modification towards pimaricin productivity advancement. To our disappointment, the frequency of PRS formation did not improve (4.37% vs. 1.64%). It was determined that E80 was coupled with R186 and R266, to pose a spatial barrier at the entrance and prevent the admission of other substrate, while functionally maintaining the closure and hydrophobicity of substrate pocket. Nevertheless, when R186 was mutated to Phe, a crack appeared (Figure 8b), and frequency of E80-R266 interaction was lowered as well (Table S2). Worse still, lacking the tying force, the distance CαF186-CαE80 also increased, implying a larger entrance (Figure 8c).
Based on our findings, pima-TE was re-modified into R186Y mutant. This time, we endowed a hydroxyl to the side chain of mutated residue to bond with E80, while the remainder stayed hydrophobic. We were more than pleased to find a significant rise in PRS ratio (4.37% vs. 18.14%), with Y186-E80 bonding partly restored (Table S3). Of particular note, a close-to-reaction PRS conformation appeared in md3 and maintained for over 10 ns, with the terminal methyl oriented towards the entrance and forcing O7 closer to C1 (Figure 8d). We thus regard R186Y as a promising modification towards pimaricin productivity advancement.

Mutation 4-S138C
As presented by Koch et al. [29], compared with pikromycin synthase (PICS)-TEWT, PICS-TES148C could promote macrocyclization efficiency by over 300%. Therefore, pima-TES138C mutant were subjected to MD simulations to study whether the superiority of Cys over Ser applied in pima-TE as well. After the clustering analysis, the dominant polyketide structure of each S138C trajectory demonstrated unbelievable similarity ( Figure S3). As seen in Figure 9, S138C frames were

Mutation 4-S138C
As presented by Koch et al. [29], compared with pikromycin synthase (PICS)-TE WT , PICS-TE S148C could promote macrocyclization efficiency by over 300%. Therefore, pima-TE S138C mutant were subjected to MD simulations to study whether the superiority of Cys over Ser applied in pima-TE as well. After the clustering analysis, the dominant polyketide structure of each S138C trajectory demonstrated unbelievable similarity ( Figure S3). As seen in Figure 9, S138C frames were significantly more concentrated in the conducive range for reaction compared to wild type ones, suggesting potential catalytic advantage. However, due to O→S atomic radius enlargement, bond length involving O/S increased by 0.3 Å, and 0.5 Å, respectively, and distance O 7 -C 1 would mostly gather around 5.5 Å in mutated complex. The density of advantageous conformations in S138C system strongly suggested the favorability of this mutation.

Study on TE's Effect on the Release of Pimaricin Product
The binding energy between pima-TE and polyketide product was calculated with MM/PBSA program (Figure 10a). In study of product (i.e., MOL) movement across the channel, distance between its mass center and Cα S138 was measured (Figure 10b). As a result, the product migrated towards the exit for approximately 4 Å in md1, while it hardly moved in others. Therefore, md1 was regarded to have a tendency of product release, and the other two disclosed the stabilization effect produced by the protein. Energy decomposition revealed residues around the exit to play key parts in protein-product interactions (Figure 10c), and to assume important roles in product release. (Figure S4) significantly more concentrated in the conducive range for reaction compared to wild type ones, suggesting potential catalytic advantage. However, due to O→S atomic radius enlargement, bond length involving O/S increased by 0.3 Å, and 0.5Å, respectively, and distance O7-C1 would mostly gather around 5.5Å in mutated complex. The density of advantageous conformations in S138C system strongly suggested the favorability of this mutation. Figure 9. (a) Diagram of representative PRS conformations in wild type (left) and S138C (right) trajectories. (b) Density map and marginal histogram indicating the distribution of all frames on the basis of distances O7-NεH261 and O7-C1 in wild type and S138C trajectories. The rectangles in lightcoral, slate-blue, and thistle highlight points with distance (O7-NεH261) ≤ 3.0Å, distance (O7-C1) ≤ 4.5Å and 4.5Å ≤ distance (O7-C1) ≤ 6.0Å, respectively.

Study on TE's Effect on the Release of Pimaricin Product
The binding energy between pima-TE and polyketide product was calculated with MM/PBSA program (Figure 10a). In study of product (i.e., MOL) movement across the channel, distance between its mass center and CαS138 was measured (Figure 10b). As a result, the product migrated towards the exit for approximately 4 Å in md1, while it hardly moved in others. Therefore, md1 was regarded to have a tendency of product release, and the other two disclosed the stabilization effect produced by the protein. Energy decomposition revealed residues around the exit to play key parts in proteinproduct interactions (Figure 10c), and to assume important roles in product release. (Figure S4 Next, a careful analysis was conducted on the disengagement of product in md1, and three patterns of hydrogen bonding between Q29 and product were generalized (Figure 11). For the first 25 ns, the product remained its original state and O6 from the product ring continued bonding to Q29 (I). Afterwards, in cooperation with Q29's side chain turnover (II), the ring lied down a little and interacted with Q29 from the right (III). Then, the free hydroxyl on Q29 (OE1) grasped O2 from the Next, a careful analysis was conducted on the disengagement of product in md1, and three patterns of hydrogen bonding between Q29 and product were generalized ( Figure 11). For the first 25 ns, the product remained its original state and O 6 from the product ring continued bonding to Q29 (I). Afterwards, in cooperation with Q29's side chain turnover (II), the ring lied down a little and interacted with Q29 from the right (III). Then, the free hydroxyl on Q29 (OE 1 ) grasped O 2 from the other side of the molecule, further enabling the molecule to lie flat (IV). In a word, steps II-III and III-IV played decisive roles in altering the product's layout and pulling the product farther away from the active site. Due to distinguished distribution of hydrophilic and hydrophobic areas on protein, rotation of the product might partly attenuate its interaction with peripheral residues, and impel the product's departure. On the other hand, H187 seemed to provide thrust towards the release of product ( Figure 11). Distance between the mass center of H187 and product ring shrank along the simulation, revealing established hydrophobic interaction between the imidazole of H187 and terminal methyl on the product ( Figure 11).
Taken together, after cyclization, the product would stay in the vicinity of active site for a while due to van der Waals (VDW) and electrostatic interactions from peripheral residues. Later on, the product layout was altered by molecule rotation, varied hydrogen bonding, etc., which impaired the spatial constraint, and caused the ring to gradually migrate towards the exit, with Q29 hydrogen bonding as a driving force and H187 as a rear helper.

Discussion
Due to the limitation of experimental instruments, present computational strategies combining homology modeling, molecular docking, MD simulation, and QM/MM calculation have been extensively utilized to provide insight into atomistic details in protein-substrate recognition and catalytic mechanism. Over recent years, packages and software [30][31][32] to study protein-substrate interaction have sprung up relentlessly, and MD simulation has become a regular routine herein [33- On the other hand, H187 seemed to provide thrust towards the release of product ( Figure 11). Distance between the mass center of H187 and product ring shrank along the simulation, revealing established hydrophobic interaction between the imidazole of H187 and terminal methyl on the product ( Figure 11).
Taken together, after cyclization, the product would stay in the vicinity of active site for a while due to van der Waals (VDW) and electrostatic interactions from peripheral residues. Later on, the product layout was altered by molecule rotation, varied hydrogen bonding, etc., which impaired the spatial constraint, and caused the ring to gradually migrate towards the exit, with Q29 hydrogen bonding as a driving force and H187 as a rear helper.

Discussion
Due to the limitation of experimental instruments, present computational strategies combining homology modeling, molecular docking, MD simulation, and QM/MM calculation have been extensively utilized to provide insight into atomistic details in protein-substrate recognition and catalytic mechanism. Over recent years, packages and software [30][31][32] to study protein-substrate interaction have sprung up relentlessly, and MD simulation has become a regular routine herein [33][34][35].
In this work, MD simulations were carried out on pima-TE-substrate/product complexes. Residues playing critical roles in product recognition, assembly, and release were uncovered through hydrogen bonding and hydrophobic interaction network analysis, which could be obtained from representative conformations of trajectories, as well as decomposition of MM/PBSA binding energy. Q29 and M210 might contribute to tight binding effect, and the structural correlation between protein and substrate was reduced once they were eliminated. R186 was uncovered to maintain pocket hydrophobicity yet distract the substrate from a proper position, and its mutation to Tyr could benefit macrocyclization by raising the proportion of advantageous conformation. The computer-aided methods could provide theoretical basis to enzyme clarification.
Since the transition states of enzyme catalysis were hard to obtain in silico, we chose pre-reaction state (PRS) as an evaluation indicator. According to our previous research [26,36], PRS was the very prior stage of macrocyclization in terms of both structure and energy, and its formation was decisive to TE cyclization. The proportion of PRS was regarded as the degree of reaction readiness. Besides, PRS proportion was used in mutated systems as well to help elucidate the functions of these residues and speculate their effect on TE activity. However, a more accurate account of the mutation required explanation in energy and experimental verification as well.
To conclude, the study approach applied in our work involved protein-substrate interaction, residue targeting, and mutation analyses with PRS occurrence as an indicator. The strategy could provide structural rationale for TE-substrate complex and guide future experiments on design of efficient protein mutants or novel compounds.

System Preparation
Given the unavailability of pima-TE crystallization data from Protein Data Bank (PDB), initial structure of pima-TE was produced through homology modeling with PICS-TE [37] (PDB: 2H7Y) as a template (sequence similarity: 48.1%). Twenty pima-TE models were generated in the discovery studio 3.5 [38]. The one with the lowest total energy was selected, and its stereochemical quality was further validated by Procheck [39], with 93.7% of its residues falling in the most favored region.
Considering the extensively-acknowledged catalytic process of pimaricin PKS, the mature pimaricin product was disconnected at carbonyl C 1 , and the lactonic ring as well as exocyclic mycosamine were removed. Furthermore, carboxyl on C 12 was also substituted by a methyl. The precursor was optimized with Gaussian09 [40] AM1 method [41], after which the buckled conformation still sustained. The energetically-stabilized substrate was then covalently bonded to active site Ser138 on pima-TE model, with a hydrogen bond forming between its terminal hydroxyl and Nε of active site His261. Protonation state of His261 was altered to HID to facilitate PRS formation. The polyketide-bound acyl-enzyme intermediate was utilized as the initial structure of MD simulations.
During the preparation of the system parameters, an N-terminal cap (-CO-CH 3 ) and a C-terminal cap (-NH-CH 3 ) were firstly added onto the Ser138 to block its ends. Conformational optimization at the level of HF/6-31G(d) was then employed on the intermediate, and its electrostatic surface potential (ESP) charge was computed. Afterwards, a two-step restrained electrostatic potential (RESP) model was applied to determine charge distribution on the substrate. Finally, two prior-added caps were removed, and parameters for the intermediate were generated by the Antechamber package, on the basis of which topology files for protein-substrate complexes were prepared with tleap module in with 12 sodium ions added to maintain charge neutralization.

Molecular Dynamics Simulation
Starting from the solvated polyketide-bound acyl-enzyme intermediate, classical molecular dynamics simulations were carried out utilizing AMBER14 [43] ff03.r1 force field [44]. The system was firstly subjected to 10,000 steps of steepest descent energy minimization followed by 1000 cycles of conjugate gradient minimization with bonds involving hydrogen constrained by SHAKE algorithm [45], and then another 10,000 steps of steepest descent energy minimization followed by 5000 cycles of conjugate gradient minimization with no constraint exerted. The system was then gradually heated from 0 to 300 K through 25000 iterations. After a 200ps-equilibrium in NPT ensemble, five 50-ns simulations (300 K, 1 atm) with different random seeds were conducted. The VDW interactions were cut off at 10 Å and long-range electrostatic interactions were calculated with particle mesh Ewald (PME) method [46]. Analyses of trajectories were performed using cpptraj in Ambertools14.

Quantum Mechanics/Molecular Mechanics) Calculation
Quantum mechanics/molecular mechanics (QM/MM) calculations were performed with a two-layered ONIOM method [47,48] in Gaussian09 program. Geometrical snapshots from the dominant MD cluster were extracted as PRS, active state, and POS, and were further subject to geometry optimization. The quantum mechanical (QM) layer consisted of side chains of active site triad (Ser138, Asp166 and His261) and the polyketide chain, which added up to 96 atoms and bore one negative charge. The optimization process was carried out under M06-2X [49] functional and basis set 6-31G(d) [50].

Simulation of Site Mutation Proteins
Based on the analyses, M210 and Q29 were selected as key residues in protein-substrate interaction. Considering the optimization of pima-TE, R186 was mutated to Phe and Tyr in succession to reduce its interference against PRS formation. In accordance with a previously published article of Koch et al. [29], S138 was also mutated to Cys to examine pima-TE S138C 's effectiveness. Single site mutation was employed directly on the initial structure of wild type pima-TE, and all mutants (M210G, Q29A, R186F, R186Y, S138C) went through 30-50 ns simulations following identical procedures as mentioned in Section 4.2.

Free Energy Calculation and Conformational Stability Analysis
The polyketide chain, which was extracted from dominant structure in wild type md4, was manually rang up and subjected to conformational optimization with Gaussian 09 AM1 method. The optimized product was then docked into the channel with C 1 adjacent to the active site. After the model construction, 3 × 50 ns MD simulations was carried out with AMBER14 program.
After clustering analysis in cpptraj, a 20 ns segment with dominant conformation was extracted from each trajectory, and was further subject to a molecular mechanics Poisson-Boltzmann surface area [51] (MM/PBSA) calculation to estimate the free energy difference (∆G tot ) between bound and detached states of product-protein complexes in solution. The MMPBSA.py program in AMBER14 was performed, and the free energy discrepancy was decomposed to peripheral residues in terms of hydrophobic and electrostatic forces. Table 1 lists the number and duration of all MD simulations utilized in the study. pima-TE R186Y + polyketide chain 3 30 S138C pima-TE S138C + polyketide chain 3 50 Product ring pima-TE WT + product 3 50

Conclusions
In this paper, MD simulations were utilized as a primary tool to explore pimaricin TE catalysis on an atomic level. Firstly, 5 × 50 ns trajectories on polyketide were conducted in search of pre-reaction states (PRS), and transformation between POS and PRS were examined. POS was found to bear lower energy, yet less mature conformation in comparison with PRS. Protein-polyketide hydrogen bonding and hydrophobic interactions were deciphered, with several key residues subjected to mutations. As discovered, Q29 was responsible for holding a polyketide hydroxyl and controlling the substrate position, and M210 contributed to favorable protein-ligand interaction by virtue of its hydrophobicity. R186Y might promote productivity by reducing the interference on PRS formation, and S138C could effectively enhance the proportion of required conformations. Ultimately, the MM/PBSA program was employed to unveil residues mediating product release, and the postulation of a mechanism of polyketide product departure from the active site was proposed. We gave a comprehensive overview on pima-TE catalysis, with computational methods, and offered opinions for protein engineering.