Modeling and analysis of substrate bound holo-ACP:DH complex
The modeling of substrate bound holo-ACP:DH complex was carried out in several steps involving protein-protein docking of apo-ACP:DH complexes, docking of P-pant group, docking of substrate moiety and finally formation of required covalent bonds and energy minimization/refinement of the substrate bound complex as described in methods section (Figure

Additional file
1: Figures S1 and Figure S2). During docking of apo-ACP on DH domain, 10,000 complexes obtained from FTDOCK program [
26] were re-ranked using residue-residue pair potentials, RPscore [
27]. Out of these 10000 complexes, 355 complexes with a positive RPscore were selected for further analysis. These 355 complexes were further filtered using functional constraints. The functional requirement for catalysis by DH domain involved P-pant group attached to the Ser of ACP to reach the active site pocket of DH domain, thus bringing Ser in proximity with the residues which are located near the entrance to the active site cavity of the DH domain. Therefore, a distance cutoff of 5

Å between any atom of catalytic Ser of ACP and any atom of Phe 227 and Arg 275 of DH domain was used as functional constraint. In fact, site directed mutagenesis studies reported earlier have suggested involvement of these two residues from DH domain in recognition of ACP [
19]. Interestingly, R275D mutation in erythromycin DH is known to significantly lower the yield of 6-dEB and based on these results, role of Arg 275 in ACP recognition has been proposed [
19]. These functional constraints based filtering resulted in only 10 complexes satisfying the criteria. Detailed analysis of interface residues in these 10 apo-ACP:DH complexes indicated that, complex 2 (shown in Additional file
1: Figure S3) has not only high RPScore value, but also optimal orientation of Ser 46 for extension of the covalently attached P-pant group into the DH binding pocket. It also had most favorable interactions involving Asp 45 (ACP): Arg 275 (DH) and Leu 47 (ACP): Phe 227 (DH) pairs. Hence complex 2 was selected for the next step of P-pant docking to obtain holo-ACP:DH complex. The top panel of Figure

depicts the 10 solutions obtained from protein-protein docking of DH and apo-ACP domain. Complex 2 is represented with helices of ACP domain in blue color while in other solutions ACP domain is shown in orange color. The Ser of ACP is placed at the entrance of a deep cavity which leads to the catalytic His 44 and Asp 206 of DH domain [
19].
After obtaining the biologically meaningful structural model for apo-ACP:DH complex, the next step was to model the holo-ACP:DH complex by docking of P-pant moiety as described in the methods section. The total of 250 bound conformations of P-pant obtained from AutoDock were clustered using a RMSD cut off of 2

Å and this resulted in 113 clusters with different binding energy values. The central panel in Additional file
1: Figure S4 shows the binding energy values and number of conformations for these 113 clusters obtained from P-pant docking. Since, functionally meaningful P-pant conformation should form covalent bond between Ser of ACP and its thiol group should be close to the catalytic residues of DH domain, the cluster corresponding to the lowest energy and highest population might not always be biologically meaningful. Therefore, distance between Oγ atom of Ser 46 (ACP) and phosphate of P-pant as well as distance between –SH group of P-pant and Oδ atom of Asp 206 (DH) was computed for representative conformations from each of these 113 clusters. The representative conformations from each of the clusters having the values of above-mentioned distances less than 7

Å but more than 3.5

Å have been depicted in middle panel of Figure

, while the corresponding values of distances as well as binding energy are shown in the bottom panel. The cases where the distance was less than 3.5

Å between Ser of ACP and phosphate group of P-pant showed a steric clash between them (Clusters 3 and 8 shown in lower panel of Figure

). After removal of such solutions, the cluster 5 was chosen in terms of minimum binding energy as well as distance constraint required for formation of covalent bond. As can be seen from the middle and lower panels of Figure

, all the solutions satisfying the functional constraints fall into the same deep cavity leading to the catalytic residues and binding energy of these solutions varies in the range of −3.19 to −5.24

kcal/mol. The Additional file
1: Figure S4 shows these two distances for all clusters (upper panel) and the relative orientation of the conformation of P-pant group from cluster 5 with respect to Ser 46 (ACP) and Asp 206 (DH) (lower panel). Therefore, docked conformation corresponding to cluster 5 was chosen for modeling of holo-ACP:DH complex by formation of required covalent bond between Ser 46 and P-pant moiety. The middle panel of Figure

shows in stick representation the covalently attached P-pant conformation (larger thickness and deep purple color) obtained after minimization of the complete complex containing DH, ACP and P-pant moiety from cluster 5.
The next requirement was to dock the substrate moiety onto holo-ACP:DH complex such that the carboxyl group of the substrate will be at a position to form covalent bond with the S atom in P-pant and also the beta hydroxyl group of the substrate will be in close proximity of the catalytic residues (Asp 206 and His 44) in the DH domains. As described earlier in the methods section, the substrate for DH domain was docked as two separate fragments ( Additional file
1: Figure S1) and the required covalent bonds were made after selecting suitable solutions from these two AutoDock runs. The docking of the first fragment of the substrate ( Additional file
1: Figure S5) yielded 250 conformations belonging to 32 clusters (RMSD within cluster <2

Å) with different binding energy values. Additional file
1: Figure S5a shows the number of conformations, distance between Oδ atom of Asp 206 (DH) and beta hydroxyl of the substrate and also the distance between the carboxyl carbon of the substrate fragment and S atom of P-pant for each cluster. Inset to Additional file
1: Figure S5a shows the conformation of the substrate fragment 1 from the cluster 1 showing minimum values of these two distances also. Additional file
1: Figure S5b shows similar results from docking of the second fragment and distribution of distance between Cγ-Cδ in various docked clusters. The whole complex, hence obtained, was minimized after building the required bonds between ACP, P-pant and substrate fragments to obtain the final substrate bound holo-ACP:DH complex. Additional file
1: Figure S6a shows the final energy minimized structure for the substrate bound holo-ACP:DH complex. As can be seen, the substrate moiety as well as a portion of the P-pant group enters into a deep cavity which harbors the catalytic residues of the DH domain. In order to further validate the results from the current study, the bound conformation of the substrate obtained from the current study was compared with the bound conformation of the mechanistic inhibitor 3-decynoyl-N-acetylcysteamine which has been crystallized in complex with FabA (PDB ID 1MKA), the type II FAS DH enzyme from E. coli [
9]. Additional file
1: Figure S6b shows the superposition of 1MKA on the DH domain in substrate bound holo-ACP:DH complex obtained from the current study. As can be seen, in both cases the ligands bind into the narrow tunnel of DH leading to the catalytic residues ( Additional file
1: Figure S6b) and the bound conformation of the mechanism-based inhibitor 3-decynoyl-N-acetylcysteamine is also very similar to the conformation of the bound substrate in the ACP:DH complex. Additional file
1: Figure S6b shows the complete docked substrate backbone in green while the backbone atoms of the mechanism-based inhibitor have been depicted in blue color. The only difference in the conformation lies at the α-β positions because the inhibitor covalently binds to the catalytic His residue with cis conformation of the α-β bond[
9], while the substrate for EryDH4 binds non-covalently in trans conformation of the corresponding α-β bond ( Additional file
1: Figure S6b). It shows that the binding site obtained by docking results is same as that observed by earlier experimental studies.
Detailed analysis of the substrate bound holo-ACP:DH complex indicated that the charged residues in DH domains interact with hydroxyl moieties of the substrate while the hydrophobic residues were in contact with methyl groups. The P-pant moiety was exposed to solvent to a large extent and therefore, fewer charged or polar residues had interaction with the P-pant moiety in the apo-ACP:DH complex. Table

shows the list of residues in DH domain, which interact with the ACP, P-pant and the substrate moiety. We also wanted to investigate whether the residues of DH domain which interact with ACP, P-pant and substrate moieties remain conserved in other modular PKS clusters. Therefore, the evolutionary conservation of these residues was checked. Additional file
1: Figure S7 shows the multiple sequence alignment (MSA) of DH domains from different modular PKS cluster with the interacting residues marked using arrow symbol. As can be seen from the MSA in Additional file
1: Figure S7, 10 out of 13 residues (His44, Pro53, Gly54, Ser55, Leu161, Tyr168, Asp206, Ala209, Pro226, Val274) interacting with substrate, 9 out of 13 residues (His44, Leu51, Val52, Gly54, Pro81, Leu82, Leu225, Pro226, Phe227) interacting with phosphopantetheine and both residues (Phe227 and Arg275) interacting with ACP show good conservation. The conservation of these residues indicates that the interactions of DH domain with ACP, P-pant and substrate obtained in holo-ACP:DH complex modeled for eryPKS can also be used to predict the interacting residues in case of DH domains from other clusters whose structures are not available yet.
| Table 1Contacts of DH domain with bound substrate and ACP in holo-ACP:DH complex |
Since, the protein-protein as well as protein-ligand docking used to model this substrate bound holo-ACP:DH complex did not incorporate flexibility of protein residues, MD simulations were carried out in the explicit solvent environment for a period of 50

ns to further refine the substrate bound holo-ACP:DH complex. This helped in incorporating protein as well as ligand flexibility and also to check the stability of the solutions obtained by docking studies. Detailed structural analysis of the conformations of the ACP:DH complexes sampled during the 50

ns indicated that the entire complex showed a backbone RMSD in the range of 3.0 to 4

Å with respect to the starting structure used for MD simulations (Figure

a). This suggested that there were no large scale changes in the relative orientation of ACP domain with respect to the DH domain and the structure obtained from docking was stable. Analysis of the P-pant attached substrate conformations over the entire 50

ns trajectory also indicated that Oδ atom of Asp 206 (DH) remained within 4.0-5.0

Å of the beta-hydroxyl group of the substrate, while the distance between alpha carbon of the substrate and Nϵ of His was within 6.0-7.0

Å (Figure

b). This suggests that bond between Cα and Cβ which is dehydrated to a trans double bond during catalysis by DH indeed remains in close proximity of the catalytic residues of the DH domain throughout the 50

ns simulation. Figure

shows the stereo view of the final complex obtained after 50

ns MD simulation and the contacting residues with P-pant and substrate moiety are depicted in orange color. Additional file
1: Figure S8 shows the comparison of the substrate binding tunnel in the crystal structure of the substrate free DH domain ( Additional file
1: Figure S8a) with the tunnel in the substrate bound DH domain ( Additional file
1: Figure S8b) after removal of the P-pant bound substrate moiety. As can be seen, there are subtle changes in the structure of the DH domain which facilitates substrate binding by opening of the tunnel. In fact, similar substrate binding tunnel is also seen in other DH domains. Akey
et al[
18] have suggested that substrate binding to curacin DH would also involve such opening of tunnel as was seen in the current docking and MD simulation studies on erythromycin DH domain.
Detailed analysis of the 50

ns MD trajectory indicated that the P-pant attached substrate moiety indeed sampled multiple conformations within the active site tunnel of the DH domain. Therefore, it was interesting to investigate whether the binding pocket residues identified by docking persisted throughout the 50

ns simulations and if new contacts were formed during the MD simulation. Figure

a shows the percentage of time for which the various residues of DH domain remained in contact with the P-pant as well as substrate group during the 50

ns simulation. Thus, different residues interact with P-pant and the substrate moieties for different durations. It is interesting to note that, only 26 residues showed interaction with P-pant or substrate for more than 80% of the simulation time. On the other hand, out of 18 residues of DH domain showing contact for less than 60% of the time, 14 were involved in contacts with the substrate while only 4 residues showed contacts with P-pant moiety. The residues showing higher persistence times are likely to be more critical for substrate and P-pant binding. These results have interesting implications for identifying substrate binding residues of DH domains by experimental studies. Figure

b shows the snapshots of conformations taken up by P-pant and substrate moiety extracted from the 50

ns trajectory at an interval of 1

ns. The residues interacting with the P-pant-substrate moiety for less than 40% but more than 10% of the time have also been depicted. These residues might be interacting with only a small population of the conformations of the P-pant and substrate moiety. Thus incorporation of flexibility by MD simulations has given a dynamic picture of the interactions between the substrate and DH domain.
Modeling and analysis of substrate bound holo-ACP:KR complex
The modeling of substrate bound holo-ACP:KR complex was carried out using a similar protocol (Figure

) as that used for ACP:DH complex involving modeling of apo-ACP:KR complex by protein-protein docking, generation of holo-ACP:KR complex by docking of P-pant group and finally docking of the substrate group. However, the substrate moiety was docked in a single step because of its smaller size unlike the case of docking of substrate on DH domain.
Docking of apo-ACP and KR domains by FTDOCK yielded 10000 complexes which were re-ranked using residue-residue pair potentials (RPscore). The set of 267 complexes obtained with a positive RPscore were further analyzed to identify the complexes satisfying functional constraints. The functional requirement for catalysis by KR domain necessitates the 20

Å long P-pant group attached to the Ser of ACP to reach the active site pocket of KR domain. Therefore, the complexes were filtered with a distance constraint of 20

Å between catalytic Ser 46 of ACP and the catalytic Tyr 1813 of KR domain. It may be noted that, in case of apo-ACP:DH docking study, distance constraints between Ser 46 of ACP and residues on DH domain were used because experimental studies had suggested involvement of corresponding residues in ACP recognition. However, no such information was available about the interacting residues involved in binding of ACP to the KR domain. Therefore, only the 20

Å distance constraint between Cα atom of Ser 46 of ACP and Cα atom of catalytic Tyr 1813 of KR was used to filter the complexes having positive RPscore. This functional constraint yielded six complexes that satisfied the criteria and the best scoring complex was chosen as the final apo-ACP:KR complex. The top panel in Figure

shows the six complexes satisfying functional constraints obtained from apo-ACP:KR docking, while Additional file
1: Figure S9a shows the final apo-ACP:KR complex selected for P-pant docking to obtain holo-ACP:KR complex. The analysis of interface residues of the apo-ACP:KR complex obtained from this docking study revealed favorable interactions involving Arg 1857 (KR): Asp 45 (ACP) and Phe 1856 (KR): Leu 47 (ACP) pairs ( Additional file
1: Figure S9b). It is interesting to note that interacting residue pairs in ACP:KR interface are also identical to the interacting pairs in ACP:DH complex, as discussed earlier. Similar interaction interface involving Asp:Arg pair is also seen in complexes of ACP with ACP synthetase (AcpS) [
7,
28]. This suggests that the functional constraint based docking study has successfully identified interaction interface for apo-ACP:KR complex. Analysis of apo-ACP:KR complex also revealed that Ser 46 of ACP was placed at the entrance of the shallow solvent exposed cavity on the surface of the KR domain while the opposite end of this cavity harbored the NADPH binding site adjacent to the catalytic residues.
In order to model the holo-ACP:KR complex, the P-pant moiety was docked on the final apo-ACP:KR complex obtained from protein-protein docking. Docking of P-pant moiety by AutoDock on the apo-ACP:KR complex resulted in 250 bound conformations belonging to 150 different clusters (RMSD cutoff of 2

Å). The lower left panel in Additional file
1: Figure S10 shows the number of conformations as well as binding energy values for each of these 150 clusters, while the top panel shows the distance between phosphate of P-pant and Ser 46 of ACP as well as the distance of the thiol group of P-pant from the KR bound NADPH for each of these clusters. Since, functionally meaningful P-pant conformation should form covalent bond between Ser 46 of ACP and its thiol group should be close to the catalytic residues of KR domain, the cluster corresponding to the lowest energy and highest population need not be biologically meaningful. Therefore, distance between Oγ atom of Ser 46 of ACP and phosphate of P-pant as well as distance between –SH group of P-pant and NADPH was used as functional constraints for filtering functionally meaningful bound conformation for P-pant group. The solutions with the values of above mentioned distances less than 8

Å but more than 3.5

Å have been depicted in middle panel of Figure

while the corresponding values of distances as well as binding energy are shown in the bottom panel. As can be seen from Figure
, all these solutions fall into the same shallow cavity leading to the catalytic residues and binding energy of these solutions varies in the range of −1.75 to −6.13

kcal/mol. Even though cluster 8 was best solution in terms of minimum binding energy as well as distances, it had a steric clash between P-pant group and catalytic Tyrosine. Thus, the next best solution cluster 14 was chosen for modeling the holo-ACP:KR complex by formation of required covalent bond between Ser 46 and P-pant moiety. The middle panel of Figure

shows the stick representation (larger thickness and cyan color) of the P-pant conformation after minimization of the complete complex containing KR, ACP and P-pant moiety from cluster 14. Although, in the docked conformation of the P-pant moiety, the distance between NADPH and thiol group of P-pant was around 7

Å ( Additional file
1: Figure S10), after minimization the distance was reduced to 3.6

Å. In the holo-ACP:KR complex, the P-pant group was located in the shallow cavity which also contained NADPH binding site.
In order to build the substrate bound holo-ACP:KR complex, the diketide substrate moiety was docked onto the holo-ACP:KR complex subject to the constraint that, the terminal carboxyl carbon of the substrate will be covalently bonded to the S atom in P-pant and also the beta keto group of the substrate will be in close proximity of the catalytic residues in the KR domain. The various bound conformations of the substrate obtained from AutoDock were grouped into 15 clusters (RMSD within a cluster <2

Å) with different binding energy values. Additional file
1: Figure S11 shows the number of conformations as well as distances corresponding to the functional constraints in each of the 15 clusters. The inset in Additional file
1: Figure S11 shows the docked conformation of the substrate from cluster 12 which had minimum values for the two distance constraints. The lowest energy docked conformation from cluster 12 was selected for modeling the substrate bound holo-ACP:KR complex by formation of required covalent bond between P-pant and substrate. The whole complex, hence obtained, was energy minimized and the final energy minimized substrate bound holo-ACP:KR complex is shown in Additional file
1: Figure S12. As can be seen, the substrate moiety and the P-pant group binds into ashallow cavity exposed to the surface and span the region from entrance of the cavity till the NADPH binding site, which is located at the other end of this cavity. Figure

shows the conformation of the P-pant bound substrate and its orientation with respect to the key stereo specificity determining residues in the catalytic pocket of the KR domain. In addition to the catalytic residues, it has also been proposed that LDD motif interacts with substrate and is responsible for determining stereo-specificity in case of B-type KRs [
20]. In addition, the second aspartate of the LDD motif as well as the Phe1801 residue is proposed to determine stereo specificity of keto-reduction by experimental studies. In contrast, the conserved Tryptophan is responsible for guiding the substrate into the binding cavity in case of A-type KRs. In case of B-type KRs the conserved Tryptophan is replaced by a Phenylalanine residue (F1805). It has been proposed that, F1805 blocks the substrate binding in orientations similar to A-type KRs and hence the substrate interacts with LDD motif and enters from other side of the cavity. As can be seen from Figure

a, in the substrate conformation obtained from our docking studies, the Leu of this LDD motif (Leucine 1756 depicted in red) is in close proximity of the hydrophobic region of the P-pant moiety and is away from F1805. Keatinge-Clay & Stroud have also proposed that Tyr 1813 is liberated out of the helix αF by the helix-breaking residue Pro 1815. This imparts higher mobility to Tyr 1813 which allows it to abstract hydrogen from alpha position on the polyketide to bring about epimerization [
20] at alpha-position. After epimerization Val 1852 from lid helix as well as Leu 1810 can interact with epimerized methyl group via hydrophobic interactions and hence, stabilize it after the KR domain brings about the epimerization (Figure

b) [
20]. As can be seen from Figure

b, Tyr 1813 and Val 1852 are located on either side of the substrate. Thus, the substrate binding site identified by the current docking and MD simulations is also consistent with the mechanistic details of the proposed mechanism of epimerization by KR domain at alpha carbon of the substrate.
Analysis of the apo-ACP:KR complex indicated that the charged residues in KR domains formed contacts with hydroxyl or carboxyl moiety of the substrate while the hydrophobic residues were observed to interact with methyl groups. In case of P-pant moiety, the charged residues in KR domain interact with the charged phosphate group as well as carboxyl atoms while the non-polar residues interact with the methyl groups of P-pant. Table

shows the corresponding interacting residues and they have been marked on the multiple sequence alignment (MSA) of KR domains to observe their evolutionary conservation ( Additional file
1: Figure S13). As can be seen, 5 out of 8 residues (Ser1800, Leu1810, Tyr1813, Trp1839, Gly1840) interacting with substrate moiety are evolutionarily conserved ( Additional file
1: Figure S13). On the other hand, only 2 out of the 11 residues (Leu1810, Tyr1813) interacting with the phosphopantetheine group show good conservation, while none of the residues interacting with ACP was conserved. In fact the binding site on the KR domain is located on the surface and many of these interacting residues in KR domain belong to loop regions on the structure. It is possible that the residues on KR domain would have co-evolved with corresponding interacting residues on the cognate ACP domains. The recent studies on the KR domains from hedamycin and actinorhodin clusters from Type II PKSs also showed that the residues interacting with substrate as well as ACP differed within these two domains [
14,
15].
| Table 2Contacts of KR domain with bound substrate and ACP in holo-ACP:KR complex |
The substrate bound holo-ACP:KR complex obtained from docking and energy minimization studies was further refined using molecular dynamics simulations in the explicit solvent environment for a period of 50

ns. Detailed analysis of the conformations obtained from 50

ns simulations indicated that the entire complex showed a backbone RMSD in the range of 3.0 to 5.0

Å from the starting structure. The conformation of the bound substrate and its orientation with respect to the catalytic residues of KR domain was also analyzed to investigate whether the binding pose of the substrate identified by docking and energy minimization was stable during the explicit solvent MD simulation. As can be seen from the lower left panel in Figure

, the distance between keto group at beta position of the substrate and hydroxyl of Tyr 1813 as well as the side chain oxygen of Ser 1800 remains within 3.0 to 6.0

Å throughout the 50

ns trajectory. It has been proposed that during catalysis by KR domain, Lys 1776 is involved in activation of Tyr 1813 to a general acid so that it can donate its proton to the carbonyl oxygen at beta position after hydride transfer from NADPH [
5]. Analysis of the orientations of Lys 1776 over the 50

ns trajectory indicated that it remained in proximity of Tyr 1813 throughout the simulation. The lower right panel in Figure

shows the variation of the distance between NADPH and beta carbon of the substrate over the MD trajectory. As can be seen this distance increases from 3.0

Å in the starting structure to as high as 7.5

Å within 5

ns, but towards the last half of the simulation it comes to a value of 5.0-6.0

Å. These results suggest that, the substrate binding site identified by the current docking and MD simulation studies is consistent with the proposed mechanism for catalysis by KR domain. Figure

shows stereo view of the final complex obtained after 50

ns MD simulation and the contacting residues are depicted in orange color. Figure

a shows the persistence of contacts between the KR domain and the P-pant as well as substrate moiety during the 50

ns simulation. As can be seen, only 17 contacting residues showed interaction with P-pant or substrate for more than 80% of the simulation time. Out of the 21 residues showing contact for less than 60% of the time, 14 were involved in contacts with the substrate while only 7 showed contact with P-pant moiety. Figure

b shows the snapshots of conformations taken up by P-pant and substrate moiety extracted from the 50

ns trajectory at an interval of 1

ns.
The holo-ACP:DH as well as holo-ACP:KR complexes obtained by our study gave the position of ACP domain with respect to DH and KR domains but in a modular PKS the ACP would bind to a complete module structure comprising other domains also. Thus, it was necessary to check if the predicted site of ACP binding on each of these domains is also available in the three dimensional structure of the complete PKS module or occluded by other catalytic domains. Since, no complete module structure is available for PKS, the structure for complete module 4 of erythromycin PKS was modeled using FAS as template using SBSPKS web-server [
29]. The apo-ACP:KR and apo-ACP:DH complexes obtained from protein-protein docking studies were transformed on to this structural model for complete module 4 of erythromycin PKS. As can be seen from Figure

a and b
, even in the complete module structure, ACP can bind to the same sites as predicted by docking and there are no steric clashes of ACP with any other domain in the module. It is interesting to note that the ACP binding site on KR domain is located in a void surrounded by KR, DH, AT and KS domains.