Hydration thermodynamics of cytosolic phospholipase A₂ GIVA predict its membrane-associated parts and its highly hydrated binding site

Abstract

During biological events, the water molecules associated with the protein are re-oriented to adapt to the new conditions, inducing changes in the system’s free energy. The characterization of water structure and thermodynamics may facilitate the prediction of certain biological events, such as the binding of a ligand and the membrane-associated parts of a protein. In this computational study, we calculated the hydration thermodynamics of cytosolic phospholipase A₂ group IV (GIVA cPLA₂) to study the hydration properties of the protein’s surface and binding pocket. Hydrophobicity scales and the Grid Inhomogeneous Solvation Theory (GIST) tool were employed for the calculations. The hydrophobic areas of the protein’s surface were predicted more accurately with the GIST method rather than with the hydrophobicity scales. Based on this, a model of the protein-membrane complex was constructed. In addition, the calculation revealed the highly hydrated binding pocket that further contribute to our understanding of the ligands’ binding.

Communicated by Ramaswamy H. Sarma

Keywords:

• Cytosolic phospholipase A₂ GIVA
• hydrophobicity scales
• grid inhomogeneous solvation theory
• GIST

Introduction

Phospholipases A₂ (PLA₂s) is a super-family of enzymes that cleaves a fatty acid from the sn-2 position of phospholipids (Dennis, Cao, Hsu, Magrioti, & Kokotos, 2011). The products of this enzymatic reaction are metabolized through several biological paths (e.g. COX, LOX), promoting the synthesis of inflammatory mediators. The exact conditions that promote the upregulation of these enzymes and the interactions between the isoforms are still under investigation. The selective inhibition of PLA₂s has long been a promising approach for the elucidation of their contribution to the mechanism that leads to inflammation and the discovery of an effective treatment. However, both the high homology of the enzymes belonging to the same group and the similarities of their catalytic mechanism, make this target challenging.

In this work, we focused on group IVA of cytosolic phospholipase A₂ (GIVA cPLA₂) (Leslie, 2015). This enzyme is the most well studied among the cytosolic PLA₂s and its Ca²⁺-dependence for activation and substrate preference for arachidonic acid (Dennis et al., 2011) characterizes it. The enzyme is located in the cytosol. There is a translocation of the enzyme close to the cell membrane upon its activation by two calcium ions, which bind in the sub-domain of the enzyme. The translocation is necessary for the enzyme to reach the substrate and therefore its association with the cell membrane is a crucial step for its catalytic activity. Hydrogen/deuterium exchange mass spectrometry (DXMS) has been employed in the past for the description of the GIVA cPLA₂ orientation on the membrane and the interactions formed in the interfacial binding region (Mouchlis, Bucher, McCammon, & Dennis, 2015). However, this type of analysis is not available for every membrane-associated enzyme and theoretical calculations are usually employed in order to predict the parts of proteins that interact with the cell membrane.

Known GIVA cPLA₂ inhibitors are characterized by a diversity of the chemical groups they belong to, as it has been described in a number of reviews (Kokotou, Limnios, Nikolaou, Psarra, & Kokotos, 2017; Nikolaou, Kokotou, Vasilakaki, & Kokotos, 2019). More specifically, pyrrolidine-based inhibitors (Seno et al., 2000) and indole-based inhibitors (Lee et al., 2007; McKew et al., 2006) are potent inhibitors of GIVA cPLA₂ and some of them advanced to clinical trials (http://ClinicalTrials Identifier: NCT00396955 (2006) and NCT02795832 (2016)). In addition, novel inhibitors belonging to indole series have been designed and studied by Tomoo et al. (2014) and demonstrate some promising results. Activated carbonyl-based indol-1-yl-propan-2-ones inhibitors (Bovens et al., 2010; Ludwig, Bovens, Brauch, Elfringhoff, & Lehr, 2006) have also been systematically studied by Lehr and co-workers (Bovens et al., 2010), while Kokotos and co-workers have developed a series of 2-oxoamides derivatives (Kokotos et al., 2004; Six et al., 2007) thiazolyl ketones (Kokotos et al., 2014) and 2-oxoesters (Kokotou, Galiatsatou, et al., 2017), which demonstrate potent and selective inhibition against PLA₂s. However, the development of optimized lead compounds for the inhibition of GIVA cPLA₂ action is still an active field of research.

The only available crystallographic structure of GIVA cPLA₂ is of the apo form of the enzyme (Dessen et al., 1999). Due to the lack of any co-crystalized X-ray ligand structure, the binding mode of its inhibitors is unknown, which make the design of new inhibitors challenging. Additionally, the lid region is closed in the apo form and the binding mode of few inhibitors has been defined using DXMS, a method that does not provide details about the conformation of the ligand (Burke et al., 2009). Therefore, only computational simulations (Burke et al., 2008) may describe the structural rearrangement of the binding pocket upon ligand’s insertion.

Both events, the GIVA cPLA₂’s binding on the cell membrane and the formation of the ligand-protein complex, are partly driven by hydration effects. Here, we have addressed both events by computationally calculating the hydration properties of the protein’s surface and binding region. We aim to understand the role of hydrophilic and hydrophobic forces lead to the assembly of the protein-membrane complex and the hydration properties of the binding pocket.

Hydrophobicity scales have been employed in the past to assign hydrophobicity values to amino acids and predict the hydrophobic domains of proteins. Various hydrophobicity scales have been published (Cornette et al., 1987; Eisenberg, Schwarz, Komaromy, & Wall, 1984; Engelman, Steitz, & Goldman, 1986; Hopp & Woods, 1983; Janin, 1979; Kyte & Doolittle, 1982; Rose, Geselowitz, Lesser, Lee, & Zehfus, 1985), but the complete range of an amino acid behavior is difficult to be described by one parameter, and therefore the results are discordant. Lately, progress has been made to improve the accuracy of such scales by using modern computational methods (Nicolau, Paszek, Fulga, & Nicolau, 2014; Schauperl, Podewitz, Waldner, & Liedl, 2016). Increasing the resolution by replacing the amino acid values by atomic values is one parameter that has been introduced in the development of these scales. In addition, the environment of each amino acid, which includes other amino acids and their effect on the surrounding water molecules, has also been considered.

In this work, we assigned hydrophobicity values on the amino acids based on the calculation of the thermodynamic properties of the surrounding water molecules. We implemented this using the GIST (Nguyen, Cruz, Gilson, & Kurtzman, 2014; Nguyen, Gilson, & Young, 2011; Nguyen, Kurtzman, & Gilson, 2012; Raman & MacKerell, 2013) tool in AmberTools, an application of the inhomogeneous fluid solvation theory (IFST) on biomolecules (Green, 1952; Huggins, 2012; Lazaridis, 1998; Li & Lazaridis, 2006). GIST divides a defined solvent area in a 3D grid, which consists of small grid boxes. It discretizes the equations of IFST to approximate the localized thermodynamic properties of these small grid boxes. The calculation of the localized thermodynamic properties, including entropy and enthalpy, is an important feature of GIST.

Herein, we report the application of hydrophobicity scales and a method based on GIST analysis to map the hydrophobic areas of GIVA cPLA₂ surface in order to predict the hydrophobic areas of the protein’s surface that may interact with the membrane. We also describe the thermodynamic properties of the solvent areas in the active site of GIVA cPLA₂ that may affect ligands’ binding.

Computational methods

The only available crystallographic structure of GIVA cPLA₂ (PDB ID: 1CJY) (Dessen et al., 1999) was used for the calculations. The missing amino acids of the apo structure were added in Chimera (Pettersen et al., 2004), using Modeller 9.17 (Sali & Blundell, 1993). The charmm-gui on-line tool was used for building the lipid bilayer of 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) and placing the protein on the lipids surface (Figure 1), according to a previously published model created on the basis of DXMS data (Burke et al., 2009).

Figure 1. GIVA cPLA₂ on a lipid bilayer of POPC generated with the charmm-gui tool.

All atoms unrestrained MD simulations in explicitly solvent were performed with AMBER (Case et al., 2016). The one system we simulated was the apo GIVA cPLA₂ in solvent (named ‘the apo system’ in this paper) and the other was the GIVA cPLA₂-membrane complex in solvent (named ‘the membrane system’). The crystal water molecules and the two calcium ions were retained, because calcium ions keep the waters ordered and create interactions with the membrane, which contribute to the formation of the complex. In addition, a water molecule close to the catalytic dyad has been resolved in the X-ray structure. This water may be involved in the catalytic mechanism and therefore it was important for our calculations to keep it. tLEaP module and lipid14 (Dickson et al., 2014) and ff14SB (Maier et al., 2015) force fields were used to produce the parameter files. Systems were solvated in tLEaP with approximately 37,000 (the apo system) and 62,000 (the membrane system) water molecules. The TIP3P water model was used to describe the properties of the solvent. A truncated octahedral water box with 12.0 Å periodic boundary conditions and a nonbonded cutoff distance of 8 Å were used for the apo system. The solvation of the membrane system was accomplished with the AddToBox program and SetBox command in tLEaP to create a periodic system of {133.769 134.537 178.282} dimensions in Å. The nonbonded cutoff distance was set at 10 Å. Both systems were neutralized by adding 28 in the apo and 36 in the membrane system Na⁺ ions using the Joung and Cheatham ion parameters (Joung & Cheatham, 2008).

Consequently, we performed minimization with SANDER (Salomon-Ferrer, Case, & Walker, 2012). For the apo system, a four-step energy minimization workflow including the steepest descent algorithm followed by the conjugate gradient method was set. During the first step, the solute was kept fixed with a harmonic force constant of 500 kcal·mol⁻¹·Å⁻¹. We applied two additional minimization steps during of which the strength of the restraint was gradually reduced. For the membrane system, a four-step energy minimization workflow was also set: 10,000 steps of the steepest descent algorithm and 5,000 steps of conjugate gradient with the solute constrained and 40,000 steps of the steepest descent algorithm and 10,000 steps of conjugate gradient unconstrained. The solute was kept fixed with a harmonic force constant of 500 kcal·mol⁻¹·Å⁻¹. Minimization is necessary to reduce any bad contacts and allow the solvent and sidechains to reach an optimal orientation.

The systems were heated from 0 to 310 K in three steps under constant volume and with a reduced restraint of 10 kcal·mol⁻¹·Å⁻¹. The SHAKE algorithm was applied and the Langevin thermostat with a collision frequency of γ = 2.0 ps⁻¹ for the apo system and 1.0 ps⁻¹ for the membrane system was used. Equilibration was performed in the isothermal isobaric (NPT) ensemble at 310 K and 1 atm pressure. For the apo system, the equilibration was performed in two steps: one with the solute restrained and only the water molecules allowed to move and one with the entire system unrestrained. The total time of equilibration was 200 ps and the Berendsen barostat was utilized for pressure control with isotropic position scaling (ntp = 1). For the membrane system, equilibration was performed for 1.2 ns in four sequential steps, keeping the protein and lipid atoms fixed and 5 ns keeping only the protein fixed. The Monte Carlo barostat was utilized for pressure control with anisotropic pressure scaling (ntp = 2). The long-range electrostatic interactions were accounted for using the particle mesh Ewald (PME) method (Darden, York, & Pedersen, 1993). The pmemd.cuda (Salomon-Ferrer, Götz, Poole, Le Grand, & Walker, 2013) engine of AMBER was used for the heating, equilibration and production time. The total time of the unconstrained MD simulation was 300 ns.

The MD trajectories were clustered and the top three clusters were further proceeded for a GIST analysis. We repeated the GIST analysis three times to allow a better sampling of the analysis and increased the accuracy of the results. Clustering was performed using AmberTools. The bottom-up clustering method was used and distance between frames calculated via best-fit coordinate RMSD using residues 182-186, 214, 218, 248-251, 253, 277, 281, 382-384, 387, 390, 403, 407, 535-537, 541, 564, 664 and 666-670, which define the lid region. Each cluster was then prepared and simulated for 20 ns keeping the solute constrained (harmonic force constant of 10 kcal·mol⁻¹·Å⁻¹) and allowing only the waters to move. It is suggested (Ramsey et al., 2016) to keep the solute of interest restrained in the lab frame, according to the source particle method used in IST (Lazaridis, 1998). 10,000 simulation snapshots were generated from each simulation. The GIST analysis was performed with AmberTools17 and the grid spacing was set to 0.5 Å. The grid dimensions along each coordinate axis were set to 150 Å to cover the whole protein surface. The grid points within 5 Å distance from the solute were counted in for determining the thermodynamic properties of the hydrating water molecules. The gistpp post-processing tool (Ramsey et al., 2016) and in-house scripts were used for analyzing the GIST output files. For the calculation of the localized values for the free energy, the ΔG values of the grid points around each atom in a distance of 5 Å were integrated and normalized by the number of water molecules at each grid point.

Results

Hydrophobicity scales and GIST analysis of GIVA cPLA₂ surface

Hydrophobicity scales have been used in the past to identify the hydrophobic regions of transmembrane proteins that interact with the membrane (Eisenberg et al., 1984). Here, we employed well known hydrophobicity scales as well as a method build on the GIST tool to calculate these regions of GIVA cPLA₂, based on amino acids’ and water molecules’ thermodynamic properties, respectively. Then, we compared our results with the experimental data to justify the most accurate computational method for predicting the membrane-associated part of the protein.

The hydrophobic surface of GIVA cPLA₂ generated by various hydrophobicity scales is presented in Figure 2a–i by a colour-coded representation. The results generated by the GIST analysis of the apo system are shown in Figure 2k. Figure 3 presents the orientation of the protein on the membrane according to DXMS data.

Figure 2. Amino acids color-coded (red hydrophobic; blue hydrophilic) based on the values of hydrophobicity scales of: (a) Cornette et al. (1987) (−3.1 to 5.7); (b) Eisenberg et al. (1984) (−2.53 to 1.38); (c) Engelman et al. (1986) (−12.3 to 3.7); (d) Hopp and Woods (1983) (−3.4 to 3.0); (e) Janin (1979) (−1.8 to 0.9); (f) Kyte & Doolittle (1982) (−4.5 kcal/mol to 4.5 kcal/mol); (g) Rose et al. (1985) (0.52 to 0.91); (i) Schauperl et al. (2016) (−15 to −5) and (k) the ΔG values of the surrounding water molecules (cut-off 4 Å) calculated using GIST tool (−4kcal/mol to 1 kcal/mol).

Figure 3. Orientation of GIVA cPLA₂ in a lipid bilayer according to DXMS data. The description of the hydrophobic areas (red) of the protein surface was produced by the GIST tool. These areas penetrate the membrane and this is in accordance with the DXMS data.

According to DXMS data, the areas of the enzyme which penetrate the cell membrane are addressed by residues 35 to 39 and 96 to 98, belonging to C2 domain, including the residues around the entrance of the tunnel which leads to the catalytic dyad. Based on our GIST calculations, these areas were predicted to be hydrophobic (colored in red, in Figure 3) and so expected to penetrate the membrane and interact with the aliphatic parts of it.

To understand further these results, we studied the interactions formed between the protein and the lipids during the MD simulation of the complex. Based on the simulation, residues at the entrance of the channel form interactions with the aliphatic chains and the polar head groups of the lipids membrane (depending on the type of the amino acids), which stabilize the protein-membrane association (Figure 4). We think that the combination of the GIST results with the MD simulations may be harnessed to identify regions in proteins that are more likely to bind to a membrane than other methods.

Figure 4. A representative snap shot from the MD simulation of the protein-membrane complex. The entrance of the tunnel (yellow residues) is highlighted with a yellow circle. It leads to the catalytic dyad Ser228/Asp549 shown by the orange arrow. Met468 and Trp464 (yellow) interact with the lipophilic part of lipids (cyan) and Arg467, Lys273 and Arg274 (yellow) with the polar head groups, in accordance with the DXMS data (only some parts of the lipids are represented here (cyan ‘stick’) for clarity reasons).

GIST analysis of the solvent areas in the binding pocket of GIVA cPLA₂

The binding pocket of GIVA cPLA₂ has the shape of a tunnel, which starts at the protein’s surface and ends at the catalytic dyad. The number of water molecules there is limited, mainly because the tunnel is narrow and consists of lipophilic amino acids. We focused our study on the hydration properties of the binding pocket close to the catalytic dyad. As notice above, only one water molecule was resolved close to the catalytic dyad by the X-ray structure.

The MD simulations revealed that during the simulation time a well-organized network of waters is formed around the catalytic dyad Ser228/Asp549 (Figure 5 and Figure S4 (Supplementary material)). It is worth to mention here that the formation of this network was observed in both systems, the apo system and the membrane system. Two more water molecules are located close to Arg200, which according to the catalytic mechanism stabilizes the intermediate product of the reaction.

Figure 5. The network of water molecules around the catalytic dyad Ser228/Asp549 of GIVA cPLA₂ formed during a 300 ns MD simulation.

To investigate further the water occupancy of the binding pocket, we calculated the water occupancy around the catalytic serine within 6.5 Å using the watershell command in the AmberTools. In the graphical representation (Supplementary material, S1), the area around the catalytic serine is been occupied by 6 to 7 water molecules for most of the simulation time. These waters make interactions with the catalytic dyad of Ser/Asp and the amino acids Gly551, Asn555, Gly198, Trp393, Thr231, Gly226 and Thr330. The high occupancy map (calculated with GIST) has identified the locations of these hydrophilic centers. Using the ‘closest’ command in AmberTools we obtained the 10 closest residues to the catalytic serine (Supplementary material, S2). The results relieved that the water molecules persist in this area for a considerable simulation time. They do not exchange with the solvent frequently due to the interactions create with the amino acids of the binding pocket.

GIST analysis was applied on the binding pocket of GIVA cPLA₂ in the apo system. For the GIST analysis, three clusters, produced by the 300 ns simulation time, were used. ΔG values were calculated based on the GIST analysis: (1) $${\mathrm{\Delta }G}^{\mathit{norm}}={\mathrm{\Delta }E}_{\mathit{tot}}^{\mathit{norm}}-\mathrm{ }{\mathrm{T}\mathrm{\Delta }S}_{\mathit{six}}^{\mathit{norm}}$$(1) where the ${\mathrm{\Delta }E}_{\mathit{tot}}^{\mathit{norm}}$ is the sum of the $E_{ww}^{\mathit{norm}}$ (water-water enthalpy referenced to the bulk) and $E_{sw}^{\mathit{norm}}$ (solute-water enthalpy referenced to the bulk). The ${\mathrm{T}\mathrm{\Delta }S}_{\mathit{six}}^{\mathit{norm}}$ term represents the total entropy. The calculated values of ΔG around the binding pocket are represented as a grid in Figure 6. The positions of the water molecules are from a representative structure of the most populated cluster.

Figure 6. The $\mathrm{\Delta }{G}^{\mathit{norm}}$ values are presented as orange grid, with an energy cutoff of −5 kcal/mol.

The energy calculations could be used to identify the less favorably bound water molecules. They could be replaced by a part of the ligand, towards a higher binding affinity, as it has been addressed elsewhere (Abel, Young, Farid, Berne, & Friesner, 2008; Abel et al., 2011). The hydration properties of a binding pocket may be useful in identifying key features that a drug molecule could mimic when binding its target. Therefore, the structure of the ligand should be suitable to replace specific water molecules and mimic their interactions with the key amino acids.

Conclusions

Hydration properties of proteins’ surface play a pivotal role in protein-protein and protein-membrane interactions. Hydration scales have been used in the past in an attempt to predict hydration properties of proteins’ surface, based on the amino acid’s properties. The available hydrophobicity scales assign a fixed value for every type of amino acid, neglecting their individual environments. Hence, these scales fail to describe accurately the energetic contributors of hydration. To accurately describe the hydrophobic character of the individual amino acids, the thermodynamic properties of water molecules around every amino acid have to be investigated including also the influence of local environment e.g. other molecules or amino acids.

In order to account the local difference, we calculated the hydrophobic character per atom of every amino acid using GIST analysis and MD simulations. We apply this method on the surface of the apo structure of GIVA cPLA₂ in explicit solvent, in order to identify the hydrophobic points of the surface. We concluded that the prediction of GIVA cPLA₂ surface parts, which penetrates the membrane, better meets the previously experimental DXMS data, when made based on the GIST analysis and MD simulations rather than other hydrophobicity scales. The main difference between the methods used by the hydrophobic scales and GIST tool for the calculations is that GIST includes the influence of the surrounding molecules.

GIST uses the discrete sums of the solvation energy and entropy quantities over the voxels of a 3 D grid, in order to estimate a local free energy of solvation for the voxels. The published hydrophobicity scales, referenced in Figure 2, use a variety of approaches to calculate the amino acids properties. Our results show that by including the solvation energy and entropy terms, assisted by MD simulations, a quantified description of the environment of the amino acids is possible.

Based on the MD simulations of the GIVA cPLA₂ – membrane complex, interactions between the lipids and the surface amino acids are formed and stabilize the protein’s orientation in the membrane. These interactions may be used as a target to inhibit enzyme’s activity by inhibiting its attachment to the membrane.

In addition, we analyzed the hydration properties of the binding pocket in an attempt to understand better ligands’ binding. The MD simulation of the apo structure revealed a highly solvated catalytic center. In addition, the tunnel that leads to the catalytic center of the enzyme is narrow and it consist of, mainly, lipophilic amino acids. Given the structure of the natural substrate (membrane phospholipids), we suggest that this highly hydrated area around the catalytic center at the bottom of the tunnel, promote the right placement of the ligand near the catalytic serine via its polar groups. Although the favorable ΔG values (cutoff of −5 kcal/mol) of the waters around this area may suggest that the replacement of these waters is not promoted, we know – based on the catalytic mechanism and experimental data – that this area may be occupied by a ligand upon its insertion (Dessen et al., 1999; Kokotou, Galiatsatou, et al., 2017). Thus, we suggest that the replacement of certain water molecules should occur upon ligand’s insertion, when a favorable interaction between the ligand and the amino acids is formed.

In conclusion, hydration thermodynamics assist by MD simulations may be used to predict the interactions in the interface between enzyme-membrane complexes. In the case of GIVA cPLA₂, the predicted protein-membrane model based on the hydrophobicity properties of its surface, calculated by GIST and MD simulations, was in a good agreement with the DXMS data. Thus, we could suggest this application for the prediction of the interactions in a protein-membrane complex, in cases where no experimental data are available.

Moreover, we analyzed the hydration thermodynamics of the catalytic site of GIVA cPLA₂. The area, where the catalytic mechanism takes place, seems to be highly hydrated in contrast with the rest of the binding pocket, which consist of lipophilic amino acids. This event may act as a driving force for the right placement of the natural substrate close to the catalytic serine to promote the hydrolysis. The less favorably bound water molecules may be replaced by a part of a ligand and therefore these results may be used in the design of new GIVA cPLA₂ inhibitors.

Abbreviations
COX	=	cyclooxygenase
DXMS	=	Hydrogen/deuterium exchange mass spectrometry
GIVA cPLA2	=	cPLA2 cytosolic phospholipase A2 group IVA
GIST	=	Grid Inhomogeneous Solvation Theory
LOX	=	lipoxygenase
MD	=	molecular dynamics
POPC	=	1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported by computational time granted from the Greek Research & Technology Network (GRNET) in the National HPC facility – ARIS – under project ID: pr002030-MODYCPLA2. S.V. was supported by the EMBO Short-term fellowship. This work was supported by the Austrian Science Fund (FWF) [project P 30565].