**Shuobing Fan, Rufan Wang , Chen Li, Linquan Bai, Yi-Lei Zhao and Ting Shi \***

State Key Laboratory of Microbial Metabolism, Joint International Research Laboratory of Metabolic and Developmental Sciences, School of Life Sciences and Biotechnology, Shanghai Jiao Tong University, Shanghai 200240, China; bigyellowjar@sjtu.edu.cn (S.F.); wangrufan@sjtu.edu.cn (R.W.); l1105309790@sjtu.edu.cn (C.L.); bailq@sjtu.edu.cn (L.B.); yileizhao@sjtu.edu.cn (Y.-L.Z.)

**\*** Correspondence: tshi@sjtu.edu.cn; Tel./Fax: +86-21-34207347

Received: 21 January 2019; Accepted: 6 February 2019; Published: 18 February 2019

**Abstract:** 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.

**Keywords:** pimaricin thioesterase; protein-substrate interaction; macrocyclization; molecular dynamics (MD) simulation; pre-reaction state

#### **1. Introduction**

Produced predominantly by the genus *Streptomyces*, polyene polyketides consist of a large family of natural compounds [1,2], including pimaricin [3], amphotericin B [4], nystatin A1 [5], and candicidin/FR008 [6]. The members of this class are widely used in clinical medicine for their broad spectral antifungal properties [7]. With constant progress in scientific research, their potential pharmaceutical values on antiviral, antiprotozoal, antiprion, and anticancer activity have been progressively clarified [8,9].

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–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]. makes pimaricin appealing to researchers, and its biosynthetic pathway modification and drug design have become new science hotspots [17].

*Int. J. Mol. Sci.* **2019**, *20*, x FOR PEER REVIEW 2 of 18

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–21]. Situated in the last module, the thioesterase (TE) domain off-loads the ACP-tethered polyketide from PKS via macrocyclization or hydrolysis. 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–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 acyl-enzyme intermediate with no efficient intra-molecular nucleophile present. 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 O<sup>11</sup> and carbonyl oxygen O (Figure 1), as accompanied by the swing of C<sup>11</sup> 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 O<sup>11</sup> 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-C<sup>1</sup> for nucleophilic attack was larger than 6 Å in *active state*. Finally, an advantageous conformer (the pre-reaction state, PRS) was found [24] after ~270 ns MD simulation, which possessed both hydrogen bond O11-NεH259 and an appropriate distance between O<sup>11</sup> and C<sup>l</sup> to facilitate nucleophilic attack. Critical to mutual recognition and catalytic preparation between TE and substrate, the PRS seemed decisive in the occurrence of macrocyclization. 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.

**Figure 1.** Structures of pre-organization state (POS), active state and the pre-reaction state (PRS) of 6- **Figure 1.** Structures of pre-organization state (POS), active state and the pre-reaction state (PRS) of 6-deoxyerythronolide B synthase (DEBS) thioesterase (TE) system.

deoxyerythronolide B synthase (DEBS) thioesterase (TE) system.

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

critical residues highlighted. We also provided an explanation on the departure of product from the

#### **2. Results and Discussion 2. Results and Discussion**

#### *2.1. Key Structural Conformations in MD Simulations 2.1. 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*. 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*.

**Figure 2.** Conformational change of pima-TE system during molecular dynamics (MD) simulations. (**a**) Structural variations between post (opaque) and pre-simulation (transparent) complexes, with lid region, polyketide chain, α-helix αL2, αL3 and loop *l1*, *l2* & *l3* colored in tv\_blue, gray, yellow, cyan, green, red and orange. (**b**) Root-mean-squared fluctuation (RMSF) of five trajectories with key structural elements highlighted. **Figure 2.** Conformational change of pima-TE system during molecular dynamics (MD) simulations. (**a**) Structural variations between post (opaque) and pre-simulation (transparent) complexes, with lid region, polyketide chain, α-helix αL2, αL3 and loop *l1*, *l2* & *l3* colored in tv\_blue, gray, yellow, cyan, green, red and orange. (**b**) Root-mean-squared fluctuation (RMSF) of five trajectories with key structural elements highlighted.

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<sup>7</sup> 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* 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<sup>7</sup> was proposed to be conducive for nucleophilic attack onto carboxyl C<sup>1</sup> 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. 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 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.

*Int. J. Mol. Sci.* **2019**, *20*, x FOR PEER REVIEW 4 of 18

**Figure 3.** Classification of trajectory frames based on polyketide chain conformation. (**a**) Representative structures of PRS, *active state* and POS extracted from clustering analysis. (**b**) Proportion of PRS, *active state*, and POS in each trajectory. **Figure 3.** Classification of trajectory frames based on polyketide chain conformation. (**a**) Representative structures of PRS, *active state* and POS extracted from clustering analysis. (**b**) Proportion of PRS, *active state*, and POS in each trajectory.

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). 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<sup>7</sup> and carbonyl oxygen O<sup>1</sup> 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).

#### *2.2. Conformational Transition between POS and PRS 2.2. Conformational Transition between POS and PRS*

Next, the transformation between POS and PRS was studied using dihedral angle Cα-Cβ-Cγ-O7 as an indicator of polyketide terminal rotation. In PRS, bond Cα-Cβ ran anti-parallel against Cγ-O7 Next, the transformation between POS and PRS was studied using dihedral angle Cα-Cβ-Cγ-O<sup>7</sup> as an indicator of polyketide terminal rotation. In PRS, bond Cα-C<sup>β</sup> ran anti-parallel against Cγ-O<sup>7</sup> (−180◦ ), but in POS, the dihedral angle was altered to an acute angle fluctuating between (−30◦ , −70◦ ).

(−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 O7 firstly swung up and disassociated from H261, followed by Cβ-Cγ twisting clockwise and terminal methyl 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<sup>7</sup> firstly swung up and disassociated from H261, followed by Cβ-C<sup>γ</sup> 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). *Int. J. Mol. Sci.* **2019**, *20*, x FOR PEER REVIEW 5 of 18

**Figure 4.** Conformational transformation between PRS and POS. (**a**) Dihedral value representation of md1 along with intermediate conformation change. (**b**) Presentation of dihedral angle Cα-Cβ-Cγ-O7. **Figure 4.** Conformational transformation between PRS and POS. (**a**) Dihedral value representation of md1 along with intermediate conformation change. (**b**) Presentation of dihedral angle Cα-Cβ-Cγ-O<sup>7</sup> .

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 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-O<sup>1</sup> 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<sup>γ</sup> bond rotation.

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

through dihedral flip and conformation adjustment, and the energies on POS, *active state,* and PRS

#### were computed to understand the reaction process. *2.3. Hydrophilic and Hydrophobic Interactions in Pima-TE System*

*2.3. 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 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, O<sup>4</sup> or O<sup>5</sup> 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<sup>6</sup> and gave rise to an erected molecule (39.26%). The main chain of M210 fixed O<sup>4</sup> 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.

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. 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%).

**Figure 5.** Hydrogen-bonding and hydrophobic interaction network variations in PRS and POS. (**a**) Proportion of top-ranked hydrogen bonding interactions. (**b**) Diagram of protein-substrate interaction produced by ligplot+ [28]. Backbone of the polyketide substrate was colored in yellow, and residues providing hydrogen bonding and hydrophobic interactions in slate\_blue and brown. **Figure 5.** Hydrogen-bonding and hydrophobic interaction network variations in PRS and POS. (**a**) Proportion of top-ranked hydrogen bonding interactions. (**b**) Diagram of protein-substrate interaction produced by ligplot<sup>+</sup> [28]. Backbone of the polyketide substrate was colored in yellow, and residues providing hydrogen bonding and hydrophobic interactions in slate\_blue and brown.

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

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 *2.4. Key Residues Analyzed Via Mutant Simulations*

#### driving force (99.28%). Besides Y180 (36.17%) and L183 (57.06%) as in PRS, L213 also participated in hydrophobicity maintenance (25.43%). 2.4.1. Mutation 1-Q29A

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. *2.4. Key Residues Analyzed Via Mutant Simulations*  2.4.1. 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 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 O<sup>6</sup> 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<sup>6</sup> 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 (**2<sup>I</sup>** , **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% (**1<sup>I</sup>** , **1II**).

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**).

**Figure 6.** Instability of polyketide conformation. (**a**) Distance O6-CQ/A29 in wild type md4 and Q29A **Figure 6.** Instability of polyketide conformation. (**a**) Distance O<sup>6</sup> -CQ/A29 in wild type md4 and Q29A md1–3 along with conformation transformation. (**b**) Diagram of distance O<sup>6</sup> -CQ/A29.

md1–3 along with conformation transformation. (**b**) Diagram of distance O6-CQ/A29. From our perspective, Q29 could regulate hydrogen bond O7-NεH261 and PRS formation by From our perspective, Q29 could regulate hydrogen bond O7-NεH261 and PRS formation by binding the substrate position with a hydrogen bond, while having little effect on POS.

#### binding the substrate position with a hydrogen bond, while having little effect on POS. 2.4.2. Mutation 2-M210G

an advantageous channel exit shape.

2.4.2. 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, 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%).

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

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

**Figure 7.** Conformational change of intermediate change upon M210 mutation. (**a**) Distance RMSD (dRMSD) value of substrate backbone in wild type md4 and M210G (lightpink for wild type, **Figure 7.** Conformational change of intermediate change upon M210 mutation. (**a**) Distance RMSD (dRMSD) value of substrate backbone in wild type md4 and M210G (lightpink for wild type, lightblue, slate and tv\_blue for M210G md1–3). (**b**) Diagram of the dominant structure in each trajectory.

#### lightblue, slate and tv\_blue for M210G md1–3). (**b**) Diagram of the dominant structure in each 2.4.3. Mutation 3-R186F & R186Y

trajectory. 2.4.3. Mutation 3-R186F & R186Y As a residue containing multiple hydrophilic groups, R186 bonded with O7 for a rather high As a residue containing multiple hydrophilic groups, R186 bonded with O<sup>7</sup> 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.

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. 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<sup>7</sup> closer to C<sup>1</sup> (Figure 8d). We thus regard R186Y as a promising modification towards pimaricin productivity advancement.

**Figure 8.** Mutational trials on R186. (**a**) Radar chart indicating the proportion of hydrogen bond O7- NεH261 and O7-NR186 formation within 5 wild type simulations. (**b**) Larger entrance of pima-TE after R186F mutation. (**c**) Coupling and non-coupling states of three entrance residues (R/F186, E80 and **Figure 8.** Mutational trials on R186. (**a**) Radar chart indicating the proportion of hydrogen bond O<sup>7</sup> -NεH261 and O<sup>7</sup> -NR186 formation within 5 wild type simulations. (**b**) Larger entrance of pima-TE after R186F mutation. (**c**) Coupling and non-coupling states of three entrance residues (R/F186, E80 and R266) located on different structure elements. (**d**) A favorable PRS emerged in R186Y md3.

#### R266) located on different structure elements. (**d**) A favorable PRS emerged in R186Y md3. 2.4.4. Mutation 4-S138C

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 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 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-C<sup>1</sup> would mostly gather around 5.5 Å in mutated complex. The density of advantageous conformations in S138C system strongly suggested the favorability of this mutation.

#### towards the entrance and forcing O7 closer to C1 (Figure 8d). We thus regard R186Y as a promising *2.5. Study on TE's Effect on the Release of Pimaricin Product*

modification towards pimaricin productivity advancement. 2.4.4. 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 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)

conformation appeared in md3 and maintained for over 10 ns, with the terminal methyl oriented

demonstrated unbelievable similarity (Figure S3). As seen in Figure 9, S138C frames were

product's departure.

system strongly suggested the favorability of this mutation.

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

**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. **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 O<sup>7</sup> -NεH261 and O<sup>7</sup> -C<sup>1</sup> in wild type and S138C trajectories. The rectangles in light-coral, slate-blue, and thistle highlight points with distance (O<sup>7</sup> -NεH261) ≤ 3.0 Å, distance (O<sup>7</sup> -C<sup>1</sup> ) ≤ 4.5 Å and 4.5 Å ≤ distance (O<sup>7</sup> -C<sup>1</sup> ) ≤ 6.0 Å, respectively. *Int. J. Mol. Sci.* **2019**, *20*, x FOR PEER REVIEW 11 of 18

**Figure 10.** MM/PBSA analysis on the outward trend of macro-lactone product. (**a**) Binding energy between pima-TE and product ring. (**b**) Distance between ring mass center and CαS138. (**c**) Residues with top-ranked van der Waals (VDW) and electrostatic contributions to binding free energy. **Figure 10.** MM/PBSA analysis on the outward trend of macro-lactone product. (**a**) Binding energy between pima-TE and product ring. (**b**) Distance between ring mass center and CαS138. (**c**) Residues with top-ranked van der Waals (VDW) and electrostatic contributions to binding free energy.

interacted with Q29 from the right (III**)**. Then, the free hydroxyl on Q29 (OE1) grasped O2 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

Next, a careful analysis was conducted on the disengagement of product in md1, and three

**3. Discussion** 

35].

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<sup>6</sup> 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 O<sup>2</sup> 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. *Int. J. Mol. Sci.* **2019**, *20*, x FOR PEER REVIEW 12 of 18

**Figure 11.** Diagram of product movement in substrate channel in md1. (**a**) Three hydrogen bonding modes between Q29 and MOL and their transformation. (**b**) Spatial location of H187 with respect to MOL. (**c**) Distance H187-MOL in md1. **Figure 11.** Diagram of product movement in substrate channel in md1. (**a**) Three hydrogen bonding modes between Q29 and MOL and their transformation. (**b**) Spatial location of H187 with respect to MOL. (**c**) Distance H187-MOL in md1.

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). 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. 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.

catalytic mechanism. Over recent years, packages and software [30–32] to study protein-substrate interaction have sprung up relentlessly, and MD simulation has become a regular routine herein [33–

Due to the limitation of experimental instruments, present computational strategies combining

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

#### **3. 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–32] to study protein- substrate interaction have sprung up relentlessly, and MD simulation has become a regular routine herein [33–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.

#### **4. Materials and Methods**

#### *4.1. 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 C1, and the lactonic ring as well as exocyclic mycosamine were removed. Furthermore, carboxyl on C<sup>12</sup> 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-CH3) and a C-terminal cap (-NH-CH3) 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 AMBER 14. Through *tleap*, pima TE-substrate system was placed in an octahedral TIP3P water box [42], with 12 sodium ions added to maintain charge neutralization.

#### *4.2. 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.

#### *4.3. 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].

#### *4.4. 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-TES138C'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.

#### *4.5. 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<sup>1</sup> 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.



#### **5. 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.

**Supplementary Materials:** Supplementary materials can be found at http://www.mdpi.com/1422-0067/20/4/ 877/s1.

**Author Contributions:** Conceptualization, L.B., Y.-L.Z., and T.S.; formal analysis, S.F., R.W., and C.L.; funding acquisition, L.B. and Y.-L.Z.; investigation, S.F., R.W., and C.L.; methodology, Y.-L.Z. and T.S.; writing-original draft, S.F. and T.S.; writing-review and editing, S.F. and T.S.

**Funding:** This research was funded by the National Science Foundation of China, grant number 21377085 and 31770070.

**Acknowledgments:** The authors thank SJTU-YG2016MS42, SJTU-ZH2018ZDA26, SJTU-YG2016MS33 for financial supports and SJTU-HPC computing facility award for computational hours.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**


#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
