1 Introduction

Paracetamol (PCT) or acetaminophen is one of the most commonly used analgesic and antipyretic drugs with relatively little anti-inflammatory activity [1, 2]. The analgesic and antipyretic actions of PCT resemble those of non-steroidal anti-inflammatory drugs (NSAIDs). PCT acts as weak inhibitor of prostaglandin (PG) synthesis by reducing cyclooxygenase (COX-1 and COX-2) [3, 4]. It has some demerits depending on the type and nature of unusual physical condition of the body and on the limit of dose. Overdose and long term dose of PCT produce hepatic and/or renal injury in human and other experimental animals [5,6,7]. N-acetyl-p-benzoquinone-imine (NAPQI) is the toxic intermediary metabolite of PCT that accumulate in the liver due to overdose when glutathione storage is depleted [8]. Considering molecular mechanism of analgesic activity as well as the hepatotoxicity of PCT, there are several modified derivatives are studied to search improved analgesic and antipyretic activities and to reduce side effects [6]. Different electron donating (–CH3, –OCH3, and –SCH3) and electron withdrawing (–Br, –Cl, –F, –I, –CF3 and –CCl3) substituents are placed at X, Y and Z positions of core structure in order to rationalize the analgesic and antipyretic activity.

Drug modification is another alternative way to search better agent, which can increase the selective action and reduce the side effects of drugs. Currently, it is observed that the modification of drugs inserting halogen, alkyl, alkoxy, and hydroxyl groups which play an important role in improving drug performance [9,10,11,12]. In this study, we reported the optimization of PCT and its modified derivatives to investigate their biochemical behaviour on the basis of quantum mechanical calculation. The free energy, enthalpy, dipole moment, electrostatic potential, HOMO–LUMO gap, hardness, softness, and atomic partial charge have been calculated. Molecular docking and nonbonding interactions have been performed to investigate the binding affinity, mode(s) and interactions of drugs with the amino acid residues of receptor protein. The docking prediction was further validated by dynamics simulation. Most of the derivatives showed improved thermal stability, chemical reactivity, binding affinities and interactions. Hydrogen bonding plays crucial role in the biological system which have ultimate effect on thermodynamic and structural stability of a system [13,14,15,16]. From the regarding quantum chemical studies, we are assuming that, some of the designed compounds may have profound effect as drugs.

2 Computational details

2.1 Geometry optimization

In computer aided drug design, quantum mechanical methods has greater attention on calculation of thermodynamic properties, molecular orbital features, dipole moment, atomic partial charge, molecular electrostatic potential properties [17]. Initial geometry of Paracetamol (PCT) was taken from online structure database named ChemSpider [18]. Geometry optimization and further modification of all structures were carried out using Gaussian 09 program [19]. Density functional theory (DFT) with Becke’s (B) [20] three parameter hybrid model Lee, Yang and Parr’s (LYP) correlation functional [21] under Pople’s 6-31G+(d,p) basis set was used for geometry optimization [22]. Initial optimization of all compounds was performed in the gas phase then re-optimized in water phase under polarizable continuum model (PCM) [11, 23] at same level of theory. The difference of electronic energy (∆E), enthalpy (∆H), Gibb’s free energy (∆G) and dipole moment (∆μ) were calculated by subtracting gas phase energy from water phase energy.

Molecular orbitals were calculated at same level of theory. Hardness (η) and softness (S) of all compounds were also calculated from the energies of HOMO (highest occupied molecular orbital) and LUMO (lowest unoccupied molecular orbital) considering Parr and Pearson interpretation [24, 25] of DFT and Koopmans theorem [26] on the correlation of ionization potential (I) and electron affinities (E) with HOMO and LUMO energy (). The following equations are used for the calculation of hardness (η) and softness (S):

$$\upeta = [\varepsilon {\text{LUMO}} - \varepsilon {\text{HOMO}}]/2,\quad {\text{S}} = 1/\upeta$$

2.2 Preparation of target protein

The 3D crystal structure of aspirin acetylated human cyclooxygenase-2 (PDB ID: 5F19) was obtained in pdb format from online protein data bank (PDB) database [27]. The structure was prepared by erasing water molecules, hetero atoms, inactive chain, and adding hydrogen atoms by PyMol (version 1.3) software package [28] and proceed with chain A of 5F19 for further energy minimization with Swiss-PDB Viewer software package (Version 4.1.0) [29].

2.3 Docking, analysis and visualization

Molecular docking is an important tool in computational drug design by which one can predict the predominant binding mode(s) of a ligand with a target protein. The optimized structures were subjected for molecular docking study against human prostaglandin synthase protein 5F19 (chain A) considering the protein as macromolecule and the drug as ligand where the parent drug considered as control and derivatives binding affinity was compared with the parent drug, which is a known Cox-2 inhibitor. The centre grid box was set at 64.84 Å, 73.29 Å and 57.94 Å in x, y and z direction respectively. Flexible docking was performed using PyRx software (version 0.8) [30]. Finally we saved both the proteins and ligand structures in pdb format for further non-bonding interactions and hydrogen bond surface calculation. Accelrys Discovery Studio (version 4.1) software was utilize to analyze and visualize the docking result [31].

2.4 Molecular dynamics (MD) simulations

To validate the predicted molecular docking, MD simulation was performed on the docked PCT-5F19 and D6-5F19 complex using AMBER14 force field implemented in YASARA dynamics program. Water molecules were added and system was neutralized by adding NaCl salt at 0.9% concentration. MD simulation was equilibrated for 100 ps followed by 20 ns production at 298 K with time step of 5 fs.

2.5 ADMET prediction

AdmetSAR online database was utilized to predict absorption, distribution, metabolism, excretion, and toxicity (ADMET) of PCT and its modified derivatives [32].

3 Results and discussion

3.1 Thermodynamic properties

Simple modification of drug can improve the physicochemical and binding properties [33]. With the modification of PCT, significant change observed in structural properties as well as energies, electrostatic potential, partial charge distribution, dipole moment. Free energy is an important criterion to represent the interaction of binding partners where both the sign and magnitudes carry the important information. Negative sign indicate spontaneous binding interactions, where as positive sign opposite to that [34]. Here, all calculated energies were found negative which indicate preferable binding interactions. The results (Table 1) demonstrate that all the derivatives of PCT are chemically stable and few of them are more stable than parent drug. By considering the structure of PCT and its derivatives D1–D9 in the gaseous phase to water solvent it has been observed that there is a gradual decrease of enthalpy, Gibbs free energy, thermal energy. The free energy difference of PCT and D4 are − 8.59 kcal/mol and − 8.84 kcal/mol respectively. Where D4 has the highest free energy (− 9.04 kcal/mol) due to the addition of higly electronegative fluorine atom. Moreover, increasing of negative values from PCT to D9, hence suggesting –F, and -SCH3 substituted PCT molecules (D3, D4, D6 and D9) energetically and configurationally preferable.

Table 1 Molecular formula, difference of electronic energy (∆E), enthalpy (∆H), Gibb’s free energy (∆G) in kcal/mol and dipole moment (∆μ, Debye) of paracetamol and its derivatives
Full size table

Improved dipole moment can enhance the polar nature of molecule, binding affinity and interactions with the amino acid residues of receptor protein [35]. Significant results are observed from gaseous phase to the water solvent. The change of dipole moment of PCT is 0.99 Debye where D5 shows the largest value (2.10 Debye) respectively and most of the modified drugs (D1, D3, D4, D5, D9) show higher value than PCT. Resulting, increased dipole moment causes an increasing of binding affinity against 5F19 protein (Table 2).

Table 2 Energy (eV) of HOMO, LUMO, Gap, hardness (η), and softness (S) of all drugs
Full size table

3.2 Frontier molecular orbital analysis

The HOMO–LUMO gap is related to the chemical hardness and softness of a molecule [36, 37]. Larger HOMO–LUMO gap related to high kinetic stability and low chemical reactivity where as small HOMO–LUMO gap is important for low chemical stability, addition of electrons to a high-lying LUMO and/or removal of electrons from a low-lying HOMO is energetically favourable in any potential reaction [38, 39]. In this analysis, it is found that, brominated (D7), iodinated (D8) and fluorinated (D6) derivatives show the lower HOMO–LUMO gap and higher softness, which may contribute to their chemical activity and polarizability than PCT.

3.3 Molecular electrostatic potential analysis

Molecular electrostatic potential (MEP) was calculated at B3LYP/6-31G+(d,p) level of theory to forecast the reactive sites for electrophilic and nucleophilic attack of all optimized structures [40]. Red colour represent maximum negative area which favourable site for electrophilic attack, blue colour indicate the maximum positive area which favourable site for nucleophilic attack and green colour represent zero potential area. MEP displays molecular size, shape as well as positive, negative and neutral electrostatic potential regions simultaneously in terms of colour grading. It is seen from MEP map, region having the negative potential are over electronegative atom (oxygen atoms) and having positive potential are over hydrogen atoms. Here, D3 have the maximum negative potentiality (− 0.2313 a.u., deep red) and D6 have the highest positive potentiality (+ 0.2769 a.u., deep blue) which suggesting the maximum possibility for the nucleophilic and electrophilic attack to the respected region (Figs. 1, 2).

Fig. 1
figure 1

Chemical structures of Paracetamol (PCT) and its designed derivatives

Full size image
Fig. 2
figure 2

Frontier molecular orbitals (HOMO, LUMO) and related energy of PCT and D6

Full size image

3.4 Atomic partial charge

Mulliken and NBO methods have been utilized to compute the atomic partial charges of all structures [41]. Dipole moment and polarizability influenced by the atomic partial charge [42]. In this study (Figure S5), all the hydrogen atoms show positive charge and other electronegative atoms (N,O) exhibit the negative charge in both methods. In all compounds, C-1 and C-10 have the negative charge but in D6, C-1 has the positive charge due to the attachment of highly electronegative fluorine atoms. Among the halogen substituents, fluorine and bromine exhibit negative charge (in D4, D6, D7) but chlorine and iodine exhibit the positive charge in D5, D9, and D8 respectively (Figs. 3, 4).

Fig. 3
figure 3

Molecular electrostatic potential map of PCT, D3 and D6

Full size image
Fig. 4
figure 4

A Docked conformation of PCT, D4, D6 and D8 at inhibition bounding site of receptor protein 5F19 (Chain A). B Superimposed view of them after flexible docking

Full size image

3.5 Binding affinity and non-bonding interactions

Binding affinities and non-bonding interactions are summarized in Table 3. The addition of halogen not only increase the physicochemical properties but also binding affinity and speciality. The incorporation of –CF3 group increased inhibiting and medicinal properties [43,44,45,46]. Carbon trifluoride group has a great significant applications in the field of agronomical dyes, pigments, pharmaceuticals, polymers and material science when it is incorporated to different organic molecules due to the strong electronegative and hydrophobic characteristics, which can be used in drug design to improve the selective functionality [47,48,49,50]. In our study, all the compounds show several significant hydrogen bonds. Some recent studies reported that, halogen bonding similar to hydrogen bonding plays crucial role for biological and chemical system [51, 52]. The binding affinity of D6, D7, D8, and D9 have considerably increased to − 7.5, − 7.1, − 7.1, − 7.0 and − 7.3 kcal/mol respectively from − 6.4 kcal/mol of PCT. Improved hydrogen bond is observed in D9 and D6 not only contribute in increasing binding affinity but also enhance the binding specificity [53]. All the molecules show multiple nonbonding interactions after the docking with 5F19. Two strong hydrogen bonds with His386 (2.36 Å) and Thr206 (2.84 Å), along with increasing number of halogen bonds and one hydrophobic interaction are found in D6-5F19 complex. Thus strong hydrogen bonding and halogen bonding have important role in increasing binding affinity of D6 with 5F19. The D7-5F19 complex is stabilize by two hydrogen bonds and one hydrophobic interaction (Table 3). The most important fact is that the -F atoms of –CF3 group in D6-5F19 complex interact with amino acid to form strong halogen bond interactions. We also observed six fluorine bond (halogen bond) in D6-5F19 and one in D4-5F19 complex (Fig. 5) which may have positive effect on protein-ligands stability, as well as in the binding affinity and selectivity. D6 has higher softness than PCT which may promote the polarizable nature of the drug that can stimulate greater nonbonding interactions with the receptor. Nonbonding interactions of (PCT, D6, D7, D8, and D9) with 5F19 involved the following amino acid residues Ala199, Thr206, His386, His388, Tyr385, Leu531, Leu534, Gly533 and Met522 (Table 3). Hydrogen bond surface having residues such as Ala202, Val349, Gin461, Gly135, Pro156, His39 and His386 help in creating strong donor regions where as residues such as Ala199, Thr206, Tyr385, Met522, Gly45, Gln461, Pro154, Phe529, Gly526, Asn382, His207, His388 and Tyr385 help in creating strong acceptor regions on the drug protein interaction surface (Fig. 6). This observation helped to confirm that PCT and its selected derivatives (D6, D8 and D9) are binding at the desired binding site of receptor protein after molecular docking. Ala202, Ala203, Ala528His388, Leu353 and Val524 are found to donate their π-electrons cloud toward the alkyl chain and the carbon attached.

Table 3 Binding energy and nonbonding interaction of Paracetamol (PCT) and its modified derivatives after flexible docking
Full size table
Fig. 5
figure 5

Nonbonding interactions of PCT, D6 and D8 with 5F19 (Chain A) generated by Discovery Studio

Full size image
Fig. 6
figure 6

Hydrogen bond surface of 5F19 (Chain A) with PCT, D6 and D8

Full size image

In case of D1, D2, D3, D4 and D5, X and Y positions of PCT are substituted by following functional groups –CH3, –OCH3, –SCH3, –F, –Cl and they their binding affinity with 5F19 are − 6.9, − 6.6, − 5.7, − 7.0, − 6.8 kcal/mol respectively. Three strong hydrogen bonds are found in D2 with Gly45 (2.00Å), Gln461 (2.04 Å), Gln461 (1.97 Å) and one hydrophobic interaction with Leu153. Where, D1, D3 and D4 also show strong hydrogen bond Tyr385 (2.25 Å), Gly135 (2.28 Å) and D5 and Ala199 (2.27 Å). D5 shows the maximum number of hydrophobic interactions with two amide-pi stacked (Ala527 and Gly526). All the compounds show non-conventional hydrogen bond except D1, D5 and D9. Generally, CH···O bonding term as non-conventional hydrogen bond, slightly weaker than its classical OH···O hydrogen bond [54, 55]. Those interactions may be the factors of improved binding affinity of D5 with 5F19 protein.

3.6 Molecular dynamics calculation

Parent compound (PCT) along with the best compound from docking analysis (D6) are analysed furthermore for conformational and binding stability through molecular dynamics study. 20 ns MD simulation was performed for PCT-5F19 and D6-5F19 complex. Root-mean-square deviation (RMSD) values of α-carbon and backbone atoms are investigated (Fig. 7). The RMSD values of PCT-5F19 and D6-5F19 complexes are plotted as a time-dependent function of the MD simulation. Very little RMSD deviation between alpha-carbon and backbone carbon atoms throughout the whole simulation process was evident and it indicates that the protein structure remains stable in both drug–protein complexes. The maximum RMSD value for the PCT-5F19 complex is 2.3 Å at around 16 ns whereas for the D6-5F19 complex the maximum RMSD value is same but at around 4.8 ns and 19.8 ns. Evaluation of molecular docking can also be performed through molecular dynamics simulations. In this study, both the drugs remain in the binding pocket after MD simulation similar to the docked pose of PCT-5F19 and D6-5F19 complexes (Fig. 8).

Fig. 7
figure 7

RMSD values of C-alpha and Backbone atoms of PCT and D6

Full size image
Fig. 8
figure 8

Superimposed structure of PCT-5F19 and D6-5F19 complex before and after 20 ns MD simulation (green and blue = before MD simulation; yellow and red = after MD simulation)

Full size image

3.7 ADMET analysis

Selected pharmacokinetics parameters are reported in Table 4. From ADMET results, it is found that all the structures show positive response for blood brain barrier (BBB) criteria, predicting that drugs can pass through the BBB. All the compounds are non-carcinogenic and show III category acute oral toxicity so, the modified drugs expected to be safe for topical biological use. All drugs are P-glycoprotein non-inhibitor and show positive response for human intestinal absorption. P-glycoprotein inhibition can block the absorption, permeability and retention of the drugs [56]. However, PCT and its modified derivatives show weak inhibitory property for human ether-a-go-go-related gene (hERG). Inhibition of hERG can lead to long QT syndrome [57], so more study of this aspect is necessary. The rat acute toxicity level of all drugs is higher than the PCT. So, modified derivatives have higher median lethal dose (LD50) than PCT.

Table 4 Selected pharmacokinetic parameters of paracetamol and its derivatives (probability values related to each of the parameters are given in the parenthesis)
Full size table

4 Conclusion

In this investigation, PCT and its modified derivatives are studied to explore their physicochemical properties and binding affinities with receptor protein. All the compounds are thermally stable and most of them have lower HOMO–LUMO gap with higher softness values. All the derivatives (except D3) have better binding affinity and interactions with the receptor protein compared to parent drug. All the derivatives have improved pharmacokinetic properties and safe for biological use. Moreover, D4, D7, D8, and D9 also show higher binding affinity and important non-bonding interactions than PCT. Where, D6-5F19 complex shows the highest binding affinity than others and remain inside the binding pocket of target protein after dynamics simulation. From the above results, it is observed that the insertion of halogen improved the physicochemical and binding properties. Considering the above investigation, D4, D6, D8, and D9 can be potential candidates for better performance.