Quantum chemical calculations on calcium oxalate and dolichin A and their binding efficacy to lactoferrin: An in silico study using DFT, molecular docking, and molecular dynamics simulations

1.   Introduction

Urinary calculi, commonly known as kidney stones, are solid particles developed within the urinary system. Kidney stone disease frequently leads to severe pain, and it has the possibility of kidney failure in certain cases [1]. Currently, people have access to treatment procedures and preventive measures against urinary calculi; however, the prevalence of this disease is rapidly increasing in all developed and developing nations [2]. Urinary calculi have different chemical constituents, with the majority of them being calcium oxalate [3], [4]. Fourier transform infrared (FTIR) spectroscopy is a potential method for predicting the presence of calcium oxalate in these calculi [5]. The calcium oxalate molecule, consisting of one calcium and four electronegative oxygen atoms, has a possibility of hydrogen bond (H-bond) and nonbonded interactions with the proteins through the electric dipole and charge transfer process. Furthermore, these characteristics can be analyzed by computing various quantum chemical properties of compounds using the density functional theory (DFT) [6].

Figure 1.  (a) Crystallographic structure of lactoferrin liganded with iron oxalate in red color (1BKA.PDB) and (b) representation of active residues of the lactoferrin with ligand iron oxalate in red color and hydrogen bonds in yellow broken lines using molecular visualizing software PyMOL.

DownLoad: Full-Size Img PowerPoint

In the process of forming calcium oxalate type kidney stones, which encompasses nucleation to aggregation, numerous proteins serve as promoters; however, the precise mechanism of urolithiasis is still unknown [7], [8]. Lactoferrin, a member of the transferrin family, is found in the kidney and matrix of urinary calculi [8], [9]. Experimental investigations of the kinetic rate of calcium oxalate monohydrate growth in the presence of lactoferrin have revealed its role as a growth promoter through progression of layers on the surfaces of crystals [8]. Lactoferrin also shows the internal flexibility to bind oxalate ions without significantly changing the polypeptide conformations. The first bounded anions create the proper binding site for the metal ions and there exits the synergistic relationship between anions and metal ions for the tight binding [10]. Arginine121 (ARG121), Histidine253 (HIS253), Aspartic acid60 (ASP60), Tyrosine192 (TYR192), Serine191 (SER191), Arginine210 (ARG210), Phenylalanine208 (PHE208), Threonine117 (THR117), Glycine124 (GLY124), Alanine123 (ALA123), and Threonine122 (THR122) exhibit the active residues in the binding site of the lactoferrin (Figure 1). Therefore, calcium oxalate is one of the metal oxalates that can also bind with the active residues of lactoferrin and promotes the formation of calcium oxalate type kidney stones. The different phytocompounds play the inhibitory roles against the urinary proteins promoting urolithiasis [11]–[13]. Macrotyloma uniflorum (M. uniflorum) has been traditionally used for the treatment of kidney stone disease. Experimental studies identified that the legumes of M. uniflorum show the antilithiatic properties [14]–[16]. The potential bioactive compounds of M. uniflorum can inhibit the lactoferrin by the strong interactions with the residues in its oxalate binding site. The interactions between protein and ligand can be investigated by means of molecular docking and molecular dynamics (MD) simulations following molecular mechanics energies combined with the generalized Born and surface area continuum solvation (MM-GBSA) approach [17]–[19].

This research aims to investigate the quantum chemical properties of calcium oxalate and dolichin A (an effective phytochemical from M. uniflorum), along with their possible interactions within the active region of lactoferrin. The druglikeness of dolichin A and its conformational stability in complex with lactoferrin have been analyzed. The present study explores the understanding of the promoting role of lactoferrin and contributes to the development of a drug candidate against calcium oxalate induced kidney stone disease.

2.   Materials and methods

2.1.   FTIR characterization

Kidney stones were obtained from the Tribhuvan University Teaching Hospital following the receipt of approval from the Institutional Review Committee of the Institute of Medicine, Tribhuvan University, Kathmandu, Nepal. Surgically removed eleven samples were collected from 1 September 2022 to 1 March 2023, ensuring an unbiased selection process. Informed consent was obtained from all participants, including parents or legal guardians of minors, to allow for the collection and analysis of the kidney stones. The collected samples were meticulously cleaned with distilled water and subsequently air-dried at room temperature. Subsequently, the kidney stones were crushed using a stainless steel mortar and pestle. FTIR spectroscopy was employed for sample characterization in the spectral range of 400-4000 cm⁻¹ using a SHIMADZU IRAffinity-1S spectrophotometer.

2.2.   Gaussian calculations

The initial structures of calcium oxalate (Compound ID 33005) and dolichin A (Compound ID 44257432) were obtained from the PubChem database. These molecules were optimized using the DFT method at the Becke 3-parameter; Lee, Yang, Parr (B3LYP)/6-311++G(d,p) level of calculation in Gaussian 16W software [20]. From the optimized molecular structures, various molecular properties including the dipole moment, Mulliken charges, molecular electrostatic potential (MEP), highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO), global reactivity descriptors, and density of states (DOS) were computed in the gas phase. Additionally, infrared intensities in both the gas phase and water solvent using the integral equation formulation of the polarizable continuum model (IEFPCM) and thermodynamic properties in the gas phase of the calcium oxalate molecule were calculated using the DFT method at the same level of calculation. Furthermore, ultraviolet-visible (UV-Vis) spectra were also generated using the time-dependent (TD)-DFT method at the same level of calculation [20]. For visualization and in-depth analysis, GaussView 6 [21], GaussSum 3.0 [22], and VEDA4 software were utilized [23].

2.3.   Molecular docking

Crystallographic structure of lactoferrin (1BKA.PDB) was downloaded from https://www.rcsb.org [10] and refined using SWISS-MODEL [24], resulting in a template with global model quality estimate (GMQE) and distance constraints applied on model quality estimation (QMEANDisCo) global scores of 0.93 and 0.88 ± 0.05, respectively. Employing literature sources [25]–[28], twenty-four distinct chemicals present in M. uniflorum were primarily selected as ligands. These chemical compounds were downloaded in sdf format from the PubChem database [29]. Notably, among the selected compounds, dolichin A exhibited the best binding score within the active domain of lactoferrin as determined using PyRx software [30]. In the preparation of molecular docking using AutoDock4, removal of water molecules, addition of Kollmann's charges, integration of polar hydrogen atoms, and conversion of protein data bank (PDB) to protein data bank with partial charge and atom type (PDBQT) format for the proteins and ligand were performed using the graphical user interface of AutoDockTools [31], [32]. In this conversion process, partial atomic charges and atom types were added to the receptor proteins [33]. The atomic affinity potential of each atom of ligands with proteins was calculated using the AutoGrid process. Molecular docking was performed using grid points (34 × 30 × 30) and (50 × 50 × 50) for calcium oxalate and dolichin A along the x, y, and z axes, respectively. Out of several search methods in molecular docking, the Lamarckian genetic algorithm was applied for the protein and ligand [34]. A large number of conformations were developed and the best solution was obtained after analyzing the different predicted energies of the conformations [34]. To ensure exhaustiveness in the conformational search, a maximum of 2,500,000 energy evaluations were performed, and a number of individuals in the population was set to 300. Additionally, the maximum number of generations was set to 27,000, and the number of docking procedures was extended to 100 [32], [34]–[36]. For validation purposes, the DockRMSD online tool was utilized to calculate the root mean square deviation (RMSD) between the generated coordinates of the docked ligand and the original coordinates of the native one [37]. The protein-ligand interactions were analyzed using PyMOL 2.5.2 [38], LigPlot+ v.2.2 [39], and BIOVIA Discovery Studio Visualizer v21.1.0.20298 [40].

2.4.   Molecular dynamics simulations

The protein-ligand complex having the best binding scores from molecular docking was subjected to MD simulations using GROMACS 2019.6 [41]. The ligand dolichin A was parameterized using the SwissParam online tool [42]. The simulations employed the chemistry at Harvard macromolecular mechanics (CHARMM) force field, using the CHARMM27 parameter set [43] and the complexes were solvated using the transferable intermolecular potential with a 3 points (TIP3P) water model [44]. Neutralization was achieved by adding Na⁺ and Cl⁻ ions followed by a 50,000-step minimization using the steepest descent method. Prior to the production run, the system was equilibrated in the constant-temperature, constant-volume (NVT) ensemble with temperature control using the V-rescale (modified Berendsen thermostat). Again, a subsequent equilibration was performed in the constant-temperature, constant-pressure (NPT) ensemble, employing the V-rescale thermostat and isotropic Berendsen pressure coupling. Finally, production runs for MD simulation were performed using the V-rescale thermostat and Parrinello-Rahman barostat, maintaining a constant temperature of 300 K and pressure of 1 bar for a duration of 100 ns with a time step of 1 fs. Trajectories were recorded every 0.1 ns, and conformational changes were analyzed through RMSD, root-mean-square-fluctuation (RMSF), radius of gyration (R_g), solvent accessible surface area (SASA), and H-bond interactions [19].

2.5.   Druglikeness and toxicity prediction

The evaluation of physicochemical, lipophilic, pharmacokinetic, and druglikeness parameters and toxicity of dolichin A, based on its canonical simplified molecular input line entry system in PubChem, was performed using the SwissADME tool, admetSAR 2.0, and ProTox-II [45]–[47].

2.6.   MM-GBSA calculations

The open source packages for Groningen machine for chemical simulations (GROMACS), gmx_MMPBSA and gmx_MMPBSA_ana were used to analyze the binding free energy contributors in lactoferrin-dolichin A conformers during MD simulations [48]. The binding free energy of the protein-ligand complex (ΔG_bind) is calculated using Equation 2.1 [18]

where, G_PL is the bound free energy of the protein-ligand complex, whereas G_P and G_L are the unbound free energies of protein and ligand, respectively.

3.   Results and discussion

3.1.   Geometry optimization

Calcium oxalate and dolichin A were optimized by the DFT method at the B3LYP/6-311++G(d,p) level of calculation. The optimized structures have been shown in Supplementary Figure 1S. The calculated dipole moments of calcium oxalate and dolichin A are found 14.880 Debye and 4.264 Debye, respectively. Dipole-dipole interactions between protein and ligand can affect the binding energy of the system [49]. So, these compounds can also interact through dipole-dipole interactions with proteins related to nephrolithiasis.

3.2.   Frontier molecular orbitals

Electronic absorption spectra were obtained using the TD-DFT method at the B3LYP/6-311++G(d,p) level of calculation (Figure 2a,b). The wavelengths of 335 nm, 406 nm, and 507 nm are responsible for the electronic transition of calcium oxalate. Among these transitions, the maximum oscillatory strength is observed for wavelength 507 nm due to major contribution from the electron donor (HOMO) to electron acceptor (LUMO) transition of electrons (Table 1). Again, wavelengths of 268 nm, 265 nm, and 278 nm are observed for dolichin A. The wavelength of 271 nm is due to the transition of HOMO to LUMO+1, and HOMO to LUMO+2 shows maximum oscillatory strength (Table 1).

Figure 2.  Calculated electronic absorption spectra of (a) calcium oxalate and (b) dolichin A.

DownLoad: Full-Size Img PowerPoint

Figure 3.  The frontier molecular orbitals significantly contributing to the electronic transitions in (a) calcium oxalate and (b) dolichin A.

DownLoad: Full-Size Img PowerPoint

Figure 4.  DOS in both occupied and virtual orbitals in (a) calcium oxalate and (b) dolichin A.

DownLoad: Full-Size Img PowerPoint

In the frontier molecular orbitals (Figure 3), the green and red regions surrounding the atoms correspond to negative and positive charge phases, respectively [50]. In the case of calcium oxalate, a significant portion of the charge density is predominantly localized around the carbon and oxygen atoms, as shown in Figure 3a. In dolichin A, there is a noticeable shift in charge density from one molecular orbital region to another, as visualized in Figure 3b. This visual representation effectively conveys the nature of charge distribution within the molecule [51]. The energy gaps, represented as the difference between E_HOMO and E_LUMO, are found to be 3.07 eV and 4.87 eV for calcium oxalate and dolichin A, respectively (Figure 3). The number of states per unit energy interval at a given energy level in both occupied and virtual orbitals for calcium oxalate and dolichin A are illustrated with the DOS spectrum in Figure 4. The HOMO-LUMO orbitals (Figure 3) and the spectra of DOS (Figure 4) show good agreement with each other. These results indicate that both calcium oxalate and dolichin A can facilitate charge transfer process and have enhanced potential for biological activities [52].

Table 1.  Calculated electronic properties of calcium oxalate and dolichin A.

Compounds	Maxiumum absoption wavelength (nm)	Oscillatory strength	Major contributions
Calcium oxalate	507	0.014	HOMO → LUMO (99%)
	406	0.007	HOMO-1 → LUMO (99%)
	335	0.000	HOMO-2 → LUMO (99%)
Dolichin A	289	0.008	HOMO→LUMO (97%)
	271	0.049	HOMO→LUMO+1 (64%), HOMO→LUMO+2 (10%)
	268	0.009	HOMO→LUMO+2 (64%), HOMO→LUMO+3 (27%)

---

 | Show Table

DownLoad: CSV

3.3.   Global reactivity descriptors

The ionization potential (I = −E_HOMO) and electron affinity (A = −E_LUMO) of calcium oxalate were calculated using Koopmans' theorem [53]. The global reactive descriptors of the chemical can be predicted by the chemical hardness (η = (I − A)/2), chemical softness (S = 1/η), electronic chemical potential (µ = −(I + A)/2), and global electrophilicity index (ω = µ²/(2η)) [53]–[55]. Calcium oxalate and dolichin A show positive values for chemical hardness and electrophilicity index (Table 2). This observation suggests that the compounds are suitable for charge transfer processes and can impact the binding energy in protein-ligand interactions [56]. Furthermore, the chemical softness of calcium oxalate (0.65 per eV) is found to be higher in comparison with dolichin A (0.41 per eV), which suggests calcium oxalate is more reactive than dolichin A [52].

Table 2.  Calculated global reactivity descriptors of calcium oxalate and dolichin A.

Compounds	Ionization potential (I) (eV)	Electron affinity (A) (eV)	Chemical hardness (η) (eV)	Chemical softness (S) (eV)−1	Electronic chemical potential (µ) (eV)	Global electrophilicity index (ω) (eV)
Calcium oxalate	6.22	3.15	1.54	0.65	−4.68	7.15
Dolichin A	5.56	0.68	2.44	0.41	−3.12	1.99

---

 | Show Table

DownLoad: CSV

3.4.   Molecular electrostatic potential and Mulliken atomic charges

MEP map and Mulliken atomic charge are helpful for the interpretation of reactive behavior of the molecules related to electrophilic and nucleophilic attacks [57]. The visual representation of MEP shows the relative polarities within the calcium oxalate and dolichin A molecules (Figure 5). In the MEP map, an increase in electrostatic potential regions corresponds to the order red < yellow < green < blue, with the color codes ranging from −0.273 a.u. to 0.273 a.u. for calcium oxalate and from −0.121 a.u. to 0.121 a.u. for dolichin A, respectively. Within the calcium oxalate molecule, relative positive potential region is observed in the vicinity of calcium atom, while relative negative potential regions are concentrated around oxygen atoms. In dolichin A, relative positive potential regions are found around hydrogen atoms, while relative negative potential regions are found around oxygen atoms. Both positive and negative potential regions are also present around the carbon atoms. The distribution of Mulliken atomic charges of calcium oxalate and dolichin A molecules are shown in Figure 6. In calcium oxalate, the Ca7 atom shows the maximum positive Mulliken atomic charge. Ca7 is bonded to O1 and O4 oxygen atoms, resulting in the induction of more negative charges on the oxygen atoms O1 and O4 than that of O2 and O3 (Table 1S). In dolichin A, all hydrogen atoms exhibit positive Mulliken atomic charges, while all oxygen atoms show negative Mulliken atomic charges (Table 1S). Carbon atoms in calcium oxalate show positive Mulliken atomic charges, whereas in dolichin A, carbon atoms exhibit both positive and negative Mulliken atomic charges. From the MEP and Mulliken atomic charge distributions, calcium oxalate and dolichin A molecules show the polar properties, and these electrophilic and nucleophilic sites interact with the lactoferrin (Figure 2S).

Figure 5.  MEP map with contour lines (a) calcium oxalate and (b) dolichin A.

DownLoad: Full-Size Img PowerPoint

Figure 6.  Mulliken atomic charge distributions of calcium oxalate and dolichin A.

DownLoad: Full-Size Img PowerPoint

3.5.   Vibrational analysis

Calcium oxalate has seven atoms with fifteen modes (3N-6) of vibration [58]. The calculated vibrational frequencies (scaled and unscaled) compared with the experimental results in the gas phase and the water solvent along with potential energy distribution (PED) are shown in Tables 3 and 4. Due to the presence of a calcium atom, the bond strength between (O1-C5) and (O4-C6) becomes weaker than (O2-C6) and (O3-C5). So, the frequencies of (O1-C5) and (O4-C6) have been downshifted by 533 cm⁻¹ in symmetric and by 651 cm⁻¹ in asymmetric stretching, with respect to corresponding modes of (O2-C6) and (O3-C5) in the gas phase (Table 3). Furthermore, frequencies of (O1-C5) and (O4-C6) have been downshifted by 255 cm⁻¹ in symmetric stretching and by 347 cm⁻¹ in asymmetric stretching, with respect to corresponding modes of (O2-C6) and (O3-C5) in the water solvent (Table 4). The calculated frequencies are uniformly scaled by using the wavenumber linear scaling procedure as: ν_observed = ν_calculated × (1.0087 - 0.0000163 ν_calculated) for the DFT method at the B3LYP/6-311++G(d,p) level of calculation [59].

Table 3.  Calculated vibrational analysis of calcium oxalate along with potential energy distribution in the gas phase, available experimental FTIR results of calcium oxalate monohydrate, and FTIR results of kidney stones.

Calc. FTIR intensity	Calc. Freq. (unscaled) cm−1	Calc. Freq. (scaled) cm−1	Exp. FTIR freq. cm−1 [5]	Exp. FTIR freq. cm−1	Assignment [PED (> 10%)]
1901.32	1800	1763			ν[O2-C6 & O3-C5](93)
261.07	1782	1746	1618	1609	νa[O2-C6 & O3-C5](91)
1498.56	1245	1230	1317	1383	ν[O1-C5 & O4-C6](75)
374.46	1105	1095		1314	νa[O1-C5 & O4-C6](76), δa[(O2-C6-O4) & (O1-C5-O3, C5-O1-Ca7 & C6-C5-O1)](14)
348.09	856	852			ν[O1-C5 & O4-C6](11), νa[C5-C6 & Ca7-O1](17), δ[O2-C6-O4 & O1-C5-O3](62)
0.00	833	829			γ[O2-C5-O4-C6 & O3-C6-O1-C5](97)
522.08	784	781			νa[O1-C5 & O4-C6](11), δa[(O2-C6-O4) & (O1-C5-O3, C5-O1-Ca7 & C6-C5-O1)](72)
435.81	554	553			ν[Ca7-O1](80)
30.10	543	543			δa[C6-C5-O1 & C5-C6-O4](79), νa[O1-C5 & O4-C6](11)
433.46	493	493			γa[O2-C5-O4-C6 & O3-C6-O1-C5](90)
416.55	370	371			νa[C5-C6 & Ca7-O1](24), δ[C5-O1-Ca7, C6-C5-O1, C5-C6-O4 & O2-C6-O4](50)
472.06	364	365			νa[C5-C6 & Ca7-O1](36), δ[O2-C6-O4 & O1-C5-O3](17), δ[C5-O1-Ca7, C6-C5-O1, C5-C6-O4 & O2-C6-O4](30)
46.43	321	322			δa[C5-O1-Ca7 & C5-C6-O4](85)
779.45	118	119			τ[Ca7-O1-C5-C6](90)
0.00	52	52			τ[O4-C6-C5-O1](97)

---

Symbols: ν, δ, τ and γ represent the stretching, bending, torsion and out-of-plane modes of symmetric vibrations, respectively, and the subscript "a" represents the corresponding asymmetric mode for each mode of vibration. Abbreviations: Calc., Exp. and Freq. represent calculated, experimental and frequency, respectively.

 | Show Table

DownLoad: CSV

Table 4.  Calculated vibrational analysis of calcium oxalate along with potential energy distribution in the water solvent, available experimental FTIR results of calcium oxalate monohydrate, and FTIR results of kidney stones.

Calc. FTIR intensity	Calc. Freq. (unscaled) cm−1	Calc. Freq. (scaled) cm−1	Exp. FTIR freq. cm−1 [5]	Exp. FTIR freq. cm−1	Assignment [PED (> 10%)]
4719.19	1645	1615			ν[O2-C6 & O3-C5](83), ν[O1-C5 & O4-C6](12)
253.03	1637	1608	1618	1609	νa[(O4-C6 & O1-C5) & O3-C5](86)
1540.01	1379	1360	1317	1383	ν[O1-C5 & O4-C6](72), ν[C5-C6](13)
1044.70	1277	1261		1314	νa[O4-C6 &(O1-C5, O3-C5)] (89)
172.24	871	866			δ[O2-C6-O4 & O1-C5-O3] (49)
					ν[O1-C5 & O4-C6](12), ν[C5-C6](17)
0.22	859	855			γ[O2-C5-O4-C6 & O3-C6-O1-C5](98)
904.65	772	769			δ[ O1-C5-O3](82)
128.24	573	573			νa[(O4-C6 & O1-C5) & O3-C5](11), γa[O2-C5-O4-C6](21), δ[C6-C5-O1](60)
556.00	481	481			γ[O2-C5-O4-C6](80), δ[C6-C5-O1](12)
372.20	476	477			δ[O1-C5-O3 & O2-C6-O4] (30), ν[C5-C6](59)
417.29	344	345			δ[C5-C6-O4, C6-C5-O1, O2-C5-O4 & O1-C5-O3](82)
2679.31	275	276			νa[C5-C6 & Ca7-O1](44), δ[C5-O1-Ca7](36)
911.14	256	257			νa[C5-C6 & Ca7-O1](41), δ[C5-O1-Ca7](44)
1318.69	102	103			τ[C6-C5-O1-Ca7](95)
6.68	82	83			τ[O1-C5-C6-O4](97)

---

 | Show Table

DownLoad: CSV

In the gas phase, symmetric stretching of (O2-C6) and (O3-C5) and that of (O1-C5) and (O4-C6) with a minor contribution from the bending modes of (O2-C6-O4 and O1-C5-O3) (Table 3), are responsible for the most intense FTIR intensity at frequencies 1763 cm⁻¹, and 1230 cm⁻¹, respectively (Figure 3S). In the water solvent, the most intense FTIR spectra are observed at frequencies of 1615 cm⁻¹, 1360 cm⁻¹ and 1261 cm⁻¹ (Figure 7c). The frequency of 1615 cm⁻¹ is due to the symmetric stretching of all (O-C) bonds. The frequency of 1360 cm⁻¹ is due to the symmetric stretching of (O1-C5) and (O4-C6) with a minor contribution by stretching of (C5-C6). The frequency of 1261 cm⁻¹ is contributed by the asymmetric stretching of (O1-C5) with (O2-C6) and (O4-C6) (Table 4).

3.6.   Experimental FTIR analysis of kidney stones

Out of eleven samples, FTIR spectra of eight samples having similar patterns have been considered for further analysis (Figure 7). Average values with standard deviations of major peaks are observed at 1609 ± 1 cm⁻¹, 1383 ± 1 cm⁻¹, and 1314 ± 1 cm⁻¹. FTIR spectra of calcium oxalate in the gas phase and in the tetrahydrofuran (THF), carbon tetrachloride (CCl4), and water solvents were calculated at the B3LYP/6-311++G(d,p) level of calculation (Figure 4S). The calculated FTIR spectrum in the water solvent nearly matches with the experimental FTIR peaks (Figure 7) and also with the peaks reported for calcium oxalate monohydrate [5].

3.7.   Thermodynamic properties

The thermodynamic properties of compounds are intricately linked with their interactions with the target proteins [60]. The variation of thermal energy correction (E), molar heat capacity at constant volume (C_V), and entropy (S) of calcium oxalate in the gas phase with temperatures ranging from 50 K to 500 K with an increment of 50 K are presented in Table 2S. The correlation equations of thermal energy correction, molar heat capacity at constant volume, and entropy are found as:

Thermal energy (E) = 16.04170 + 0.00567 T + 0.00003 T²

Heat capacity (C_V) = 5.66920 + 0.05303 T - 0.00001 T²

Entropy (S) = 48.52010 + 0.16623 T - 0.00018 T²

Based on the correlation graph best fitted with the coefficient of determining of nearly one, each of the three calculated thermodynamic parameters shows an increasing trend with temperature (Figure 8). This behavior is due to the rise in the vibrational intensities of molecules with increasing temperature [61].

3.8.   Molecular docking

The protein-ligand complexes of lectoferrin with calcium oxalate and dolichin A were generated using molecular docking. For the validity of the molecular docking, the native ligand FeC₂O₄ was extracted from the lactoferrin receptor (1BKA.PDB) and redocked in the vicinity of its original position using AutoDock4.2 (Figure 5S). The RMSD of the generated coordinates of docked ligand from its original coordinates is found at 0.208 Ǻ, and such validated setup was used in molecular docking of calcium oxalate and dolichin A into the lactoferrin.

The molecular docking of calcium oxalate in the vicinity of the metal-oxalate binding regions of the lactoferrin results the binding energy of −3.86 kcal/mol (Table 5). Specifically, residues ARG121 and ASP60 form H-bonds, and the residues PHE63, PHE183, Leucine59 (LEU59), and Proline42 (PRO42) engage in hydrophobic interactions with calcium oxalate (Table 6).

Figure 7.  (a) Experimental FTIR spectra of kidney stone samples, (b) FTIR spectrum of the average intensity of 8 kidney stone samples, and (c) calculated FTIR spectrum of calcium oxalate in the water solvent.

DownLoad: Full-Size Img PowerPoint

Figure 8.  Variation of thermal energy correction (E), molar heat capacity at constant volume (C_V), and entropy (S) with temperature of calcium oxalate in the gas phase.

DownLoad: Full-Size Img PowerPoint

Figure 9.  Three-dimensional representation of the binding domain of lactoferrin interacting with (a) calcium oxalate and (b) dolichin A using PyMOL. The ligands and interacting residues are shown in yellow and red colors, respectively.

DownLoad: Full-Size Img PowerPoint

Table 5.  Molecular docking results of calcium oxalate and dolichin A with lactoferrin at temperature 298.15 K.

Compounds	Final intermolecular energy [vdW + H-bond + desolv+electrostatic energy] (kcal/mol)	Final total internal energy (kcal/mol)	Torsional free Energy (kcal/mol)	Unbound system's energy (kcal/mol)	Estimated free energy of binding (kcal/mol)
Calcium oxalate	−3.86	0.00	0.00	0.00	−3.86
Dolichin A	−9.57	−1.69	+1.79	−1.69	−7.78

---

 | Show Table

DownLoad: CSV

Figure 10.  Two-dimensional visualization of the interactions of lactoferrin with (a) calcium oxalate and (b) dolichin A using BIOVIA Discovery Studio.

DownLoad: Full-Size Img PowerPoint

Table 6.  Active residues of lactoferrin and their interactions with calcium oxalate and dolichin A in the best docked pose after molecular docking.

	Active amino acids	Types of interaction	Interacting atoms (Receptor − ligand)	Bond length (Å)
Calcium oxalate	ARG121	H-bond	NH1 − O4	3.02
	ASP60	H-bond	N − O2	3.04
	PHE63, PHE183, LEU59, PRO42	Non-bonded
Dolichin A	HIS253	H-bond	ND1 − O3	3.08
	HIS253	H-bond	ND1 − O4	2.95
	GLU80	H-bond	OE2 − O4	3.03
	THR58, ASP297, ARG121, PHE300, LYS301,GLU15, VAL11, VAL57, LEU59, PRO42, PHE190, SER191, TYR192, LEU299	Non-bonded

---

 | Show Table

DownLoad: CSV

The binding energies of twenty-four chemicals found in M. uniflorum [25]–[28] and two-dimensional docking images at the oxalate binding site of lactoferrin are presented in Table 3S and Figures 6S and 7S. Among all compounds, dolichin A exhibits the strongest binding affinity with the lactoferrin resulting in the binding energy of −7.78 kcal/mol (Table 5). The residues HIS253 and Glutamic acid80 (GLU80) form the H-bonds whereas the residues THR58, ASP297, ARG121, PHE300, LYS301, GLU15, ASP60, Valine11 (VAL11), VAL57, LEU59, PRO42, PHE190, SER191, TYR192, and LEU299 perform hydrophobic interactions with dolichin A (Table 6). Three-dimensional and two-dimensional visualizations of calcium oxalate and dolichin A in complex with lactoferrin are shown in Figures 9, 10, and 2S. The results of molecular docking show that dolichin A is a potential inhibitor of the lactoferrin.

3.9.   Molecular dynamics simulations

Figure 11.  Distributions of RMSD (a, b), R_g (c, d), and SASA (e, f), along with their relative frequencies for both the apo form and dolichin A binding complex of lactoferrin during 100 ns MD simulations.

DownLoad: Full-Size Img PowerPoint

Stability of a protein-ligand system can be examined by the trajectory analysis during its MD simulations. The time evolution changes of RMSD and its relative frequency in the cases of lactoferrin in apo form as well as liganded with dolichin A have been shown in Figures 11a,b. During a 100 ns trajectory of the apo form, the average RMSD of the protein backbone is 1.9 ± 0.4 Å with the maximum value of 4.0 Å. The most frequently observed RMSD in the apo form is 2.0 Å. The average RMSD of the protein backbone in dolichin A bound lactoferrin is 1.9 ± 0.3 Å with a maximum value of 3.0 Å, and the most prevalent RMSD is 1.8 Å. Again, the average RMSD of heavy atoms of the ligand with respect to the protein backbone is found to be 2.1 ± 0.3 Å. These outcomes show the stability of both apo and ligand binding forms. Furthermore, dolichin A fits well within the binding pocket of lactoferrin, and there is slightly less deviating conformations of the dolichin A bound protein system in comparison with the apo structure.

The compactness of a protein system is evaluated by calculating R_g of its backbone atoms. Throughout the simulation, both the apo and ligand binding forms consistently exhibit the R_g values of 29.7 ± 0.2 Å and 29.8 Å ± 0.2 Å with the maximum values of 30.5 Å and 30.4 Å, respectively (Figures 11c,d). These results show the compactness of the lactoferrin is nearly consistent after binding with dolichin A.

The impact of solvation on biomolecules is assessed through the calculation of SASA. The SASA of both apo and dolichin A binding forms remain relatively stable, having the values 305 ± 4 nm² and 307 ± 3 nm², respectively (Figures 11e,f). The maximum SASA values are 316 nm² for the apo form and 317 nm² for the ligand binding form. There is no significant change in solvent accessibility of lactoferrin after binding with dolichin A.

Figure 12.  (a) RMSF of the amino acid residues in apo form and dolichin A bound lactoferrin where the metal-oxalate binding active residues are marked with the dotted lines and (b) distribution of H-bonds between lactoferrin and dolichin A during MD simulations.

DownLoad: Full-Size Img PowerPoint

The flexibility of residues in the protein systems is obtained through RMSF calculation. For the apo and ligand binding forms, the average RMSF of the backbone is observed to be 1.3 ± 0.5 Å with maximum values of 5.4 Å and 5.8 Å, respectively (Figure 12a). Among the active residues, LEU59 and ASP60 in apo form, whereas ARG121 and THR122 in ligand binding form, show the lowest RMSF of 0.56 Å (Table 4S). From the 100 ns trajectory analysis, one can conclude that the lactoferrin appears to be stable throughout the MD simulations.

The distribution of H-bonds between dolichin A and active amino acids of lactoferrin during the 100 ns of MD simulations is illustrated in Figure 12b. In the initial 60 ns, a single H-bond has been observed over the majority of time. During the time interval of 60 to 100 ns, two H-bonds are predominantly observed. The maximum of five H-bonds have been formed between dolichin A and lactoferrin, which is observed at 0.9 ns of equilibration run.

3.10.   MM-GBSA calculations

The post MD free energy of dolichin A in complex with lactoferrin has been analyzed by MM-GBSA calculations. The major contributors to the total binding free energy are van der Waals energy (vdW), electrostatic energy (EEL), polar solvation energy (EGB), and nonpolar solvation energy (ESURF). The total binding free energy contributors of dolichin A in complex with lactoferrin taking all 1000 frames under consideration are illustrated in Figure 13a. The total MM-GBSA free energy written as average ± SD is −17.61 ± 4.03 kcal/mol, where the major contributors are vdW energy (−37.13 ± 2.44 kcal/mol) and electrostatic energy (−23.41 ± 5.56 kcal/mol) are depicted in Table 7. The binding free energies contributed by active residues of lactoferrin and dolichin A are illustrated in Figure 13b. The ligand, dolichin A contributes −10.10 ± 2.57 kcal/mol to the total MM-GBSA free energy. Among all residues, ARG121 and HIS253 show the lowest binding free energies of −1.58 ± 0.76 kcal/mol and −1.14 ± 0.43 kcal/mol, respectively (Table 8). Additionally, Figures 13c-d depict variations in vdW energies and total binding free energies across different trajectories in 100 ns MD simulations. During the simulations, the 33^rd frame exhibits the minimum vdW energy of −44.25 kcal/mol, whereas the 862^nd frame shows the lowest total binding energy of −29.67 kcal/mol.

Table 7.  Average binding free energies of lactoferrin-dolichin A complex using the MM-GBSA method (all values are in kcal/mol).

Van der Waals energy	Electrostatic energy	Polar solvation energy	Non-polar solvation energy	Total binding free energy
−37.13 ± 2.44	−23.41 ± 5.56	48.05 ± 3.48	−5.12 ± 0.20	−17.61 ± 4.03

---

 | Show Table

DownLoad: CSV

Table 8.  The binding free energies contributed by active amino acid residues of lactoferrin and dolichin A (LIG).

Residues	Average energy ± SD (kcal/mol)
VAL:11	−0.89±0.29
GLU:15	1.88±1.01
PRO:42	−1.06±0.30
THR:58	−0.51±0.30
LEU:59	−0.72±0.19
ASP:60	1.15±0.56
PHE:63	−0.79±0.21
GLU:80	0.71±0.26
ARG:121	−1.58±0.76
PHE:190	−0.60±0.35
TYR:192	−0.73±0.43
HIS:253	−1.14±0.43
LIG:692	−10.10±2.57

---

 | Show Table

DownLoad: CSV

Figure 13.  MM-GBSA results of dolichin A in complex with lactoferrin during 100 ns MD simulations: (a) binding free energy contribution by different interactions, (b) binding free energy contributions by active residues and ligand, (c) vdW interaction energy contribution in different frames with moving average, and (d) total binding free energy in different frames with moving average, where vertical lines in bars represent the standard deviations of the corresponding values from 1000 frames.

DownLoad: Full-Size Img PowerPoint

3.11.   Druglikeness and toxicity

Incorporating in silico investigation into pharmacokinetics and toxicity assessments of potential drug candidates is helpful for optimizing cost-efficiency, time management, and resource utilization within the realm of drug development [62]. The physicochemical and pharmacokinetic properties, lipophilicity, druglikeness, and toxicity of dolichin A are detailed in Tables 5S - 10S. Dolichin A possesses a molecular weight of 340.37 g/mol (< 500 g/mol), number of H-bond acceptors 5 (< 10), and number of H-bond donors 3 (< 5) (Table 5S). Additionally, dolichin A exhibits lipophilicity (MLogP) value of 1.90 (< 4.15). These findings indicate that dolichin A adheres to Lipinski's rule of five, underscoring its druglike properties [63]. The bioavailability radar (Figure 8S) comprehensively evaluates six critical physicochemical attributes, namely, lipophilicity, size, polarity, solubility, flexibility, and saturation (Tables 5S and 6S), where dolichin A falls entirely within the pink region. This positioning implies a promising potential for dolichin A as a candidate molecule suitable for the drug development [45]. In the BOILED-Egg model, considering lipophilicity (WlogP) and total polar surface area (TPSA), dolichin A is situated within the white region with a blue marker (Figure 9S). This placement implies that the compound possesses high gastrointestinal absorption potential accompanied by active efflux from the central nervous system and the gastrointestinal lumen [45]. Toxicity assessments for dolichin A reveal inactivity with respect to carcinogenicity, hepatotoxicity, nephrotoxicity, and mutagenicity; however, it exhibits activity concerning immunotoxicity (Table 10S). Additionally, dolichin A falls into category IV with a median lethal dose (LD₅₀) of 500 mg/kg.

4.   Conclusions

This research provides a comprehensive understanding of calcium oxalate (the primary constituent of kidney stones) and dolichin A (a chemical found in M. uniflorum), employing a multifaceted approach that encompasses FTIR spectroscopy, DFT analysis, molecular docking, molecular dynamics simulations, MM-GBSA, druglikeness, and toxicity prediction. The experimental FTIR spectra of the kidney stone samples exhibit distinct peaks at 1609 ± 1 cm⁻¹ and 1314 ± 1 cm⁻¹, which are in close agreement with the calculated FTIR spectrum of calcium oxalate in a water solvent as well as the experimental FTIR spectrum of calcium oxalate monohydrate. Quantum chemical properties including dipole moment, HOMO and LUMO orbital energies, and electrophilicity index suggest that calcium oxalate and dolichin A exhibit their potential for H-bond and nonbonded interactions with lactoferrin. The protein-ligand interactions and conformational stability have been studied using molecular docking and molecular dynamics simulations. Notably, calcium oxalate exhibits a binding energy of −3.86 kcal/mol with lactoferrin. Among the various compounds identified in M. uniflorum, dolichin A emerges as a promising inhibitor for lactoferrin, demonstrating the binding energy of −7.78 kcal/mol. Lactoferrin in complex with dolichin A exhibits dynamical stability within the trajectory of 100 ns. The MM-GBSA calculation after molecular dynamics simulations of lactoferrin-dolichin A complex shows the binding energy of −17.61 ± 4.03 kcal/mol with mostly contributed vdW energy of −37.13 ± 2.44 kcal/mol. These results suggest that dolichin A stands as a potential inhibitor of lactoferrin, making it a prospective candidate for the development of a drug targeting calcium oxalate type kidney stone disease. The preclinical experiments are recommended to validate its inhibitory capabilities.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.