Cation–π and hydrophobic interaction controlled PET recognition in double mutated cutinase – identification of a novel binding subsite for better catalytic activity

Accelerated hydrolysis of polyethylene terephthalate (PET) by enzymatic surface modification of various hydrolases, which would not degrade the building blocks of PET in order to retain the quality of recycled PET, is a promising research area. Many studies have been reported to identify mutations of different hydrolases that can improve PET degradation. Recently, the mutation of glycine and phenyl alanine with alanine in cutinase was found to improve the activity of PET degradation 6-fold. Yet, a deep insight into the overall structural basis as well as the explicit role played by the amino acid residues for PET degradation is still elusive, which is nevertheless important for comparative analyses, structure–function relations and rational optimization of the degradation process. Our molecular dynamics simulations coupled with quantum mechanical study demonstrate that mutations of anchor residue phenyl alanine to alanine at the PET binding cleft of cutinase unveiled a distal yet novel binding subsite, which alters the nature of dispersive interaction for PET recognition and binding. The phenyl alanine engages in π–π interaction with the phenyl ring of PET (−8.5 kcal mol−1), which on one side helps in PET recognition, but on the other side restricts PET to attain fully extended conformations over the entire binding cleft. The loss of π–π interaction due to mutation of phenyl alanine to alanine is not only compensated by the favourable cation–π and hydrophobic interactions from the arginine residues (−17.1 kcal mol−1) found in the newly discovered subsite, but also favours the fully extended PET conformation. This subsequently impacts the overall increased catalytic activity of mutated cutinase.


Introduction
Polymeric esters of terephthalic acid or PET, which is a polyester comprising terephthalic acid (TPA) and ethylene glycol (EG) linked via ester bonds, constitute a good portion of the plastics industry due to their massive production at low cost. [1][2][3][4] The chemical structure of these plastics makes them durable and in turn also resistant to natural processes of degradation, 5 meaning harsh chemical conditions as well as high temperature are required to degrade them. 6,7 Hence, they are a potential threat to human health. 8 However, it is nearly impossible to replace plastics with other materials throughout the various industries. Therefore, it is essential to adopt the concept of reduce, recycle and reuse for a sustainable future, but not much recycling and reuse happens. 9 PET has hydrolysable ester bonds, which can in-principle be cleaved by esterase enzymes, and this makes PET an achievable target for biological recycling. PET can be bio-degraded to monomeric products such as bis(hydroxyethyl)terephthalic acid (BHET), mono(hydroxyethyl) terephthalic acid (MHET), terephthalic acid (TPA) and ethylene glycol (EG), which can be reused for the reproduction of recycled PET of high quality. 10 Thus, enzymatic degradation of PET offers a green and promising alternative with regeneration of monomers, thereby adding to the economy.
The century long exposure of PET to the ecosystems has provided a long time for the bacteria in the nature to evolve and utilize PET as a novel carbon source. 10 In 2016 for the rst time, Yoshida et al. reported an enzyme, IsPETase from Ideonella sakaiensis, which can cleave PET at 30 C. 10 IsPETase, belonging to hydrolase superfamily, shows high substrate specicity to PET, though it suffers from low reactivity and low thermal stability. 11,12 Since its discovery, a lot of research has been focusing on improving the efficiency and thermal stability by introducing mutations [13][14][15][16][17] and other strategies. 18,19 There are approximately 69 IsPETase like enzymes from phylogenetically different organisms, which can in principle degrade the ester linkage present in PET. 14,[20][21][22] Studies to improve the performance of these enzymes is indeed a valuable research area. [23][24][25] In the year 2022, a novel enzyme PHL7 was discovered by Zimmermann and co-workers, which shows promising activity for large scale degradation of PET. 26,27 In the same year, application of machine learning techniques are also used to design FAST-IsPETase, 28 which shows higher activity. 28 Thus, unsurprisingly, the hydrolase family of enzymes from evolutionary diverse sources is greatly explored for the reactivity towards PET degradation. Among them, the bacterial thermophilic cutinase family have emerged as one of the suitable candidate for further research owing to its structural similarity and stability at high temperature between 60-84.7 C, 29,30 which is much higher than the melting temperature (T m ) of IsPETase (46.8 C) 14 and above the glass transition temperature (T g ¼ 65-70 C) of PET. 31 The earlier reports pointed out, that the most reactive bacterial cutinases towards PET are leaf compost 32 (LCcutinase or LCC) and TfCut2 from Thermobida fusca. 33 The reactivity calculated by weight loss of low crystalline PET at 65-70 C with LCC and G62A mutated TfCut2 (ref. 15) aer 24 h are 24% and 20% respectively. 11 Structural similarity of TfCut2 (51%) with IsPETase is slightly more than that of LCC (49%), which makes TfCut2 a valid target for our study. 23 Moreover, the amino acid residues involved in the cleavage of ester bond in TfCut2 are the same as that of IsPETase. 12 Muller et al. in 2005 rst reported the PET degradation by TfCut2 with a degradation rate of 2.5 nmol À1 min À1 cm À2 . 33 Later in 2016, Barth et al. reported TfCut2 could be inhibited by degradation products MHET and BHET, with MHET being the major inhibitor. 15 Though, degradation action of IsPETase is higher than that of WT TfCut2 at their optimum temperatures on PET lm, 10 a notable increase in activity was seen with the G62A mutation on TfCut2. 15 The mutant G62A decreased the binding constant by $5.5 fold leading to relieved product inhibition by MHET than in WT TfCut2. In 2019, Furukawa and co-workers added another mutation F209A to the single mutant resulting in G62A/F209A double mutant (DM) TfCut2 with a PET degradation rate of 15.0 nmol À1 min À1 cm À2 . 11 Interestingly, in 2018 Austin et al. reported a double mutant W130H/S209F IsPETase, which mimics the active site of TfCut2 cutinase showing a higher PET degradation rate. 13 However, the enzymatic PET degradation process is yet to be used in a large scale. 34,35 Therefore, to improve the bio-degradability of PET, it is very important to identify the important amino acid residues present in the enzyme responsible for the interaction with PET and at the same time classify the nature of interaction these amino acid residues exhibits. This process would help engineering of PET degrading enzymes with higher reactivity.
In general, it is observed that mutations can increase the activity of all these enzymes. However, a systematic approach in predicting appropriate mutations to achieve optimum PET degradation is elusive. Thus, our study aims to correlate the increased PET degradation activity of DM TfCut2 with the structure-based role of the mutation, which not only is instrumental for understanding the catalytic mechanism but it will also offer avenues for the applications of genetic engineering strategies to improve the efficiency of enzymes. In doing this, the major challenge is the complexity in characterizing all plausible approaching pathways that can lead to different PET binding conformations on the target protein. In order to address this, the detailed structural information as well as different types of interactions on different PET binding states at the molecular level was analysed based on Molecular Dynamics Simulations and Quantum Mechanical calculations to gain a deep understanding and rational control of PET binding. We have chosen the TfCut2 (WT) and its double mutated (DM) variant G62A/ F209A and a tetrameric PET, 2-hydroxyethyl-(monohydroxyethyl terephthalic acid) 4 (2-HE(MHET) 4 ) for our theoretical study.

Structure preparation of proteins and PET
The PDB structure 4CG1 (resolution 1.40 A) was used as starting structure for WT TfCut2. 36 The missing residue, Asp246 was added using SWISS-PROT. 37 Hydrogens were added using AMBER20 38,39 leap program corresponding to the protonation states predicted by H++ 40 at pH 7. The protein structure thus obtained was minimized for 10 000 steps using the steepest descend-algorithm 41 using the CUDA version of pmemd program in AMBER20. 38 G62A/F209A TfCut2 (ref. 11) structure was prepared from the above WT TfCut2 structure by introducing mutations via PyMOL. 42 2-Hydroxyethyl-(monohydroxyethyl terephthalic acid) 4 or 2-HE(MHET) 4 was taken as a representative of PET polymer. 14 Note that, earlier reports indicate that the tetrameric PET is long enough to establish contact with the entire binding cle of TfCut2. 14 The tetrameric PET structure is minimized using the B3LYP 43 functional and 6-31+g(d) 44,45 basis set using Gaussian 09 program package. 46 The molecular electrostatic surface potential (MESP) of the protein was calculated using PDB2PQR 47,48 online server followed by continuum electrostatics calculations employed in the Poisson-Boltzmann solver (ABPS). 49

Molecular docking
The minimized tetrameric PET and protein structure was used for performing 11 sets of docking on each protein viz. WT TfCut2 and G62A/F209A TfCut2 with PET tetramer using Autodock Vina. 50 The 99 PET conformers thus generated were clas-sied using a distance cutoff of 5 A between any ester carbonyl carbon-atom of PET to the hydroxy-oxygen atom of the catalytic serine. These structures were further grouped using a PET RMSD cut-off of 2 A. This procedure resulted in 19 and 16 different poses of PET with WT and DM TfCut2 respectively, which were subsequently used for Molecular Dynamics Simulations to investigate the conformational landscape.

Parametrization of PET and molecular dynamics simulations
PET was parametrized using the antechamber 51 module of AMBER20 38 using AM1-BCC 52 charge model and the Generalized AMBER Force Field (GAFF). 53 The protein part was treated with the ff14SB 54 forceeld. Proteins with docked tetrameric PET was placed in 0.15 M NaCl with TIP3P 55,56 water in a truncated octahedron box. Non bonded interactions were restricted to 12 A cut-off and long range electrostatic interactions were treated using the particle mesh Ewald method. 57 Solvated protein-PET system was relaxed by minimization of 10 000 steps using the steepest descend-algorithm. 58 The system was gradually heated in 100 K increments from 0 K to 300 K with positional restraints on the protein and tetrameric PET system for 120 ps. Positional restraint was then slowly reduced from 12 kcal mol À1 A À2 to zero for the protein, but by retaining a 0.5 kcal mol À1 A À2 restraint on PET, which was then released followed by a 200 ns production run on each of the structures obtained from docking study. The 4 trajectories where PET moved out of the binding site were not considered further, while the rest was run for 1 ms simulation. This resulted in 15 and 12 simulations each of one ms duration for WT and DM TfCut2, respectively. Also, a total of 1 ms MD simulation on WT and DM proteins without PET were performed.
All MD simulations were performed at a constant pressure of 1 bar in the NPT 59 ensemble. Hydrogen-mass repartition 60 along with SHAKE 61 algorithm has been applied on bonds involving hydrogen atoms. VMD, 62 CPPTRAJ 63 and the Qt version of xmgrace 64 were used for trajectory analysis. The simulations were clustered into different groups by the average-linkage hierarchical clustering algorithm 65 in AMBER20. Clustering using this algorithm is reported to be a reliable approach for analysing the MD trajectories. 65 The trajectories obtained aer clustering were subjected to binding energy calculations using the Molecular Mechanics-Poisson Boltzmann Surface Area 66 (MM-PBSA) method. Though this binding energy calculation does not include the entropic contribution to ligand-binding affinity, it is still reported to give a good correlation with the experimental values. 67 Normal-mode analysis was carried out on extended and folded conformers of PET with WT and DM protein systems using the nabnmode module of AMBER to estimate the conformational entropy contributions to the binding free energy using Kongsted method (Fig. S11 †). 68 This follows the same trend as the MM-PBSA energies, indicating an steady entropy contribution for all the conformers. Further to highlight the most important residues involved in the ligand recognition, MM-PBSA per residue energy decomposition analysis was also carried out.

Quantitative estimation of PET binding energy
The lower order symmetry adapted-perturbation theory calculations, SAPT0 69,70 which comes with a good error cancellation with jun-cc-pvDZ 71 basis was employed with density tting techniques and core electron freezing 70 using the quantum chemistry package PSI4. 72 The total interaction energy in SAPT is a sum of energy terms arising from electrostatic, exchange, induction and dispersion effects. 70 Further, PET binding energy was calculated including dispersion 73 corrections at the B3LYP-D3/6-31+g(d) level of theory by correcting basis set super position error using Gaussian 09 program package. 46

Results and discussion
The amino acid residues involved directly in the ester bond cleavage of PET are conserved in wild type (WT) and double mutated (DM) cutinase, yet the catalytic activity differs signicantly. Hence, the catalytic activity of these enzyme is hugely dependent upon the molecular recognition of the substrate to the PET binding site through multiple weak non-covalent interactions. Concomitantly, the extent of these weak interactions depends upon the orientation of the amino-acid side chain residues and the binding poses of the substrate. A deeper understanding of the complementary molecular reactivity and molecular recognition of PET by cutinase can be obtained by the detailed analysis of the plausible PET binding poses on the protein surface.

PET binding poses on WT and DM TfCut2
In order to gain a preliminary understanding, we have extensively explored 99 possible binding poses of tetrameric PET the with the PET binding cle of each of WT and DM TfCut2 by molecular docking using AutoDock Vina. 42 The DM TfCut2 has an average binding affinity of À6.5 kcal mol À1 (range: À5.0 to À7.5 kcal mol À1 ) which is a little higher than the average binding affinity of WT TfCut2, À5.5 kcal mol À1 (range: À4.2 to À7.2 kcal mol À1 ). 46 conformations show a docking score higher than the average binding affinity in WT TfCut2, while 72 conformations of tetrameric PET in DM TfCut2 have a docking score higher than the average. However, all the binding poses with higher docking score are not identical, which warrants more thorough analysis of the different binding poses to enable structural basis for enhanced degradation due to mutation.
The observed PET binding poses from docking can be broadly classied into three categories, (i) extended, where the entire tetrameric PET maintains contact with the entire binding cle (ii) semi-extended where the tetrameric PET does not bind throughout the entire region of binding cle and (iii) folded, where the tetrameric PET forms a loop like structure, binding only to a very small region (Fig. 1). Austin et al. 13 and Joo et al. 14 in two separate works had reported that tetrameric PET adopting an extended conformation in IsPETase shows higher docking score as well as higher reactivity. Likewise, we have identied the extended PET conformers among the docked structures for both WT and DM TfCut2. Only 15 docked poses exhibit extended conformation in WT TfCut2 (Fig. 2a), while as much as 51 docked poses exhibit extended conformation in DM TfCut2 (Fig. 2b). Thus, the binding site cle of DM TfCut2 seems more suitable to extended PET conformations, whereas folded PET conformations are more common for WT TfCut2. A comprehensive analysis of the docked conformations of PET indicates that PET binds to three-fourths of the active site (up to F209) in a similar fashion in both WT and DM TfCut2. The visual inspection revealed that DM TfCut2 uses an unprecedented terminal region of the protein for binding PET in the extended conformation as compared to WT TfCut2 (Fig. 2). Earlier binding studies did not indicate PET binding to such residues. 13,14,36,[74][75][76][77] At this point we might reason that, the mutation of the aromatic Phe209 (F209) in WT to hydrophobic Ala209 (A209) in DM cutinase might have affected the conformation of PET binding on the surface of cutinase. However, further analysis is required to elucidate this striking difference in the mode of PET binding upon mutation.
From theoretical mechanistic studies on the enzymatic PET hydrolysis, [78][79][80][81][82] it was proposed that the hydroxyl group of Ser130 of the catalytic triad initiates the attack to the carbonyl carbon of the ester bond, followed by the stabilization of the negative charge by the two back bone nitrogen atoms of Tyr60 and Met131 of the oxyanion hole. 83 Thus, in order to identify the best binding poses of PET with WT and DM TfCut2 in terms of catalytic activity, it is very important to analyse the proximity of the carbonyl carbon-atom of PET to Ser130 that lie in the subsite S1. Thus, a distance cut off of 5 A is used to identify 62 and 65 PET poses respectively for WT and DM TfCut2. These structures were further grouped together by a PET RMSD cut-off of 2 A, and one representative pose with the highest docking score was selected from each group. This resulted in shortlisting of 19 and 16 PET complex of WT TfCut2 and DM TfCut2 respectively ( Fig. S2 and S3 †), which were further subjected to detailed analysis to identify what might be responsible for PET binding to novel sub-site in DM cutinase.

Analysis of structure and electrostatic potential of WT & DM TfCut2
A closer look at the entire PET binding cle of cutinase reveals that it constitutes of 36 amino acid residues (Fig. 3). These amino acid residues can be divided into two categories, the ones that are directly involved in the catalytic cleavage of the ester bonds or catalytic residues, the others that are involved mainly in the recognition and binding of PET on the protein surface or anchor residues. One can classify the entire cle into four major subsites viz. S1, S2, S3 and S4, where each subsite can roughly accommodate a monomeric PET unit. A similar description of the entire binding cle was also provided for PET degrading enzyme IsPETase. 14,77 However, our docking studies have iden-tied an unprecedented binding subsite, named as S4 0 that seems an alternative to S4. Note that, the binding of PET into S4 0 was not previously reported. Also, the presence of Arg18 at the extreme end of S4 0 have close resemblance with the binding sub-site S4 of IsPETase, where Arg251 marks the end of the binding cle. 14 For easy identication purpose, the amino acid residues in each subsite are colour coded as green (S1), yellow (S2), red (S3), blue (S4) and cyan (S4 0 ) respectively for both WT and DM TfCut2 (Fig. 3a and b). Fig. 3c and d indicates that, the conventional binding cle constituted by subsites S1, S2, S3 and S4 has a narrow and inverted V-shape (Q92, W155, F209/A209, E251, K216, S66), while the binding cle constituted by subsites S1, S2, S3 and S4 0 appears to be wider and linear (Q92, W155, F209/A209, K216, S19, S66). Also, the distance between the C-alpha atoms of the terminal residues shows that the mutation does not alter the overall size and shape of the cavity. The plot of the molecular electrostatic potential (MESP) on the van der Waals surface clearly indicates that the S1, S2 and S3 subsites are mostly neutral with slight negative region in S1 and S3 due to the presence of aromatic Trp155 (W155), Tyr60 (Y60) and Phe209 (F209) residues while the S4 subsite is mainly polar with both electrophilic and nucleophilic regions near Lys216 (K216) and Glu251 (E251) respectively (Fig. 3c). Similar to S4, the alternate subsite S4 0 also contains both nucleophilic (Trp69 and Glu72) and electrophilic (Arg73 and Arg18) residues, indicating similarity in the molecular electrostatic potential. The mutations of Gly62 by Ala62 in S2 and Phe209 by Ala209 in S3 did not change the overall MESP. However, local decrease in the electron density near the Phe209 in WT TfCut2 is observed aer the mutation (Fig. 3d).
Molecular-dynamics simulations on WT and DM (G62A/ F209A) TfCut2 with and without PET 1 ms Molecular Dynamics simulations was carried out on WT and DM TfCut2 with and without PET to identify the amino acid residues and their corresponding nature of interactions responsible for PET recognition. The average RMSD value of the WT and DM TfCut2 without the tetrameric PET is 2.60 A and 1.93 A respectively (Fig. 5a), depicting overall stability of the proteins. The RMSD of the average structures generated aer 1 ms Molecular Dynamics simulations on WT and DM TfCut2 without PET, showed an average RMSD of 1.33 A of the DM TfCut2 with respect to the WT variant. The RMSF analyses indicates higher uctuation of the N-terminal and C-terminal residues ( Fig. 4b and S1 †). The low RMSD of the average structures (0.83 A) by excluding the N-terminal and C-terminal Fig. 3 Cartoon representation of (a) WT TfCut2 and (b) DM TfCut2 with amino acid residues in subsites S1, S2, S3, S4 and S4 0 shown in green, yellow, red, blue and cyan licorice representations respectively. Plot of molecular electrostatic potential (MESP) for (c) WT TfCut2 and (d) DM TfCut2 with potential over the range À8.0 to +8.0kT/e. residues highlights the importance of the local conformational changes for PET recognition.
To explore the binding interaction of PET with the side-chain in WT and DM TfCut2, 1 ms MD simulation was run on each of the selected 15 and 12 PET conformers for the WT and DM TfCut2 respectively. The analysis of the RMSD (Fig. S4 and S5 †) shows that though the RMSD of the protein-PET complexes is signicantly higher, the RMSD of protein alone in the protein-PET complex is similar to that of the proteins. Thus, the signicant change in the RMSD of the protein-PET complex is mainly contributed from the different conformations adopted by PET in the PET binding cle. Looking at the trajectories, conversion between the folded and extended conformers of PET is widely observed, validating the possible binding of both the conformers in the binding cle. In addition, the RMSF plot of the cavity residues with and without PET indicates that the subsites S1 and S3 undergo greater extent of conformational changes in WT TfCut2, while subsite S1, S4/S4 0 undergo major changes in DM TfCut2 for PET binding (Fig. S6 †).
We have further used hierarchical average-linkage clustering implemented in AMBER20 to identify the major conformational isomers of PET from the MD simulation trajectories of the PET bound WT and DM TfCut2. Clustering of all the trajectories for WT and DM TfCut2 produced 78 and 67 clusters respectively ( Fig. 5 and 6). The clustering analysis of WT TfCut2 clearly indicates that, majority of the trajectories yield a single preferred conformation, which is corroborated well with the minimal uctuation observed in the RMSD and RMSF (Fig. S4-S6 †). The WT TfCut2 shows a clear preference for PET conformers that are interacting with the Phe209 of the S3 subsite leading to either a semi-extended or a folded conformation. There are only a few instances where PET attain a fully extended conformation by binding throughout all the four subsites S1-S4/S1-S4 0 . The binding affinity of all the PET conformers in each cluster was calculated using MM-PBSA analysis, which ranges from À15.2 kcal mol À1 to À38.1 kcal mol À1 (Fig. 5). Among the prominent clusters, the higher binding energy is mostly observed for the semi-extended PET conformers, which indicates the prominent role of Phe209.
On the contrary, the clustering analysis on the trajectories of DM TfCut2 reveals conformational exibility of PET in the binding cavity leading to more than one prominent clusters. This corroborates with the higher RMSD and RMSF as compared to WT cutinase (Fig. 6). Note that, majority of them indicates binding to the hitherto unidentied subsite S4 0 leading to fully-extended conformations more oen than WT TfCut2. Also, the binding energy of PET in DM TfCut2 ranges from À10.2 to À48.8 kcal mol À1 and the binding energies for the fully extended PET conformers binding to S1-S4 0 subsites are noticeably higher. Therefore, the higher PET affinity of DM TfCut2 than WT variant can possibly result from binding to the novel subsite S4 0 , which also leads to extended PET conformations. Note that, the extended PET conformations in the binding cle has been linked with higher activity of IsPE-Tase. 13,14 Thus, extended PET conformations with relatively higher binding energy might lead to higher activity in TfCut2 as well.
In order to identify the residues that might play important role in energetically favouring extended PET conformations over the folded ones, we have selected the trajectories corresponding to extended and folded conformations having the ve highest and ve lowest binding energies in both WT and DM TfCut2. The average structures of these 20 conformers are given in Fig. S7 and S8. † Note that, fully-extended PET conformations are only observed for DM TfCut2, while the WT TfCut2 exhibits semi-extended conformations. The Fig. S9, S10 and S12 † show the contact plot analysis, MM-PBSA per-residue binding energy and comparison of average per residue binding energy with the average close contact respectively calculated using MM-PBSA method for each subsite on the extended and folded PET for WT and DM TfCut2. The pair-wise binding energy for PET with all the residues of the entire PET binding cle of WT and DM TfCut2 (S1, S2, S3, S4 and S4 0 ) was plotted in Fig. 7. We can clearly see different interaction energy pattern for the WT and DM TfCut2 as well as for the extended and folded PET. In general, the residues with highest number of close contacts have higher interaction energy with PET. We have also performed zero order Symmetry Adapted Perturbation Theory (SAPT0) analysis to get a quantitative picture of different energy components contributing to the total energy.
The S1 residues Ser130 from the catalytic triad and Tyr60, Met131 from oxyanion hole that are essential for the ester bond  cleavage of PET as well as the anchor residues such as hydrophobic Ile178 and aromatic Trp155 that are essential for PET binding exhibit more close contact and higher energy contribution irrespective of the conformational change and mutation (Fig. S9, S10 and Movie S1 †). Thus, the residues in S1 can be considered as the most important for the PET recognition and cleavage by cutinases. Hence, any mutation at S1 might compromise the recognition and degradation of PET. Indeed studies on IsPETase, which also contain the same residues, shows decline in activity upon mutation of S1 residues. 14,84 Tyr60 of the oxyanion hole is a non-polar aromatic residue with a PhOH functional group. Three main interactions viz., displaced p-stacking ( PET Ph/Ph Y60 ), T-stacking ( PET C-H/ Ph Y60 and Y60 C-H/Ph PET ) are mainly observed between Tyr60 and PET along with hydrogen bonding between amide hydrogen of Tyr60 and carbonyl oxygen of PET ( PET O/H-N Y60 ) (À10.7, À7.8 and À7.0 kcal mol À1 in Fig. 8a-c and Movie S1 †).
Note that, the hydroxyl group of Tyr60 does not interact with non-polar functional groups in PET. The other two residues, Met131 and Ser130, responsible for the ester bond cleavage show lesser number of interactions and less energy contribution as compared to Tyr60. The non-polar Met131 is involved in weak M131 C-H/p PET , hydrophobic M131 C-H/H-C PET and PET O/H-N M131 hydrogen bonding (À4.3 kcal mol À1 , Fig. 8d) and the polar Ser130 show occasional hydrogen bonding interaction S130 O-H/O]C PET with the carbonyl group of the PET (À6.2 kcal mol À1 , Fig. 8e) in the neutral state, which resembles the state prior to proton abstraction by His208 in the catalytic cycle. 78,80 The aromatic residue Trp155, which is responsible for anchoring the phenyl ring of PET to S1, shows PET C-H/p w155 (À8.1 kcal mol À1 , Fig. 8f) interaction and weak T-stacking interaction with PET (À5.0 kcal mol À1 , Fig. 8g). The hydrophobic residue Ile178 shows mainly weak I178 C-H/p PET interaction with one or both phenyl rings of PET Fig. 6 The output of clustering on 12 different DM TfCut2-PET systems with the population percent of each cluster shown in Y axis. The average structure of each cluster along with its binding energy (kcal mol À1 ) obtained via MM-PBSA calculations is shown above the bar corresponding to each cluster. A zoomed image of PET is shown on top of each average structure.
(À6.5 kcal mol À1 , Fig. 8h). The interaction energy of Tyr60, Ile178 and Trp155 are found to be highest in the folded PET conformer with double mutated TfCut2, explaining the overall higher interaction energy for the folded PET conformer in DM TfCut2 as compared to WT TfCut2.
Other than Ser130, the major contribution to the total interaction energy with PET comes from dispersion component indicating the importance of hydrophobic interaction (Fig. S13 and S14 †). However, the polar Ser130 have electrostatic component as major contributor to the total interaction energy ( Fig. S13 and S14 †). Note that, even though the energy corresponding to the individual interaction energy is very less, a signicant contribution arises from the cooperative stabilizing effect of multiple non-covalent interactions leading to preferential PET recognition.
Weak hydrophobic C-H/H-C interactions of PET are mainly observed with Gly59 and Gly62 (À4.9 and À1.1 kcal mol À1 , Fig. 9a and b) from S2. The mutation of Gly62 with Ala62 (À1.5 kcal mol À1 , Fig. 9c) does not affect the interaction energy, hence might not play important role in PET binding and cleavage. Indeed, the major role of G62A mutation was reported to assist in release of product. 74 Polar hydroxyl group of Thr61 shows hydrogen bonding PET Fig. 9e). Note that, mutation of Trp159 in WT IsPETase with His159 increases the activity signicantly. 13 Thus, S2 residues in general show less contact and hence lower interaction energy with PET as compared to S1, irrespective of the conformational change of PET as well as G62A mutation of TfCut2. Also, the presence of non-polar residues in S2 is corroborated with the dispersion interaction being the major energy component (Fig. S13 and S14 †).
The next subsite that shows signicant close contact with PET is S3, where the aromatic residue Phe209 in WT TfCut2  exhibit greater contact with both the semi-extended and folded PET conformers through displaced p and T-stacking interactions with the phenyl ring of PET (À8.5 kcal mol À1 , Fig. 10a and Movie S2 †). When F209 is mutated to aliphatic Ala209, the p-p interactions are no longer present rather weak hydrophobic C-H/H-C and hydrogen bonding PET C]O/H-C A209 interactions are observed, which is also reected in the signicant lowering of the interaction energy upon mutation (À1.9 kcal mol À1 , Fig. 10b and Movie S3 †). In addition, semi-extended PET in WT TfCut2 exhibit major interaction with polar Asn212 and hydrophobic Ile213 from S3, which are not found in folded PET conformers (Movie S4 †). Polar Asn212 shows hydrogen bonding  Fig. 10c), while nonpolar Ile213 is involved in I213 C-H/p PET and I213 C-H/H-C PET interactions (À6.8 kcal mol À1 , Fig. 10d). However, mutation of aromatic residue Phe209 to hydrophobic Ala209 had reduced the close contact of both extended and folded PET with S3 of DM TfCut2, while maintaining the close contact with Asn212 and Ile213 in the extended form.
In addition, Pro214, which is situated at the end of S3 shows signicantly high close contact and interaction energy contribution due to P214 C-H/p PET interaction, indicating the preference of the extended PET in DM TfCut2 (À8.2 kcal mol À1 , Fig. 10e) as compared to semi-extended PET in WT TfCut2. Ser66 on the other hand, though a polar residue, mainly shows weak hydrophobic interactions S66 C-H/H-C PET in both WT and DM cutinase (À4.2 kcal mol À1 , Fig. 10f). Likewise, the contribution from the dispersion energy is major in all residues other than polar Asn212, where electrostatic component is higher ( Fig. S13 and S14 †).
Note that, no residues from S4 in both WT and DM TfCut2 show any close contact with either semi-extended, extended or folded PET conformers. Accordingly, the energy contribution of S4 residues is minimum (Fig. S12 †). The striking difference in terms of close contact with the extended PET in DM TfCut2 as compared to semi-extended PET in WT Tfcut2 is observed for the novel S4 0 subsite residues. The MESP analysis ( Fig. 3c and d) revealed that the subsite S4 0 contains both polar and non-polar residues. The energy contribution from nonpolar aliphatic residues Ile67, Ala68, Leu70 and Gly71, as well as polar Ser19 is not signicant. However, the aromatic residue Trp69 shows Tstacking and W69 C-H/p PET interactions with the phenyl ring of PET (À5.6 kcal mol À1 , Fig. 11a). The carboxylic group of Glu72 interacts with PET via E72 C]O/H-C PET and hydrophobic C-H/H-C interactions (À2.0 kcal mol À1 , Fig. 11b). Arg73 shows strong hydrogen bonding interaction with the hydrogen of guanidino group and carboxyl oxygen of PET ( PET C]O/H-N R73 ) and PET C-H/H-N R73 interactions À17.1 kcal mol À1 , Fig. 11c). The terminal Arg18 having major contribution shows cation-p interaction with the phenyl ring of PET in addition to the PET C]O/H-C hydrogen bonding (À9.9 kcal mol À1 , Fig. 11d and Movie S5 †). The SAPT0 breakdown of the total interaction energy indicates that the major contribution comes from the dispersion components except for polar Arg residues, which have high electrostatic component. The signicantly higher close contact and interaction energy with polar S4 0 residues like Arg18, Trp69, Arg73 and Glu72 indicates that, the S4 0 plays a vital role in attaining favourable fully extended PET conformation in DM TfCut2, which might be prevented by the p-p interaction of the aromatic Phe209 in S3 of WT TfCut2.
From these results, one may infer that S1 is important for both binding and cleavage of PET and mutation in S1 might cause reduction in the PET degradation activity either due to poor binding or poor cleavage. In accordance with that, S1 shows very high binding energy irrespective of the conformational change as well as mutation (Fig. S14 †). The next important subsite is S3, which plays very important role in anchoring the aromatic phenyl ring of PET with Phe209 by p-p interaction in WT TfCut2. However, note that, this interaction also restricts PET to attain a favourable extended conformation as well as might delay the sequential approach of the ester bonds close to the catalytic residues situated in S1 (Scheme 1). The loss of p-p interaction in DM cutinase due to F209A mutation is compensated by cation-p, hydrogen bonding and hydrophobic interactions from the Arg18, Trp69, Arg73 and Glu72 residues in S4 0 , which not only gives higher binding energy but also assist the Fig. 11 Snapshots showing major interactions in subsite S4 0 (color code same as above). Cation-p interactions which are unique in S4 0 are shown in orange color. The single letter code indicating the amino acids are Trp (W), Glu (E) and Arg (R).
Scheme 1 Schematic representation of plausible PET degradation pathway in WT and DM cutinase. The subsites are color coded as green, yellow, red, blue and cyan respectively for S1, S2, S3, S4 and S4 0 . The TPA and EG moieties of the PET polymer are presented with hexagons and lines, respectively. fully extended PET conformation. A combination of these factors may in turn facilitate higher degree of PET degradation by DM TfCut2. Therefore, an optimum balance between the two opposing effect viz. the intrinsic PET affinity to the enzyme surface and PET binding to the active site seems to control the extent of degradation. While our manuscript was under review, Kari and Westh have showed that Sabatier principle might be helpful in rationalizing the optimum enzymatic PET degradation processes, 25 which supports our analysis.

Conclusions
In this study, we have explored the different conformational states and binding affinity of the tetrameric PET with the amino acid residues from the PET binding cle of wild type and double mutated cutinase using molecular dynamics simulations and quantum mechanical approaches. We demonstrate that mutations of anchor residue phenyl alanine to alanine at the PET binding cle of cutinase alters the conformation of PET on the binding surface of the mutated cutinase demonstrating structural adaptability of ligand on the protein surface. This unveiled a distal yet novel binding subsite, which alters the nature of dispersive interaction for PET recognition and binding. The phenyl alanine engages in p-p interaction with the phenyl ring of PET (À8.5 kcal mol À1 ), which at one side helps in PET recognition, but at the other side restricts PET to attain fully extended conformations over the entire binding cle. The loss of p-p interaction due to mutation of phenyl alanine to alanine is not only compensated by the favourable cation-p and hydrophobic interactions from the arginine residues (À17.1 kcal mol À1 ) found in the newly discovered subsite, but also favours the fully extended PET conformation. Thus, a perfect balance between the PET recognition and binding by the anchor residues and attainment of suitable conformation for the cleavage would likely be the key for a higher PET degrading activity. Clearly, further investigation is needed to identify the best possible mutations to be engineered for more efficient enzymatic PET degradation.

Author contributions
SD designed the study. AJ performed the simulations and calculations. SD and AJ did the analysis of the data, discussed the results and wrote the manuscript.

Conflicts of interest
There are no conicts of interest to declare.