A Computational Study of Some Cancer Causing S100 Proteins Interacting with A Ligand, Parecoxib and Investigation of Drug Likeness Potential of the Ligand

List Of Abbreviations

GROMACS = GROningen MAchine for Chemical Simulations

AMBER = Assisted Model Building with Energy Refinement

CHARMM = Chemistry at Harvard Macromolecular Mechanics

GROMOS = Groningen Molecular Mechanics

MMFF = The Marc Molecular Force Field

OPLS = Optimized Potential for Liquid Simulations

SPC = Single Point Charge

TIP3P = Transferable Intermolecular Potential 3 points

Introduction

Despite promising advancement in the field of Research and Development of drug development, Machine learning is an excellent way of predicting a molecule as a drug reducing both time and cost of research. Protein–ligand interactions are determined by electrostatic forces, hydrogen bonding, van der Waal’s interactions, and hydrophobic effects, collectively determining binding specificity and affinity [1–3]. Proteins belonging to the S100 family are low-molecular-weight calcium-binding proteins involved in diverse intracellular and extracellular regulatory functions. Dysregulated expression of several S100 proteins has been implicated in cancer, inflammation, neurological disorders, and immune modulation [4 - 6]. Many S100 proteins are active extracellularly and are easily detectable in blood or other body fluids. So, search for effective inhibitors that disturb S100 protein functions and studies on the structure and conformation of S100-ligand may leads to designing S100-targeted therapies [7-9].

Studying big molecules like polymers, proteins, nucleic acids etc. with pure quantum chemistry is impossible. In such cases theories of approximations and iterative approaches are applied. Molecular Dynamics (MD) computes the structure and energy of molecule based on molecular motions. In this method electrons are considered explicitly and assumed that they will find their optimum distribution once the position of the nuclei is known [3]. MD describes the potential energy of the system using mathematical equation termed as force fields. The potential energy is based on two types of forces – the bonded and the nonbonded interaction energies. Total energy is contributed from – bond stretching, bond angles, bond torsion, van der Waal’s energy and electrostatics energy [10,11].

Docking predicts binding orientation and affinity, while MD simulations capture the dynamic stability of protein–ligand complexes under physiological conditions [12,13]. Binding free energies estimated using Molecular Mechanics Poisson–Boltzmann Surface Area (MM-PBSA) methods correlate well with experimental observations [14,15]. The relative binding energy of the ligand with protein to form a protein-ligand complex (Figure 2) is given by the expression

Protein–ligand binding is a dynamic equilibrium process where ligands often bind to pre-existing conformations of proteins, sometimes inducing moderate structural changes [3,4]. Binding strengths vary widely, typically ranging from −2 to −12 kcal/mol, and are governed by enthalpic and entropic contributions [4]. Accurate determination of binding free energy (ΔG) is therefore essential for understanding molecular recognition and drug efficacy.

Molecular dynamics simulations rely on classical mechanics and force-field models such as AMBER, CHARMM, GROMOS, MMFF, and OPLS to describe bonded and non-bonded interactions within biomolecules [12,13]. Solvent effects are incorporated using explicit water models (SPC, TIP3P, TIP4P) or implicit continuum approaches [13,14]. In this study we utilise CHARMM-27 force field and SPC as water model [18].

Molecular docking is widely employed in drug discovery to predict ligand binding modes and affinities within protein active sites. AutoDock Vina (AutoDockTools-1.5.7) is an open-source docking tool known for its speed and accuracy [19,20]. Ligand libraries are commonly obtained from databases such as ZINC, PubChem, or ChemSpider [21,22]. We use ZINC database here (www.zinc.docking.org/).

Parecoxib is a naturally occurring compound and included in the class of isoxazole compounds. Isoxazole compounds are effective in cancer treatment and potential possibilities [23]. These are heterocyclic compound consisting of a five-membered ring with one oxygen atom and one nitrogen atom in the neighbouring position. Docking of a ligand into a selected binding site of a protein is carried out with a docking tool Autodock vina working on definite algorithm [24]. This tool calculates with different poses of the ligand and list out the result starting with highest binding affinity. Binding affinity is not sufficient parameter to predict drug efficacy. For that we go for Molecular Dynamics (http:// www.mdtutorials.com/gmx/complex/index.html)simulations and MMPBSA (gmx_MMPBSA 1.6.4) calculation to obtain root mean square deviation (RMSD) [25], H-Bond, root mean square fluctuation (RMSF), solvent accessible surface area (SASA), radius of gyration (RoG) and total free energy change during non¬bonded protein-ligand complex formation [26-28] (Figure 2).

In this study we select Parecoxib as a ligand and five different S100 proteins obtained from Protein data bank (https://www.rcsb. org/) are allowed to dock separately. S100 proteins are S100A6, S100A2, S100A4, S100A7 and S100A11. On the basis of affinity energy, we select the best pose of the ligand from each docking. Molecular Dynamic Simulation of each complex are carried out for 100 ns or less than 100 ns. We go for MMPBSA analysis for thermodynamic calculations to obtain free energy change in protein-ligand complex formation [16].

Materials and Methods

Molecular Docking is an approach used to predict how a ligand (small molecule) binds and orient (the 3D conformation) with a protein (a receptor) using computational method, AutoDock Vina [7]. By docking, physical interaction between the ligand and protein are calculated. The output of this docking is a docking score, which estimates the strength of binding between the molecule (the ligand) and protein.

Here Parecoxib is our ligand. We select cancer causing five different S100 proteins, such as S100A6, S100A2, S100A4, S100A7 and S100A11. Individual PDB files are taken from protein databank, RSCB.org.

Ligand Preparation: The ligand Parecoxib is search on ZINC database (Table 1). We record the ID of the compound and download the SDF file. This SDF file is then open on Biovia Discovery Studio (https://www.3ds.com/products/biovia). Now the corresponding PDF file is generated. This file is open on AutoDock tool. PDBQT file is prepared here. Now this file can be used for docking [10].

Protein Preparation: Protein PDB file is obtained from Protein Data bank (RSCB.org). Now the file is cleaned using ChimeraX (UCSF ChimeraX version: 1.9) software. Water molecules, ions, metals and sub-standard molecules are removed. This file is saved and opened on AutoDock tool. Polar hydrogen and Charges are added. After that this file is saved as PDBQT file.

Grid box preparation: Grid box is the allowed region of docking which is drawn opening the clean protein file on AutoDock tool. A box around the protein molecule is considered by defining the center and the distances of XYZ axes. With this dimension a conf. txt file is prepared.

Docking: With the help of PDBQT file of Ligand and Protein and grid box information docking can be performed using python base command prompt with the help of AutoDock Vina software. We get result of docking with 9 different poses and energies. We select the complex with lowest affinity energy for MD simulation [19,20].

MD simulation spanning up to 100 ns for protein along with their lowest energy ligand molecules are carried out using GROMACS–2024.5 software. On GROMACS topology files for proteins are generated using CHARMM-27 force field. The topology files of the ligand molecules are generated with the help of SwissParam web site (https://swissparam.ch/). Here MOL2 file of ligand molecule is required to generate topology file. Protein along with ligand molecules inserted in cubical box and solvated with simple point charged (SPC) water molecules followed by neutralizing the system keeping pH at 7.0 and temperature at 300K [10]. Subsequently, energy minimization is performed utilizing the CHARMM-27 force field, two types of thermal equilibrium conditions such as isochoric (NVT) and isobaric (NPT) ensembles are considered for energy minimisation. Long range electrostatic interaction and non-bonded interaction cut off is set at 10Å [26,27]. The MD simulation of the equilibrated system is run without restraining the solute atoms for 100ns. The analysis of MD simulation of the protein-ligand system is performed on the basis of structural perturbation and mobilities estimated in terms of Root Mean Square Deviation (RMSD), Root Mean Square Fluctuation (RMSF), H-Bond, Solvent Assessable Surface Area (SASA) Radius of Gyration (RoG) [17,20].

After successful completion of simulation, we have to go for MMPBSA, another type of simulation using conda environment on Ubuntu. One of the important steps is to install gmx_MMPBSA 1.6.4 software. Here we get an output file containing results related to thermal energy profile and the total free energy change in the formation of Protein-Ligand complex.

               Figure 1: Structure of Parecoxib Molecule

                                 Table 1: Some Properties of Parecoxib Molecule Obtained from ZINC Database

Table 2: Scores of S100A11 Docking with Ligand. Copied from Result of Autodock Vina Docking

Figure 2: Protein with Ligand at Highest Affinity Energy Saved from Discovery Studio Software

Results and Discussion

The affinity energy of Parecoxib ligand (Figure 1) in the formation of complex with each of the proteins are obtained from virtual screening with the help of Autodock Vina tool. Here we can comment that highest negative energy stands for strongest bonding with the protein. Of course, AutoDock Vina provide us a series of 9 different poses with increasing energies (Table 2). We consider the pose with highest negative (lowest energy) for MD simulation study. The different types of interactions such as hydrophobic, electrostatics, van der Waal’s and H-bonding are witnessed using the Biovia Discovery Studio (2025) software. After that we go for simulation and calculate RMSD (Figure 3), RMSF (Figure 4), Gyration Radius (Figure 6), No of H-bonds and SASA (Figure 5). But for the most important parameter we have to go for MMPBSA analysis to get the values of total change of Gibbs free energy (Table 4). Each parameters have specific values in the studies of drug ability score in selecting a ligand as drug. MD simulation for all complexes of the proteins with Parecoxib ligand performed for100 ns or less. RMSD of alpha carbon atoms, RMSF of all amino acid residues (or atoms), solvent accessible surface area, radius of gyration, number of Hydrogen bonds between protein and ligand and Gibbs free energy change in the formation of all the protein-ligand complexes are investigated.

For comparing these parameters among the five complexes (Table 3) we consider the average values of each parameter. these values are spread along simulation progress timeline. Only exception in case of RMSF (Figure 4) this is the average value of fluctuation of each atom of protein chain during the full period of simulation.

RMSD measures the structural deviation from a reference shape (Figure 3). Low values (<2.0Å) indicates structural stability. Root mean square deviation (RMSD) in atomic positions to quantify how much the protein has changed with respect to initial structure or how much conformation vary during the course of simulations. It can be calculated for the whole protein. The average of whole chain is considered for comparison. RMSD <2Å is favorable for complex (Protein- Ligand) stability. Here except S100A4 (Table 4) complex all other has value less than 2Å.

Average of root mean square fluctuations (RMSF) during the simulation to determine whether the simulation is physically realistic or not. RMSF measures the flexibility of individual residues. It measures the localized atom movement. Lower values have higher stability. RMSF should be less than 1.5Å for stability. Here S100A2, S100A4 and S100A7 complexes are less stable as they have higher RMSF (Table 3).

In stable protein-ligand complexes, a flat RMSD paired with a constant radius of gyration signifies a folded, stable structure, whereas synchronous spikes in both indicate overall conformational shift.

SASA measures the molecular surface touching the solvent (Figure 5). It is proportional to hydrophobic core exposure. Lower the surface area higher is the stability. S100A4 and S100A7 complexes poses moderately higher surface area.

Number of H- Bond indicates hydrogen bond formations over time (Table 3). Higher the Number of H-bonds stronger the protein-ligand interaction. Too much of H-bonding may rise problem of dissociation. Drug molecules also have to detach after a definite period of time. In the table (Table 3) it is observed that S100A4 and S100A7 complexes shows less chances of H-bonding.

Radius of gyration (RoG) measures the compactness of a molecule (Figure 6). Low values show high compact folding. A decrease in RoG indicates an increase in the packing density of the protein. S100A4 complex poses higher value hence it is not favourable (Table 3).

Change of Gibbs free energy determines the process spontaneity and system stability. Negative values mean spontaneity. Very high value (positive) of S100A6 complex rejects the possibility of stable complex formation (Table 4).

It is also explaining that there are eight types of intermolecular favourable interactions in case of S100A11 and Parecoxib molecule in the complex formation. They are such as two formal H- bonding, two pi-pi stacking, two pi-alkyl and two alkyl-alkyl van der Waal’s interactions (Figure 7). 2D picture is downloaded from Biovia Discovery Studio software.

Figure 3: Plots of RMSD vs Simulation Time Obtained from Output File of Simulation

Figure 4: Plots of RMSF of Each Atom Within Simulation Period Obtained from Output File of Simulation

Figure 5: Plots of SASA vs Simulation Time Obtained from Output File of Simulation

Figure 6: Plots of RoG vs Simulation Time Obtained from Output File of Simulation

Table 3: Comparison of Various Scores of Docking, MD Simulation and MMPBSA Results

SEM=Standard Error of The MeanTable 4: Differential Energies of Various Intermolecular Interactions of S100A11 and Ligand  Obtained from MMPBSA Results

Figure 7: Inter Molecular Interactions of S100A11 Protein and Ligand Obtained by Using Discovery Studio Software

Conclusion

In this study several computational tools are used to investigate the drug potentiality of Parecoxib with respect to the nature of interaction with proteins in the category of S100. Here we avoid more complex protein structures due to high computational cost. Simulation is a high computational task. Final validation of the result can be carried out preparing the stable complex followed by crystallographic studies, which is beyond the scope of this work. After establishing the drug potentiality of the ligand clinical studies and clinical trials are required to certify a compound as a save drug for treatment of disease. This process is a valid process in identifying drug probability. If molecular as well as thermodynamic properties are failed in computational method, that may be called as primary study, then we can avoid further study of that molecule as drug candidate. As a preliminary study it can be conclude that Parecoxib has the highest chances of being affective drug against S100A11 protein related cancer or other disease by acting as an inhibitor for the protein’s foul action.

References