In silico analysis of the action of saturated, monounsaturated, and polyunsaturated fatty acids against Echinococcus granulosus fatty-acid-binding protein 1

Computational biology methods

The FASTA sequence of the Echinococcus granulosus fatty acid-binding protein 1 (EgFABP1) protein from Echinococcus granulosus (ID: Q02970) was retrieved from the UniProt database (http://www.uniprot.org/),¹⁵ and subjected to the Blastp tool,¹⁶ followed by a sequence similarity search performed on E. granulosus (taxid:6210) within the non-redundant protein sequence database. To investigate the evolutionary relationship of EgFABP1 proteins, the retrieved sequences were uploaded in MegaX software and aligned using ClustalW.¹⁷ Utilizing the neighbor-joining technique utilized by iTOL, the resulting multiple sequence alignments was used to reconstruct and show a distance-based phylogenetic tree.¹⁸ Also, the EgFABP1-related sequences were uploaded into the Cytoscape platform,¹⁹ where the plugin StringApp¹⁹ was used to retrieve and extend the molecular network information from the STRING database²⁰ and the plugin cytoHubba was used to score and rank the nodes according to network properties, affording to the Maximal Clique Centrality (MCC) topological analysis method.²¹ The default parameters of Cytoscape were used to show the network; node size and color were manually changed according to the scores obtained from the MCC analysis. Finally, StringApp was also used to retrieve functional enrichment for Gene Ontology (GO) terms, regarding Biological Processes (BP), Molecular Functions (MF), and Cellular Components (CC), while the results were expressed in bubble plots employing the SRplot online platform (http://www.bioinformatics.com.cn/srplot).

Molecular geometry optimization

In this work, fourteen FAs were selected from PubChem database.²² The PubChem CID access code and chemical structure are described in Table 1. The geometries of these compounds were optimized using B3LYP/6-311++g(d,p) basis set level, implemented in ORCA software version 4.1.2.^23,24 The ground-state was verified by counting the imaginary frequencies of each compound.

Table 1. Fatty acids abbreviation, PubChem CID and 2D Representation.

FAs	Abbrev.	PubChem CID	2D Representation
Myristic acid	MRT	11005
Pentadecanoic acid	PDN	13849
Palmitic acid	PLM	985
Stearic acid	STA	5281
-Hexadecenoic acid	HDN	543268
Oleic acid	OLE	445639
Palmitoleic acid	PLA	445638
Trans-Vaccenic acid	TVA	5281127
Cis-Vaccenic acid	CVA	5282761
Arachidonic acid	ARA	444899
Docosahexaenoic acid	DHA	445580
Eicosapentaenoic acid	EPA	446284
Linoleic acid	LIN	5280450
Linolenic acid	LNA	5280934

The FAs optimized geometries were used for global and local reactivity descriptors calculation; the methodology is described below.

Global reactivity descriptors

Global reactivity descriptors are derived from conceptual Density Functional Theory (DFT),²⁵ this type of analysis allows us to understand the relationships among structure, stability and reactivity.²⁶ The global reactivity descriptors were calculated by the mean Koopmans’ theorem, which considers the HOMO (Highest Occupied Molecular Orbital) and LUMO (Lowest Unoccupied Molecular Orbital) energy orbitals (Table 2).

Local reactivity descriptors: Fukui function

Electrophilic, nucleophilic and radical attacks were calculated using Fukui functions.³⁵ The optimized geometries were selected compounds, and single point energy calculation was performed on the neutral, positive and negatively charged molecule using B3LYP/6-311++g(d,p) basis set implemented in ORCA software version 4.1.2.^23,24 The Fukui functions were calculated using Multiwfn³⁶ according to the following equations:

Nucleophilic attack:

Electrophilic attack:

Radical attack:

Where ${\rho }_N$ is the electron density at a point $r$ in space around the molecule. The $N$ corresponds to the number of electrons in the molecule. $N+1$ corresponds to an anion with an electron added to the LUMO of the neutral molecule. $N-1$ is the cation with an electron removed from the HOMO of the neutral.

Pharmacokinetic properties of FAs

The Molecular Weight (MW), Hydrogen Bonds (nHA, nHB), Polar Surface Area (PSA), Dissociation Constant (pKa), the logarithm of the molar solubility in water (logS), n-octanol/water partition coefficient (logP), n-octanol/buffer solution distribution coefficient at pH = 7.4 (logD7.4), and ADMET pharmacokinetic properties (absorption, distribution, metabolism, excretion, and toxicity) were performed on the ADMETlab 2.0 online platform (https://admetmesh.scbdd.com/). The prediction of pKa values was obtained by MolGpKa web server (https://xundrug.cn/molgpka).

Protein preparation

The protein structure of Echinococcus granulosus fatty-acid-binding protein I (EgFABP1) was downloaded from the Protein Data Bank^37,38 (https://www.rcsb.org) under access code PDB ID: 1O8V. This protein presents a hydroxylated cysteine (CSO63) that was considered in this work. We defined the AMBER99SB-ILDN force field parameters of the CSO63 and fourteen FAs in the ACPYPE³⁹ server online (https://www.bio2byte.be/login/?next=/acpype/). Finally, the water molecules present in the crystal structure were removed.

Docking determination and molecular dynamics simulations of complex

The coupling systems were solved in the DockThor server^40,41 (https://dockthor.lncc.br/v2) with a grid dimension of 5 × -7 × 10, the genetic algorithm settings with number of evaluations = 1000000, population size = 750, and number of runs = 24.

The molecular dynamics (MD) simulation was carried out using GROMACS v.2020.6⁴² with the AMBER99SB-ILDN force field. Fifteen MD simulations, EgFABP1 without ligand and EgFABP1 coupled to fourteen FAs was considered. The topologies were prepared considering a cubic box of 1 nm from the edge of the box to the protein surface. The TIP3P water model was used for the solvation of the different systems. 200000 energy minimization steps were settled with the steep-descent algorithm. Equilibrium dynamics simulations in the canonical ensemble at constant volume and the production simulation in the isobaric-isothermal at constant pressure (1 bar of reference pressure with Parrinello-Rahman barostat) were applied. The calculation time in isobaric-isothermal ensemble (NPT) was 200 ns. In both MD simulations, the V-rescale thermostat regulated the temperature at 309.65 K. The analysis of the MD was carried out with the Gromacs tools and plotted in the Gnuplot v. 4.5⁴³ software. Additionally, hydrophobic interactions between FAs and EgFABP1 were analyzed using the PLIP (Protein-Ligand Interaction Profiler) server.⁴⁴ The graphical visualizations were evaluated in the VMD (Visual Molecular Dynamics)⁴⁵ and UCSF Chimera v.1.15⁴⁶ software.

Binding free energy estimation by MM/GBSA analysis

The binding free energy from each coupling system (EgFABP1-FAs) was calculated by Molecular Mechanics-Generalized Born Surface Area (MM/GBSA) method. The mmpbsa.py⁴⁷ module from AmberTools 20⁴⁸ program and gmx MMPBSA v1.4.1⁴⁹ were used for these assays. The equations related to calculations of binding free energies ($\mathrm{\Delta }{G}_{\mathit{\text{bind}}}$) were the following:

Where $\mathrm{\Delta }{E}_{\mathit{MM}}$ is the variation between the minimized energy of the protein-ligand complexes (EgFABP1-FAs) that includes the van der Waals ($\mathrm{\Delta }{E}_{\mathit{vdw}}$) and electrostatic ($\mathrm{\Delta }{E}_{\mathit{ele}}$); $\mathrm{\Delta }{G}_{\mathit{GB}}$ is the electrostatic solvation energy, $\mathrm{\Delta }{G}_{\mathit{SA}}$ is the difference in surface area energies for protein and ligand and $-T\mathrm{\Delta }S$ refers to the contribution of entropy to temperature (T).

Correlation analysis

Spearman’s Rho (also called Spearman’s rank correlation coefficient) was used to estimate the correlation between protein stability and binding free energy of the systems (EgFABP1-FAs). The coefficient was calculated using Python SciPy package.⁵⁰ Spearman’s Rho is a rank-based non-parametric correlation measurement between two datasets, represented by a monotonic function. The Spearman’s Rho between feature X and Y is described by equation 7:

Where $\rho$ denotes the Spearman’s Rho and $n$ is the sample size.

Computational biology analysis

According to several studies, E. granulosus comprises a variety of strains with a comparatively high level of genetic variation. Also, numerous aspects of the parasite, such as its life cycle and transmission, pathogenicity, metabolic characteristics, and drug susceptibility, are impacted by this wide range of variety.⁵¹ Thirteen proteins with similarity to EgFABP1 were discovered by running a BLASTp search followed by phylogenetic analysis of the amino acid sequences against the E. granulosus proteome (Figure 1A). Also, understanding how successfully the proteins interfere with the complex regulatory machinery is essential, as developments in network biology have demonstrated that single-protein targeting is ineffective in treating complex disorders⁵²; the analysis of the extended protein network showed that the proteins with the highest centrality scores on the network are annotated as mediators of RNA polymerase II transcription (Figure 1B). Additionally, the analysis of the functional enrichment for Gene Ontology (GO) revealed that the terms linked to BP were retinol metabolic process (GO:0042572) and positive regulation of biological process (GO:0048518), while retinol and fatty acid-binding proteins have been associated to suppress host immunity during nematode infection.⁵³ Also, the GO terms associated with CC found in the analysis were intracellular organelle (GO:0043229) and organelle (GO:0043226), while the MF term was lipid binding (GO:0008289). EgFABP1 isoforms were found in the cytosolic, nuclear, mitochondrial, and microsomal fractions of the experiments, indicating that these molecules may be engaged in a variety of cellular activities⁵⁴ (Figure 1C).

Figure 1. Illustration of EgFABP1 isoforms discovered by running a BLASTp.

A. Phylogenetic analysis of the amino acid sequences against the E. granulosus proteome. B. Analysis of the extended protein network of E. granulosus. C. Analysis of the functional enrichment for Gene Ontology (GO).

Chemical reactivity of FAs

Global reactivity descriptors

DFT has been very successful in the last decade in predicting the chemical reactivity of various molecules.^55–58 Within the DFT theory, one of the most used to determine the stability and reactivity of a particular molecule are the frontier molecular orbitals HOMO and LUMO.^59,60 In Table 3, the energy values of the HOMO and LUMO orbitals for all the studied molecules match those reported in the literature for organic molecules.⁶¹ It is necessary to highlight that SFAs molecules (MTR, PDN, PLM and STA) showed the most negative values of the HOMO orbital, meaning that these molecules are the most likely to donate electrons in a chemical reaction.^60,61 According to the frontier orbital theory, the energy difference between the HOMO and LUMO orbitals (HOMO-LUMO gap) is a measure of the chemical stability of a molecule in a reaction.^61,62 High differences between the HOMO concerning the LUMO orbital indicate high stability of the molecule in a chemical reaction.⁶³ Analyzing the HOMO-LUMO gap energy behaviour shown in Table 3, the SFAs molecules obtained the greatest difference in this parameter compared to the other molecules, this means that SFAs compounds present greater chemical stability according to the frontier molecular orbitals theory and the global reactivity descriptors followed by the Koopmans theory,^27,56 which means that the HOMO and LUMO frontier orbitals were considered in the calculation.

Regarding the dipole moment, we observed that the value decreases depending on the unsaturation degree of the FAs due to their liposolubility.⁶⁴

To delve into the biological behaviour of the FAs under study, it is necessary to analyze the local reactivity descriptors such as the molecular electrostatic potential and the Fukui function. The molecular electrostatic potential provides information on the atoms of the chemical structure of a molecule with positive or negative charge density; and the Fukui function is a chemical observable that gives us the measure of which part of the molecule is more susceptible to a nucleophilic, electrophilic, or radical attack.

Local reactivity descriptors

As can be seen in Figure 2 ARA, DHA, and EPA molecules (PUFAs with the longest chains) had the most pronounced charge densities of all the compounds studied due to the structure of minimum energy in geometry optimization. These molecules have a greater torsion degree than the others, this implies that their unsaturations in addition to the carboxyl group, are very close. This generates a distribution of negative charges greater than the others. In general, the negative charge density (red lines) for all molecules were lower than the positive charge density (yellow lines), this is due to the FAs structure. These FAs molecules have between 14 and 22 carbon atoms, where the hydrocarbon chain predominates with a negative charge density towards the backbone. Likewise, the negative charge density is only located in the carboxyl group and in the unsaturations.

Figure 2. Contour lines of the electrostatic potential of the studied fatty acids.

The negative charge density in each atom is represented in red lines and the areas with positive charge density in each molecule are represented in yellow color.

The Fukui function is a local selectivity descriptor.^65–67 This parameter reveals the most susceptible molecule area of nucleophilic ($f^{+}$), electrophilic ($f^{-}$), and radicalary ($f^0$) attacks.⁶⁵ The Figure 3 shows that DHA, PDN, STA, and TVA molecules were the FAs with areas most susceptible to nucleophilic and electrophilic attack. These compounds have a significant concentration of positive and negative charges; the negative charge density is oriented around the unsaturations (with respect to DHA and TVA) and the carboxyl group.

Figure 3. Contour lines of the Fukui function of the studied fatty acids.

The negative isovalues in each atom are represented in purple lines and the areas with positive isovalues in each molecule are represented in light blue color.

The radicalary attack behavior was different from the $f^{+}$ and $f^{-}$. In this case, the density obtained for all molecules was low. The $f^0$ are focused fundamentally on the unsaturations from FAs. This result is due to fatty acids being susceptible to radical attack by reactive oxygen species in unsaturations.^68–72

The computed ADMET properties of FAs

The molecular weight (WM), number of hydrogen bond donors (nHBD), number of hydrogen bond acceptors (nHBA), and rotatable bonds (nRot) (see Table 4); have shown acceptable values according to the data based on the “Drug-Like Soft” rule (WM(100-600), nHBA(0-12), nHBD(0-7), and nRot(0-11)). The topological polar surface area (TPSA) values of all the compounds reported 37.30 A², which is an optimal value according to Verber’s rule ($\le$140 A²). Therefore, these results show that FAs would have good absorption and permeability. Likewise, all FAs are accepted for Lipinski’s rule of five ((MW$\le$500; logP$\le$5; hydrogen bond acceptor$\le$10; hydrogen bond donors$\le$5), indicating acceptable oral absorption and bioavailability.

The pharmacokinetic properties are reported in Table S1 from Supporting Information. The oral adsorption logP, distribution logD, and solubility logS descriptors were evaluated. The results showed that the logP values are above the desired range (0 to 3 log mol/L), indicating that these compounds are insoluble since they are trapped in the lipophilic bilayer. This could be confirmed by the resulting values of logD (1 to 3 log mol/L) that 197 are between 3 to 5, decreasing the solubility. Furthermore, all the compounds have a pKa = 4.8 except for 7-hexadecenoic acid with a pKa = 4.6. The fatty acids analyzed show optimal values for Caco-2 (>-5.15log cm/s); indicating that they have good permeability to intestinal cell membranes. They are also a substrate for P-glycoprotein (Pgp), where Pgp is a transporter efflux out of cells of xenobiotic molecules.

In general, the FAs have a high percentage (up to 90%) of plasma protein binding (PPB), contrary to DHA and EPA with ${\mathrm{PPB}}_{83.359\%}$ and ${\mathrm{PPB}}_{89.359\%}$, respectively. This union is determined by the structural properties of fatty acids that influence the oral bioavailability of the compounds, where a portion of the free compound could be absorbed through the membranes. The penetration of the blood-brain barrier (BBB) is important for those components that have activity in the central nervous system (CNS), some fatty acids have a better facility to penetrate the BBB than others; the values of lipophilicity, molecular weight, and TPSA allow us to predict a good penetration of the compounds to the BBB.

Cytochrome P450 (CYP) isoforms were also evaluated. The results indicated that fatty acids have a substrate/inhibitory effect on the CYP2C9 isoenzyme, which is important in the oxidation of xenobiotic compounds. ARA, DHA, and EPA fatty acids showed a moderate clearance of 6.011ml/min/kg, 9.988ml/min/kg, and 8.610ml/min/kg, respectively, while the other compounds had values below $\le$5 ml/min/kg, indicating low excretion. Among the results, DHA and EPA present relative toxicity due to their high lipophilicity.

Docking and molecular dynamics simulation

Best protein-ligand binding models were obtained using a genetic algorithm from the Dockthor server. Additionally, an affinity prediction scoring was obtained among the fourteen FAs. This prediction affinity is used to rank different ligands in virtual screening experiments considering the “Total Energy” of each FAs.^40,41 According to the score, DHA has a better affinity to EgFABP1 compared to the other FAs (see Table S2 from Supplementary Material).

Regarding Molecular Dynamics analysis, Figure 4 shows the plot of RMSD (root-mean-squared deviation), RMSF (root-mean-squared fluctuation), RG (radii of gyration), SASA (solvent-accessible surface area), and HBond (intramolecular hydrogen bonds from protein) during the simulation time of 200 ns in the isobaric-isothermal ensemble. The analysis of the RMSD plot showed us the MRT, ARA, and EPA with instability from 150 to 200 ns. Similar behaviour has the RMSF analysis, in which low and high fluctuations indicate the structural mobility of each residue in molecular dynamics. The RG plot indicates low compactness between EPA and ARA. While for the SASA study, we observed that the solvent access is similar in coupled systems and apparently, the number of HBonds is reduced in coupled systems.

Figure 4. Graphical representation of the MD simulation.

A. Root-mean squared deviation. B. Root-mean squared fluctuation. C. Radii of gyration. D. Solvent-Accessible Surface Area. E. Intramolecular hydrogen bonds from protein.

The Table 5 shows the average values for RMSD, RMSF, RG, SASA, and HBond for the last 20 ns (same time used to estimate binding free energy by MM/GBSA). However hydrophobic interactions and salt bridges were calculated from the last frame. One of the most outstanding systems (receptor-ligand) was EgFABP1-EPA because it presents an average value of RMSD (0.50±0.07 nm), RMSF (0.17±0.08 nm), RG (1.53±0.02 nm), SASA (91.37±2.26) higher and HBond (98.00±5.00) lower than the other systems. The PLA and VAC systems show a lower RMSD (0.14±0.01 nm) and RMSF (0.10±0.04 nm) than the other systems, indicating that the protein is more stable. Additional information about the hydrophobic interactions between EgFABP1-FA are reported in Table S3 from Supplementary Material.

Molecular Mechanics/Generalized Born Surface Area (MM/GBSA) estimation

The binding free energy estimation was analyzed by MM/GBSA method. Table 6 shows that DHA (-38.52±2.60 kcal/mol $\mathrm{\Delta }$G) and PLA (-37.07±4.55 kcal/mol $\mathrm{\Delta }$G) obtained the best energy followed by the VAC (-36.52±2.25 kcal/mol $\mathrm{\Delta }$G) and OLE (-35.58±3.64 kcal/mol $\mathrm{\Delta }$G). All these four ligands are UFAs (i.e DHA correspond to PUFAs and the others correspond to MUFAs). These results suggest that the EgFAPB1 has a preference for unsaturated fatty acids. In parallel, it was observed that the Van der Waals energies give the greatest energy contribution.

On the other hand, EPA obtained lower binding free energy than the other FAs. Figure 5, shows how DHA (better free energy) stabilizes the protein. However, EPA (lower free energy) destabilizes the protein. According to research by Undurti N. Das,⁷³ this effect could be beneficial in overcoming drug resistance against parasitic infections. Likewise, we observe that EPA has fewer hydrophobic interactions than DHA.

Figure 5. Interactions of DHA and EPA with EgFABP1.

Also, we evaluated Spearman’s Rho between RMSD and binding free energy ($\mathrm{\Delta }$ TOTAL). In the systems with better free energy, the protein tends to be more stable than without ligand, which could be inversely related. We obtain a Spearman’s Rho = 0.73 and p-value = 0.0029. The results indicate that the relationship between RMSD and binding free energy is a significant (p-value below the usual threshold of 0.05) correlation. In addition, according to most of the research reviewed,^74,75 the Spearman’s Rho obtained is classified as a strong correlation.