Lid opening and conformational stability of T1 Lipase is mediated by increasing chain length polar solvents

Preparation of protein and organic solvents structures

Crystal structure coordinates of the closed conformation of wild-type Geobacillus zalihae T1 lipase (PDB ID: 2DSN) (Matsumura et al., 2008) were obtained from the Protein Data Bank with a resolution of 1.5 Å. The starting coordinates, crystal water molecules and metal ions (Ca²⁺ and Zn²⁺) were reserved. The structures of organic aliphatic chain length solvents were taken from the PubChem Compound Database NCBI server (entries CID_6276, CID_263, CID_702, CID_887, CID_962, CID_1031), and are summarized in Table 1. To minimize the number of energetically unfavourable structural features, the protein structure was energy minimized and the hydrogen bonding network was optimized using YASARA force field (Krieger et al., 2004). The organic solvents were modelled using the automatic parameterization algorithm (AutoSMILES), which can derive force field parameters for organic molecules fully automatically (Wang et al., 2004; Jakalian, Jack & Bayly, 2002). This parametrization also considers force field parameters for small organic molecules.

Organic solvent stability assay of native T1 lipase

Organic solvent stability of T1 lipase activity was conducted by measuring the residual lipase activity quantitatively from the free fatty acid released based on the oleic acid standard according to Kwon & Rhee (1986). T1 lipase was purified by affinity chromatography with Ni²⁺ Sepharose loaded onto an XK16 column mounted on an AKTA Prime (GE healthcare). The purified enzyme (0.5 mg/mL) was incubated in five different polar solvents (methanol, ethanol, propanol, butanol, and pentanol) in a volume ratio of enzyme to organic solvent representing 25% and 50% final concentrations. The solution was thoroughly mixed and incubated at 70°C in a water bath shaker at 200 rpm for 1 h. Substrate emulsion (2.5 mL olive oil/Glycine-NaOH buffer pH 9 at 1:1 ratio, 20 μL 20 mM CaCl₂) was then added and incubated at 70°C with continuous shaking at 200 rpm for 30 min. The enzyme reaction in the emulsion system was stopped by adding 6 M HCl (1 mL) and isooctane (5 mL), and vortexed for 30 min. The upper isooctane layer containing the fatty acid was transferred to a test tube for analysis. Copper reagent (1 mL) (5% (w/v) copper (II) acetate-1-hydrate adjusted to pH 6.1 with pyridine) was added and vortexed for 30 s and allowed to stand for 2 h. The absorbance of the upper layer was read at 715 nm. One unit of lipase activity was defined as the releasing of 1 μmol of free fatty acid per minute.

Preparation of the molecular dynamics setup

The YASARA Structure software package version 12.10.3 (Krieger et al., 2004), installed on a Dell XPS PC (Intel Core TM i7-6700 3.41 GHz), was used to perform all simulations. A rectangular cubic simulation box was put around the protein with a distance of 5.0 Å  around all atoms. Boundary conditions were set to periodic in the x, y and z directions. Due to the size and amount of solvent mixture molecules, the dimensions of the simulation box were determined automatically using the Auto parameter. This was done to avoid the risk of violating the minimum image convention due to interactions of the solvent molecules as they move across the periodic boundary. To make sure that the cell is neutral and protonation states of amino acid side-chains are assigned, a cell neutralization and pKa prediction experiment were performed using the AMBER03 force field (Duan et al., 2003). This places the ions at the locations of the lowest/highest electrostatic potential until the cell is neutralized and the requested ion mass fraction (0.9% NaCl) is reached. The water molecules were removed leaving the counterions. Each simulation box was then filled with the calculated copies of water/organic solvent molecules in order to achieve a concentration of 50% (Table 1). Since the number of solvent molecules required to attain a given activity cannot be established ab-initio (Wedberg, Abildskov & Peters, 2012), preliminary experimental assays guided the selection of the number of solvent molecules added in the system. The number of solvent molecules was adjusted by considering the experimental density at 298 K and 1 bar. For simulations involving pure water as a control, TIP3P water with flexible bonds modelled with AMBER03 force field was used (Miyamoto & Kollman, 1992) and the pressure rescaled until the experimental density of 0.997 g/ml was achieved.

Molecular dynamic simulations protocol

The atoms of the protein structure in the different solvent mixtures in the simulation box were constrained (fixed). This is to equilibrate the solvent mixtures in the box and allow the relaxation of the solvent molecules and prevent large distortions to the protein. To remove bumps, delete unwanted water molecules and correct the covalent geometry, the atoms were unconstrained and the system setup was energy-minimized with the AMBER03 force field (Duan et al., 2003). The geometric parameters (dihedral angles, bond angles, and length) including Van der Waals and Coulomb interactions of the force field were truncated at a 7.86 Å  force cut-off. The Particle Mesh Ewald algorithm (Essman et al., 1995) was used to treat long-range electrostatic interactions with a grid spacing of <1 Å. After removal of conformational stress by a short steep descent minimization of 200 steps, the procedure continued by simulated annealing (time step 2 fs, atom velocities scaled down by 0.9 every 10th step) until convergence was reached. This final system setup was used for the production runs.

All MD simulations were carried out in three independent runs for each system over a 40 ns timescale using the AMBER03 force field on YASARA Structure software package. Different RandomSeed numbers based on the temperature control (343K) of YASARA documentation were assigned using a Maxwell–Boltzmann distribution to achieve different initial velocities for each replicate simulation. The velocities of the atoms were rescaled according to the Berendsen barostat as described by Krieger et al. (2004), and pH 9 was used in all solvents. The manometer pressure mode was used to control the pressure from the kinetic energy/virial and slowly rescale the simulation cell along the X, Y,  and Z-plane to smoothly approach 1 bar. Simulation snapshot coordinates were saved every 30 ps.

Analysis of data and evaluation

Unless otherwise specified, all analyses of trajectories were performed from replicates of three different simulations run at different initial velocities. Molecular graphics were created with YASARA (http://www.yasara.org) and POVRay (http://www.povray.org) (Krieger & Vriend, 2014). The structural deviations of the backbone atoms (RMSD) of the simulation runs from the initial crystallographic structure, the B-factors from the root mean square fluctuations (RMSF), accessible solvent area (SASA) and secondary structure interactions were analysed by executing the YASARA macros protocol described in detail by Krieger & Vriend (2015). The RMSD is expressed by the difference between the Cartesian atom coordinates based on the selection formula ${\displaystyle RMSD=\sqrt{\frac{\sum _{\mathfrak{i}=1}^{\mathfrak{n}}{\Re }_{\mathfrak{i}\ast }{\Re }_{\mathfrak{i}}}{\mathfrak{n}}}}$ where ‘R’ is the vector linking the ‘n’ corresponding atom pairs in space.

The B-factor as a measure of the RMSF was obtained with the following conversion factor ${\displaystyle RMS{F}_{\mathfrak{i}}=\sqrt{\sum _{j=1}^3\left(\frac{1}{N}\sum _{k=1}^N{P}_{ikj}^2-P_{ij}^{-2}\right)}}$ ${\displaystyle Bfacto{r}_i=\frac{8}{3}{\left(\pi \ast RMS{F}_i\right)}^2}$

Where the RMSF of atom i in Å  and j runs over the three Cartesian components x, y and z of the atom position vector P, and K runs over the N coordinate sets. The solvent accessible area (SASA) was calculated using the default YASARA’s ‘Numeric’ algorithm with the radius of the water probe (solvation radii) set at 1.4 Å.

Structural and conformational stability analysis

The contributions of hydrophobic and hydrogen bonds to protein stability was analysed. The hydrophobic interaction was calculated using YASARA’s knowledge-based potential. Hydrophobic atoms were identified and assigned an atom type. Type 1 carbons had ≥ 3 hydrogens (−CH₃), type 2 had two hydrogens or one hydrogen plus three carbons bound (−CH₂ −, HCC3) and type 3 had aromatic rings with only carbons and hydrogens bound. The interaction range was determined based on the minimum distance (where the energy becomes positive due to clashes), optimum distance (the energy minimum) and maximum distance (where energies become positive due to the cost of creating a vacuum between the atoms). Hydrophobicity based on the Eisenberg hydrophobicity scale (Eisenberg, Weiss & Terwilliger, 1984) was mapped on the surface using PyMOL 1.7.2.1 (DeLano, 2002).

The hydrogen bonds were analysed based on YASARA’s definition of hydrogen bond formation, whereby the analysis of bond energy as a function of the H-bond-acceptor distance in Å  and two scaling factors is derived by the formula: ${\displaystyle EnergyHBo=25\ast \frac{2.6-max\left({Dis}_{H-A^{2.1}}\right)}{0.5}\ast {Scale}_{D-H-A}\ast {Scale}_{H-A-X}.}$

The secondary structure analysis was carried out using the Kabsch & Sander (1983) algorithm incorporated in their Dictionary of Secondary Structure for Proteins (DSSP) program to analyse the variation of protein secondary structure changes. The secondary structure assignment algorithm program installed onto YASARA was used to analyse the change in secondary structure contents.

Time-correlated atomic motions

To identify correlated and concerted motions, we calculated the dynamic cross-correlation map (DCCM), which is a matrix representation of the time-correlated information between protein atoms i and j, (c_ij) (Hünenberger, Mark & Van Gunsteren, 1995). Dynamic cross-correlation maps were used to detect time-correlated motions in the protein as implemented by Krieger et al. (2004) in YASARA. Dynamic cross-correlation matrices between units i and j were obtained with the following expression: ${\displaystyle DCC{M}_{i,j}=\frac{〈{\overrightarrow{d}}_i.{\overrightarrow{d}}_j〉}{\sqrt{〈d_i^2〉}〈d_j^2〉}}$

Where d is the displacement between the current position and the average position of all selected pairs of atoms and how their movements correlate. The angle brackets indicate the average over all trajectories. The values in the DCCM range from −1 (perfectly anti-correlated) to +1 (perfectly correlated). The values along the diagonal are always +1 (because the motion is perfectly correlated to itself).

Stability property of MD simulations

T1 lipase (Fig. 1) is considered a model lipase that could meet the demands of industrial biocatalysts. The general essential dynamics in identifying the structural changes that occurred during the molecular dynamics simulations over the timescale of 40 ns describe the important contributions of polar organic solvents to the general stability of T1 lipase. The time-average variations of backbone RMSD from the initial protein structure of T1 lipase simulations in water-organic solvent environments over a timescale of 40 ns are reported in Fig. 2 and Fig. S1. For comparisons, the simulations in water are also reported. The overall RMSDs are found to have a displacement difference across all solvent mixtures as they approach a point of convergence at 40 ns for the most equilibrated parts. Variations in RMSD for the simulations in H₂O maintain a stable state and appear well equilibrated. The higher average RMSD values for simulations in solvent mixtures compared to water in Table 2 indicated structural differences in the protein backbone in these media. Accordingly, the RMSD of the backbone atoms increases with decreasing polarity of solvent mixtures. The conformational differences in RMSD are consistent with the role of solvents in having a different degree of hydration and dehydration effects on enzymes because the lower polarity of longer carbon chain length solvent mixtures can weaken interactions stabilizing the structure (Diaz-Garcia & Valencia-Gonzalez, 1995). The SASA (Fig. 2B and Fig. S2) deviates among all organic solvent mixture simulations. Higher SASA in some solvents can result in increased hydrophobic areas of residues, ultimately reducing the compactness of the hydrophobic core of the protein. Gyration radius (Rgyr) (Fig. 2C and Fig. S3) reveals the compactness of the simulations, and the increase in the Rgyr of all solvent mixtures is consistent with the solvation of charged and polar residues on the protein surface (Arcangeli, Bizzarri & Cannistraro, 2001). The SASA and Rgyr reported here are consistent with previous findings where the stable curves at equilibration keeps the T1 lipase in a suitable configuration that will favour its catalytic activity (Wang, Wei & Wang, 2010).

Flexibility of T1 in solvents

The flexibility of the residues throughout the time-course was inspected from the calculated measure of the B-factors based on the analysis of the root-mean-square fluctuations (RMSF) (Fig. 3, Fig. S4 and Table 2). The B-factor or temperature factor represents a smearing of atomic electron densities as it governs the equilibrium positions of the atoms due to the thermal motion and positional displacement (Parthasarathy & Murthy, 2000). High thermal motions of residues are proposed targets towards improving the stability of proteins. There is considerable experimental evidence of a strategy that targets residues with high B-factors and has been shown to improve the organic solvent stability of a Bacillus subtilis lipase (Reetz et al., 2010). In this regard, the atomic motion of the residues was high in response to decreasing solvent polarity indicating increased flexibility. In particular, the flexibility of residues was more confined to the loop region at the C-terminal in all solvent simulations, with decreasing flexibility and increased rigidity of the N-terminal region in H₂O, MtOH-H₂O, EtOH-H₂O, and PrOH-H₂O. However, two distinguishable significant fluctuations were by residue Arg103 and Leu277. These residues are positioned on the surface of the protein and exposed directly to the solvent environment. Furthermore, we investigated the unambiguous effect of exposure of the amino acid residues to the solvents and other interactions that could influence the flexibility. We found that the amino acid residues deployed in close proximity to Arg103 were non-aromatic side chain residues. The separation distance between the cation group of Arg103 and the aromatic ring residues is far apart, which could not have allowed Arg103 to experience a favourable cation-π interaction, thereby becoming energetically stable and increasing its flexibility (Gallivan & Dougherty, 1999). Similarly, no neighbouring acidic residues for the formation of salt bridges occur to stabilize the loop region of Arg103 to restrict it conformational flexibility. The observed flexibility of Leu277 suggest that it cannot interact strongly with other neighbouring residues that are not hydrophobic. Therefore, the direct exposure of these residues to the solvents triggered their flexibility. Strub et al. (2004) identified mutational substitution to lower the surface area of hydrophobic residues with amino acid arginine, and to significantly enhance the stability of protein in organic solvents. Highly flexible residues could trigger protein unfolding due to their large fluctuation. Moreover, the flexibility can also be used as an indicator to find unstable residues whose cooperative interactions can be reinforced to improve stability (Bae & Phillips, 2005; Joo et al., 2011; Reetz, Carballeira & Vogel, 2006).

To further validate the relationship between flexibility and surface area/hydrophobicity of residues, the calculated RMSF and SASA from the last 150 ps trajectories were also obtained (Fig. 4). This is important because the RMSF/B-factor of the averaged equilibrium conformation often served as an indicator to estimate flexibility of the protein (a property said to be related to the protein activity) (Li et al., 2010; Yang, Dordick & Garde, 2004). However, to achieve flexibility, the residues are often exposed to the surrounding media thereby increasing their solvent accessibility and contact. In Fig. 4, the higher flexibility indicates that both the hydrophobic and hydrophilic side chains reorient themselves on the protein surface in organic solvents (particularly for Arg103 and Leu277). Similar observations have been reported for other organic solvent media with variations of the solvent composition (Zhu et al., 2012). Usually, solvent-buried residues would have the lower RMSF/B-factor and solvent-exposed residues would have the higher RMSF/B-factor. The correlations in Fig. 4, derived from other reports, describe the correlation between the flexibility of protein and its total SASA, wherein the exposed solvent-accessible surface areas of proteins tend to adopt a more extended and flexible conformation (Marsh & Teichmann, 2011; Gao et al., 2010; Zhang et al., 2009).

The flexibility of static states and flexible motions can be further considered to predict the flexibility changes of local structures with accuracy in different solvent environments (Joo et al., 2011).

Residue accessible surfaces in structural conformation

Structural movements are a feature that influence the distribution of hydrophobicity of residues. A significant distribution of the hydrophobic surface residues in solvent-exposed positions occurred in all solvent simulations (Fig. 5 and Fig. S5). Factors other than solvents can induce the exposure of hydrophobic residues. Some highly active hydrophobic sites have a tendency for organic solvent denaturation (Chakravorty et al., 2012). T1 lipase crystallized at a lower temperature (20 °C) has its hydrophobic residues tightly packed against each other and profoundly influencing their hydrophobic effects. Organic solvents and temperature-induced fluctuation can have a weakening effect on the hydrophobic core of a protein structure. By the same reasons, organic solvent mixtures and higher temperature simulations induced a dramatic exposure of buried hydrophobic residues. This can significantly allow non-conservative hydrophobic surface residues with higher RMSD and higher hydrophobic area around the hydrophobic patches other than the lid domain to be computationally redesigned to lower the RMSD and hydrophobic area and enhance stability (Park et al., 2013). Substitution of hydrophobic for hydrophilic residues has been demonstrated as a strategy for enhancing protein stability (Pedone et al., 2001; Strub et al., 2004).

Lipases show a hydrophobic mobile ‘lid’ controlling access to the buried active site, which facilitates their ‘interfacial activation’. In a classical solvent-protein interaction, the lid domain will move away upon contact with solvent molecules, exposing the buried hydrophobic residues to the protein’s surface. The increased hydrophobic area of hydrophobic residues spanning the lid domain (Fig. 6; Table 3) caused a displacement of the lid helical segments in some solvent mixtures. Since conformational changes in proteins towards an activated form is often associated with interfacial activation (Grochulski et al., 1994), hydrophobic residues at the lid are proposed to largely be involved in the conformational changes observed in some solvent mixtures, which greatly expands the nonpolar side of the lid (Fig. 7 & Fig. S6). MD simulations of other lipases indicate that greater conformational change in methanol can destroy the water layer, breaking the hydrophobic packing on the surface and exposing the residues to solvents, which also weakens hydrophobic interactions (Jiang et al., 2014). The exposure of these residues at the lid sometimes can enhance hydrophobic interactions between the enzyme and the solvents as well as interfacial activation (Brzozowski et al., 1991). The high hydrophobic area of Pro^217,198 (Table 3) will not affect the secondary structure of the protein. This is ultimately due to its preferential location in the accessible turns of proteins and the profound local impacts that a proline can impart to a structure (Lins, Thomas & Brasseur, 2003; Cherukuvada et al., 2005). In addition, its nitrogen backbone is not available for the hydrogen bonding required for helix formation, which renders it immutable to solvent effects as observed in the native T1 structure.

Solvent strain effect on conformational stability

The role of organic solvents on lipase-catalysed reactions in liquid organic media can be understood based on the competitive inhibitory character of solvents, change in conformation and substrate solubility (Graber et al., 2007). Experimentally, the effect of five polar solvents in terms of the measured residual activity had a marked effect on both enzyme stability and activity. If the lipase stability in the presence of polar organic solvents is considered, the destabilizing effect of these solvents is ultimately related to their decreasing polarity and increasing chain length (Fig. 8). Following molecular simulations in the same solvent conditions, analysis of the last structures revealed structural movements of the lid helices involved in conformational rearrangements of some secondary structure elements of the α6-loop-α7 helix (Fig. 7 & Fig. S6). The movements are more significant in all solvent mixtures in contrast to the degree of movement in the H₂O solvent. Comparisons of the resulting conformational changes of the lid show that the fluctuations seem to trigger the displacement of the lid domain and also correlate with the rms deviations and increased distance between residue Asp175 and Arg230 of the lid (Figs. 9A, 9B & Figs. S7, S8). It is not surprising to observe a gradual displacement of the helical conformation of the lid domain in MtOH-H₂O solvent mixtures because methanol has also been shown to favour induced conformational transitions in a number of protein folding studies (Babu, Moradian & Douglas, 2001). In fact, these transitions have long been demonstrated in other polar solvents in which they are implicated in the disruption of rigid tertiary structures in protein molecules with the expansion of the helical conformation (Uversky et al., 1997). Although hydrocarbons have not been reported to inhibit enzyme activity experimentally, simulations of T1 lipase in water-octane mixtures strongly suggested otherwise. Large structural rearrangement of the lid domain was caused by the interaction between the hydrophobic residues of the lid with octane. However, removal of octane from the system restored the enzyme conformational state (Graber et al., 2007; Abdul Rahman et al., 2012). High activation energy barriers also control the conformational transition stability of enzymes in organic solvents, whereby the polarity of solvents seems to impact the structural integrity of proteins in different ways (Kovacs et al., 1995; Liu & Hsu, 2003).

Network of interactions and contributions to stability

The three-dimensional structure of a protein generally achieves a certain degree of packing and stability due in part to complementary weak interactions. Hydrogen bonding, electrostatic and hydrophobic interactions, and salt bridges all play roles in protein folding and establishing its final structure (Serdakowski & Dordick, 2008; Gallivan & Dougherty, 1999). Polar solvents are strongly competitive with intramolecular hydrogen bonds and are a major cause of disruption to structural and dynamical interactions between protein atoms (Knubovets, Osterhout & Klibanov, 1999). Interactions were analysed in the presence of organic solvent mixtures (Fig. 10A & Fig. S9). It is worth noting that the role of hydrogen bonds in terms of folding and stability of proteins does not necessarily lead to any ordered structural classification and thus a systematic investigation could be difficult (Sanjeev & Vishveshwara, 2002). The analysed hydrogen bonds revealed differences between the numbers of additional hydrogen bonds in all solvents. The hydrogen bonds remained well preserved in BtOH-H₂O-H₂O, EtOH- H₂O and PrOH-H₂O in the last 10 ns and decreased for BtOH-H₂O and PtOH- H₂O solvents. The decreased number of hydrogen bonds was not in response to the solvent polarity throughout the simulations. An ordered secondary structure can only be obtained where there is a large hydrogen bond formation. The introduction of stabilizing internal bonds also minimizes the solvent strain effect limiting aggregation and irreversible inactivation of an unfolded protein in a nonaqueous environment (Arnold, 1990). Hydrogen bonds and other interactions are characterized by the reorientation of the polar amino acid side chains making the protein rigid and compact in an organic solvent (Khabiri et al., 2013; Yang, Dordick & Garde, 2004). Disruption of these interactions could lead to reduced flexibility of residues. This is often due to the decreasing water-mediated network resulting from the exchange of water molecules with solvent molecules at the protein surface (Trodler & Pleiss, 2008; Zheng & Ornstein, 1996).

The contribution of hydrophobic interactions depends on how many nonpolar groups are buried during folding of a protein (Pace et al., 2012). Hydrophobic interactions stabilizing the protein structure remain weakened with increasing carbon atoms of polar solvents (Table 4). The reduced amount of interactions highlights the general role of hydrophobic interactions between a target protein and the surrounding environment. The hydrophobic and hydrophilic properties of the organic solvent in a protein environment as well as exposed surface residues of the protein govern the type of interactions that are formed between protein and solvent molecules (Patargias, Harris & Harding, 2010). Previous results from MD simulations have also elaborated on how methanol strengthens hydrogen bonding and weakens hydrophobic interactions in proteins (Hwang et al., 2011).

Geometric orientation changes in active site residues, Ca²⁺ and Zn²⁺ as it affects protein stability

The overall rms deviation of the structure does not capture the details of any implied global conformational changes to the active site or metal cofactors. Therefore, the distances of the proximal residues holding the protein metal ions in place were analysed (Table 5). Distances and bond angles between residues differed slightly from the initial structure 2DSN. The metal ions and bond angle appear to drift from their original position as shown by the increase of their distance of separation. The movements of the coordinating residues of the metal cofactors can allow accessibility to the metal binding pocket, which is also important in the interfacial activation of lipases. Molecular simulations have also revealed the geometric displacement of these metal cofactors and the coordinating residues cause major structural transitions in lipases (Invernizzi et al., 2009).

The formation of hydrogen bonds is central to the transition state in an ester synthesis mediated by the Ser113 acting as a nucleophile promoted by His358, which accepts a proton from Ser113. Although the angle of the alpha carbons of the active sites residues was not preserved in the solvents when compared with the crystallized T1 lipase (2DSN) with a bond angle of 43.16 Å, the distance between Asp317 and His358 of the active site in the last structure was preserved (Figs. 11A,11B, Fig. S10 & Table 2). However, Ser113 and His358 had more conserved distances in H₂O. The conserved distance in MtOH- H₂O and EtOH-H₂O is sufficient to an extent to allow the protonation of His358 by Ser113, which can promote binding and catalysis. This accounts for the higher activity observed for MtOH-H₂O and EtOH- H₂O in Fig. 8 in contrast to other solvents. This observation suggests that the solvents influenced the minimum distance between these residues, which is crucial in the formation of H-bonds that play significant stabilization functions in the transition state of the lipase reaction (Carlqvist et al., 2003). Similarly, widening and stability of the active site cavity mediated by the hydrogen bonds may also affect proper positioning and geometric orientation of the backbone atoms of residues, which are crucial for the effective catalytic activity in polar solvents.

Solvent distribution around proteins and the penetration effect on stability

The number of organic solvent molecules introduced into each simulation cell varied (Table 1). Water molecules will normally hydrate ionic and polar sites on the surface of the protein, which creates an exposure of the enzyme surface that is accessible by the organic solvents into the protein core (Yang, Dordick & Garde, 2004). We evaluated the consequences of these effects and found that indeed solvent molecules populate approximately the same regions of the protein surface in the last 40 ns structure (Fig. 12). Interestingly, it can be observed that the number of solvent molecules that migrated slightly down to the core of the protein were considerably higher in pentanol than other solvent simulations. The penetration of organic solvent molecules could have a strong interaction with surrounding residues. Organic solvent molecules appeared arranged in spatial proximity to the protein surface, which could strongly influence the protein stability. However, the size of the organic solvent molecules does not appear to govern solvent penetration. Clustered networks of the solvent molecules around the immediate vicinity of the enzyme-surrounding surface can elicit remarkable effects with the tendency to establish a solvation effect on the residues. Similarly, significant penetration of more PtOH solvent molecules further into close proximity to the active sites and the surrounding residues (despite its size) is governed by its decreased polarity, and will ultimately alter compositions of electrostatic interactions, increase the hydrophobic area and play a prominent role in reduced secondary structure contents. In fact, solvent penetration not only constrains the stability of the protein but also causes dramatic effects on the overall folding stability and dynamics. This is in good agreement with simulations in polar organic solvent-water mixtures in which the accumulation of ethanol molecules in the core of protein causes large conformational changes destroying its structure (Lousa, Baptista & Soares, 2012).

Secondary structure content analysis

The amino acid residues and the protein structure ultimately share a relationship that characterizes the secondary structure assignment of an entire protein. The preserved secondary structure assignment analysis is presented in Fig. 13 & Fig. S11, and the residue secondary structure structural analysis is presented in Fig. S12. There were alterations to the structure and dynamics of the T1 lipase in all environments, mostly to the helices and coils. The helix and helix 3–10 were preserved in all solvent mixtures implying preserved structural flexibility and rigidity (Khabiri et al., 2013). Variations in residue structural elements (Fig. S12) is a common effect of such polar solvents on lipases, whereby residue positions at the protein surface expose them to such rearrangements as observed in the α-helices contents of a Candida albican lipase (Li et al., 2010). It follows then that most protein structures have the propensity to undergo a conformational transition in certain organic solvent-induced environments (Hirota, Mizuno & Goto, 1997; Shiraki, Nishikawa & Goto, 1995).

Dynamic correlated and concerted motions of protein residues

The functional flexibility of proteins and the relatively concerted motions between groups of protein atoms, could determine large-scale conformational transitions and provide a wide range of protein conformations. This is important in understanding the biological functions of proteins based on the co-operative and coherent movement of amino acid residues (Arcangeli, Bizzarri & Cannistraro, 2001). We analysed the collective dynamic characteristics exhibited by the protein in the different solvent mixtures based on the DCCM, which reports the time-correlated motions of atomic pairs (Fig. 14). Generally, positive correlations involve neighbouring groups, which move together, and anti-correlated motions are of particular interest (Luo & Bruice, 2002). In all maps, region Asp175-Arg230 spanning the lid has strongly correlated movements to residues Ile107-Gln114 (which comprise the helix and sheet on which the active site Ser113 is located) and perfectly correlated to the loop residues Ala109-Glu219 (along the diagonal) connecting the α6 and α7 helix of the lid in all solvents. The helix region of Met287-Ser301 shielding the active site residues Asp317 and His358 were strongly correlated to the lid region in BtOH-H₂O and PtOH- H₂O and the region on which the active site Ser113 is located (Ile107-Gln114) in all solvent mixture systems. The region of the front part of the active site (His14-Ala52) has some weak anticorrelated motions with the lid in H₂O and correlated movements in MtOH-, EtOH-, PrOH-, and BtOH-H₂O with a strongly correlated motion in PtOH-H₂O. Strongly correlated motions between residues Ala52-Pro95 comprising two helices below the lid occurred in all solvent mixture systems. Checking residues Arg103 and Leu277 in all systems indicates they were weakly correlated in H₂O, MtOH-, EtOH- and PrOH-H₂O, but strongly correlated with surrounding residues in the BtOH- H₂O and PtOH-H₂O solvent mixture systems.

The negatively correlated motions of DCCM occurs mostly among residues that are far apart in the secondary structure of the protein. Correlation patterns within the protein are more pronounced in an approximately concerted fashion, but the correlation coefficients obtained from the DCCM method do not bear any information about the magnitude of the motions. It does, however, highlight the consistency with simulated B-factors corresponding to the increasing chain length of polar solvents.