Spectroscopic Investigations, Computational Analysis and Molecular Docking to SAR-Cov-2 Targets Studies of 5,8-Quinolinedione Attached to Betulin Derivatives

Abstract

The 5,8-quinolinedione-betulin hybrids were investigated using spectroscopic methods as well as a variety of quantum chemical calculations in order to characterize their molecular structure. We used FT-IR and NMR spectroscopy supplemented by the density functional theory (DFT) calculations, molecular electrostatic potential (MEP) and molecular orbital (HOMO, LUMO) analyses. The experimental and calculated FT-IR spectra showed a good correlation for all compounds. Analysis of carbonyl band showed that the compounds are the 7-mono substituted. The calculated ¹H NMR and ¹³C NMR spectra of hybrids reproduced well the experimental ones. Identification of C-6 and C-7 carbon atoms of 5,8-quinolinedione revealed the position of betulin moiety at the C-7 of 5,8-quinolinedione. Molecular electrostatic potential maps of hybrids allowed to recognize the electrophilic and nucleophilic regions within the molecules. The molecular docking study was used to examine the interaction between the 5,8-quinolinedione-betulin hybrids and the SARS-CoV-2 protein, like: Mpro and PLpro. The obtained results showed that compounds with the highest Dock Score are good anti-SARS-CoV-2 potential drug candidates.

1. Introduction

Compounds with 1,4-quinone moiety and their derivatives are a large class of compounds found in the plant kingdom and as products of secondary metabolism of some microorganisms and fungi. The occurrence of heterocyclic scaffolds in chemical compounds facilitates interaction with different molecular targets, which cause a high biological effect [1,2]. 5,8-Quinolinedione moiety is a combination of 1,4-quinone moiety (ring A) with the pyridine fragment (ring B) and was found in microorganisms from Streptomyces species (Figure 1). Chemical modification of 5,8-quinolinedione moiety at the C-6 and C-7 positions provides derivatives, which exhibit a wide spectrum of biological activity, including anticancer, antiviral, antibacterial and antimalarial activities [3,4,5,6].

Synthetic derivatives of 5,8-quinolinedione are obtained in the reaction between 6,7-dichlor-5,8-quinolinedione and a nucleophilic substance. Commonly, the product of the reaction is a mixture of two mono-substituted compounds that are difficult to separate. Moreover, the determination of the substituent position in a pure compound is complicated [7,8,9,10]. The ¹³C NMR and FT-IR spectroscopy are useful techniques for determining the substituent position in 5,8-quinolinedione moiety [11,12,13].

Continuing our research on alkoxy derivatives of 5,8-quinolinedione, we decided to study the ¹H NMR and ¹³C NMR and FT-IR spectra for the 5,8-quinolinedione-betulin hybrids. The analysis was supplemented by the density functional theory (DFT) calculations.

Respiratory tract infections are diseases cause by RNA viruses, including respiratory syncytial virus, influenza viruses, rhinoviruses, parainfluenza viruses, human metapneumovirus and human coronaviruses [14]. In late 2019, a new strain of the human coronavirus was identified in China, which was named COVID-19 or SARS-CoV-2 [15,16]. All over the world, thousands of people are infected by this virus every day and many of them die [15]. Scientists all over the world are working on a new drug and treatment for the SARS-Cov-2 disease. Antiviral compounds, such as remdesivir or dexamethasone, are used in the treatment of COVID-19 (Figure 2) [17,18]. Since the late 60s, the option to use chloroquine as an antiviral has been considered in a wide range of diseases [19] (Figure 2). Based on the recent announcements of Gao et al. [20], Wang et al. [21], and Colson et al. [22], chloroquine may be the first successful attempt to use this drug in the treatment of SARS-CoV-2 related diseases.

A library of over 50,000 compounds that could interact with the SARS-CoV-2 protein was analyzed [23,24]. The in silico study indicated that triterpenoids could be an inhibitor of the COVID-19 protein [25,26]. Zhang et al. conducted in silico research on betulinic acid and cryptotanshinone as 3CLpro and PLpro inhibitors [25]. Bevirimat and phosphate derivatives of 3-carboxyacylbetulin (Figure 3) were docked to the active center of the SARS-CoV-2 protein, like: Mpro (also known as 3CLpro), RdPp, E-protein and S-protein, and 2′-O-ribosemethyltransferase (MTase) [27,28]. Wen et al. tested diterpenoids (pinusolidic acid and forskolin) and triterpenoids (betulinic acid and betulonic acid) in vitro as inhibitors of 3CL protease (Figure 3) [26]. Moreover, Chang et al. conducted a molecular docking study of 5,8-quinolinedione derivatives, like DA3003-1 (Figure 3). Chang et al. obtained the crystal structure of HCoV-OC43 bound to the DA300-1 molecule, showing that this compound can be used as an inhibitor against the Nucleocapsid (N) of SARS-CoV-2 [29].

Inspired by literature data, the molecular docking study was used to examine the interaction of Mpro and PLpro proteins with the 5,8-quinolinedione-betulin hybrids.

2. Materials and Methods

2.1. Physical Measurements

Compounds 1–4 (Figure 4) were synthesized according to the methods described previously [30]. The Fourier Transform Infrared Spectroscopy (FT-IR) spectra were recorded using the IRAffinity-1 Shimadzu spectrophotometer. The resolution of spectrum was 0.5 cm⁻¹. The sample of each derivative was prepared as KBr pellet and measured in the range of 400–3500 cm⁻¹ at room temperature. The ¹H and ¹³C NMR spectra were measured using the Bruker Avance 600 spectrometer. The 10 mg of compounds were dissolved in 0.6 mL of CDCl₃.

2.2. Computational Details

The optimized chemistry structure of compounds 1–4 was calculated using the DFT (B3LYP/6-311G+(d.p)) method implemented in the Gaussian 16 program package [31]. All of the local minima of energy were confirmed by the absence of imaginary mode invibrational calculations. In our calculations, we applied the basis set with the diffuse function on heavy atoms (+) to obtain a better description of lone pair electrons with orbitals occupying a larger region ofspace. The geometry optimization and frequency calculation for a bigger basis set (6-31++G(d.p) and 6-311+G(d.p)) were also applied but significant changes in geometry parameters and fundamental frequencies were not observed.

The obtained geometries of compounds 1–4 were used for the calculated the IR, ¹H and ¹³C NMR spectra. The harmonic vibrational frequencies were calculated at the same level of theory (B3LYP/6-311G+(d.p)) and the theoretical wavenumbers were scaled by a factor of 0.964 [32]. The NMR spectra were calculated using the Gauge-Independent Atomic Orbital (GIAO) method implemented in the Gaussian 16 program package [31,33]. The electrostatic potential V(r) and the electron density maps, which define the molecular surfaces, were determined [34]. All obtained results were visualized using the GaussView, Version 5 software package [35].

2.3. Molecular Docking Study

The three-dimensional structures of studied compounds were generated in their low-energy conformation using the Gaussian 16 (revision A.03) computer code [31] in the density functional theory (DFT, B3LYP) and 6-311+G(d.p) basis sets. The target macromolecule for molecular docking studies was obtained from the Protein Data Bank (https://www.rcsb.org/). We used crystal structures of SARS-CoV-2 main protein (PDB ID: 5R7Z) and papin-like protease (PDB ID: 6W9C).

In addition, AutoDock Vina [36] was used for the docking analysis. The region of interest used for Vina docking was fixed as X = −33,784, Y = 20,933, Z = 33,306 for papain-like protease and X = −11,631, Y = 2201, Z = 23,194 for the COVID-19 main protein. The volume was set as 25 × 25 × 25 Å. After calculations, only the 9 highest-scored poses were returned as a docking result for ligand-cavity configuration. All the obtained results were ranked according to their score values and presented in kcal/mol. Molecular docking details were visualized using the BIOVIA Discovery Studio virtual environment [37].

MD simulations for the modeled protein-ligand complexes were performed using the program NAMD (NAnoscale Molecular Dynamics program; v 2.13) [38] and all files were generated using visual molecular dynamics (VMD) [39]. Parameterization of the ligands was performed in the CGenFF server. The parameterized ligands were inserted into the protein and saved in the form of a protein-ligand complex by VMD with the binding pocket residues. Then, the protein-ligand complex was immersed in the center of a 155.3 Å box of water molecules for Mpro and 147.8 Å box of water molecules for PLpro with a TIP3P water box. The 0.15 M ions (Na⁺ and Cl⁻) were added to provide charge neutralization and electrostatic screening. CHARMM (Chemistry at HARvard Macromolecular Mechanics) 36 parameter file for proteins and lipids; phi and psi cross-term map corrections were used in the force field for proteins with similar chemical structures. For the minimization and equilibration of the complexes in the water box, we assumed force-field parameters excluding scaling of 1.0 Å. All atoms, including those of hydrogen, were illustrated explicitly. The hydrogen atom coordinates of proteins were generated using the VMD Tk-Console salvation command. Integrator parameters also included 2 fs/step for all rigid bonds and nonbonded frequencies were selected for 1 Å. Full electrostatic evaluations for 2 Å were used, with ten steps for each cycle. The complexes’ preliminary energy was minimized using 2000 steps of the Powell algorithm at a constant temperature (310 K), followed by simulation of an additional 144,000 steps with Langevin dynamics to control the kinetic energy, temperature, and/or pressure of the system. Finally, the solvated protein-ligand complex system was equilibrated with 50,000 minimization steps and production runs with 5,000,000 steps (10 ns) were than performed.

Results were analyzed using the Visual Molecular Dynamics package (VMD, https://www.ks.uiuc.edu/research/vmd/) [40].

3. Results and Discussion

3.1. Molecular Geometry

Molecular structures with numbering of atoms for compounds 1–4 are presented in Figure 4. Tables S1–S4 present the calculated bond lengths and angles for 1–4. The bond lengths of the C-C in betulin moiety are typical for a single bond, while the bond lengths of C-C, C-O and C-N in 5,8-quinolinedione moiety are intermediate between single and double bond lengths. These results also showed the π delocalisation for the 5,8-quinolinedione ring [41,42]. Comparing the orientation of 5,8-quinolinedione in relation to betulin moiety, it can be seen that they have a similar arrangement in the following pairs 1–2 and 3–4 (Tables S1–S4).

3.2. FT-IR Spectra

The both experimental and calculated FT-IR spectra for compounds 1–4 are presented in Figure 5 and Figure 6 and their assignments are entered in

Table 1. Band assignments of experimental and calculated IR spectra for studied compounds 1–4.

Assignment	1		2		3		4
	Exp	Calc	Exp	Calc	Exp	Calc	Exp	Calc
v C≡C-H							3330–3258	3361
v OH	3418	3621
v aliphatic CH3, CH2, CH	3073–2870	3119–2841	3070–2872	3119–2942	3073–2870	3120–2856	2945–2870	3118–2858
v C≡C							2115	2138
v OC=O			1728	1777	1736	1785	1712	1747
vas C=O	1691	1704	1691	1700	1693	1700	1668	1704
vs C=O	1674	1688	1670	1690	1668	1685	1639	1687
vs C=C isopropenyl	1637	1670	1641	1669	1641	1670	1637	1666
v aromatic C-C, C-H	1593–1566	1583–1557	1593–1562	1584–1533	1589–1558	1577–1546	1591–1558	1574–1548
σ aliphatic CH3, CH2, CH	1458	1479	1465	1477	1462	1448	1463	1478
ν CH3, CH2	1388–1311	1392–1288	1373–1315	1369–1291	1369–1314	1373–1292	1388–1311	1409–1288
ν C-N	1278	1288	1280	1291	1278	1291	1274	1288
δ C-C, C-H	1246	1235	1251–1240	1230	1251	1238–1222	1252–1224	1236–1201
ν C-C, C-H	1182–1149	1166–1132	1182–1149	1168–1130	1184–1148	1168–1130	1182–1148	1166–1131
δ C-C, C-H	1096–1064	1053	1093–1022	1052–1006	1064–1035	1153–1052	1093–1043	1056–1029
ν C-O	1045–1020	1027–1010	1024	1021–1007	1033–1020	1037–1006	1018	1018–1001
ν C-O	945–927	923	945–933	921	910	916	974	958
δ C-C, CH3, CH2, CH	880–856	880–854	900–852	875–852	883–856	875–862	883–856	880–862
v C-Cl	815	828	813	829	812	829	812	828
v C-C C-H	754–725	741–706	734–704	736–705	726–690	737–705	754–723	732–706

. The band assignment of experimental FT-IR spectra was carried out on the basis of theoretically predicted vibrations by DFT method. Generally, for all compounds, the theoretical spectra reproduce the experimental ones well.

3.2.1. C-H and C-C Vibration

Compounds 1–4 contain aromatic and aliphatic C-H groups, which exhibit stretching vibrations in the ranges of 3073–2870 cm⁻¹ and 1388–1288 cm⁻¹. The absorption peak at the 1641–1637 cm⁻¹ is attributed to the C=C stretching vibration of the isopropenyl group at C-19 position of betulin moiety. Stretching vibration in the range of 1593–1558 cm⁻¹ refers to C-C and C-H vibrations of the aromatic group of 5,8-quinolinedione moiety (Figure 5). The strong band at 1478–1448 cm⁻¹ is assigned to the deformation vibration of C-H in the methyl groups of betulin. The C-H and C-C bend vibrations in the aliphatic group of betulin moiety occur in the ranges of 1251–1064 cm⁻¹ and 1252–1052 cm⁻¹ for experimental and calculated spectra, respectively.

Additionally, for derivative 4, the band resulting from the vibration of the acetylenic group at the C-28 position of betulin moiety is observed. According to the literature [43], the double peak in the region of 3298–3277 cm⁻¹ is assigned to the C≡CH stretch vibration. Moreover, the band at 2121 cm⁻¹ is assigned to the C≡CH aliphatic stretching vibrations (Figure 5d).

3.2.2. C-O and O-H Vibration

According to the literature [43,44], the absorption peaks in the range of 1045–1000 cm⁻¹ refer to the asymmetric stretching vibrations of C-O groups. In the calculated spectra, two separate peaks relating to the vibration of the C-O-C group at the C-28 and C-3 positions of betulin are observed in this region. In the experimental spectra of 1 and 3, two separate peaks are observed, while, for 2 and 4, broad peak at 1024 cm⁻¹ and 1018 cm⁻¹ are observed, respectively. Furthermore, the symmetric stretching vibrations of the C-O group at 980–910 cm⁻¹ are observed [43,44]. The symmetric stretching vibrations of the C-O-C group between betulin and 5,8-quinolinedione moieties at the 958–916 cm⁻¹ and 974–916 cm⁻¹ in calculated and experimental spectra are observed, respectively. In experimental spectra of compounds 1 and 2 containing 5,8-quinolinedione moiety at the C-28 position of betulin, peaks of the C-O-C vibration split into two separated peaks. According to the literature, this observation may suggest the formation of weak hydrogen bonds [43,44,45].

Additionally, compound 1 has a hydroxy group at the C-3 position of betulin, for which the absorption peaks are observed at 3418 cm⁻¹ and 3621 cm⁻¹ in experimental and calculated spectra, respectively.

3.2.3. C=O Vibration

Compound 1 contains two carbonyl groups at the C-5 and C-8 positions of 5,8-quinolinedione moiety (Table 1). According to the literature, in the region of 1638–1704 cm⁻¹, the asymmetric and symmetric vibrations of the carbonyl group at C-8 and C-5 position of 5,8-quiolinedione moiety can be observed, respectively (Figure 6a) [7,11,46]. Comparison of experimental and calculated spectra shows that experimental bands are shifted towards lower wavenumbers. As seen in Table 1, the frequency separation ∆ν = ν_as − ν_s of the carbonyl bands in the experimental and calculated spectra for 1 are equal to 17 cm⁻¹ and 36 cm⁻¹, respectively. As can be found in the literature [11,12], the splitting of the carbonyl band indicates that the compounds are 7-mono substituted.

The molecular structure of hybrid 2 differs from 1 in terms of the presence of the acyl group at the C-3 position of betulin moiety (Table 1). In FT-IR experimental and calculated spectra of 2, the carbonyl and carboxyl bands are observed in the ranges of 1670–1728 cm⁻¹ and 1690–1777 cm⁻¹, respectively. In the experimental spectrum, the band of the carboxyl group of betulin moiety is shifted towards lower frequencies compared to the calculated one, i.e., from 1777 cm⁻¹ to 1728 cm⁻¹ (Figure 6b). In the experimental spectrum, the position of carbonyl bands in 5,8-quinoliedione are similar to those in hybrid 1 and the frequency separation is equal to 26 cm⁻¹. In the calculated spectrum, the carbonyl bands are close to each other and the ∆ν is equal to 10 cm⁻¹ (Table 1, Figure 6).

Hybrids 3 and 4 contain 5,8-quinolinedione moiety at the C-3 position of betulin and acyl or propynoil group at the C-28 position. The calculated spectra of both compounds are similar, which means that the stretching vibration of the carboxyl group is attributed to peak at 1747 cm⁻¹, while the asymmetric and symmetric stretching vibrations of carbonyl groups are observed in the range of 1666–1687 cm⁻¹. In experimental spectra of both compounds, the carboxyl group band is observed at 1712 cm⁻¹ (Figure 6c,d). However, in the spectrum of 4 the carboxyl peak is broad, which could suggest the formation of a hydrogen bond.

Introduction of carboxyl group at C-28 position of betulin moiety influences the shift of the carbonyl band toward lower frequencies. For hybrid 3, the experimental frequencies of the carbonyl bands are similar to those for 2, i.e., 1668 cm⁻¹, 1693 cm⁻¹ and 1670 cm⁻¹, 1691 cm⁻¹, respectively (Figure 6b,c). As seen in Figure 6d, for compound 4, the asymmetric and symmetric stretching vibrations are localized at 1668 cm⁻¹ and 1639 cm⁻¹, respectively. For hybrid 4, the experimental and calculated frequency shifts ∆ν are equal to 29 cm⁻¹ and 21 cm⁻¹, respectively.

3.2.4. C-N Vibration

The C-N band vibration in the FT-IR spectra is difficult to identify because these bands are observed at the same region as the C-C and C-H bands of aromatic compounds [13,43,44]. According to literature, for aromatic compounds, the C-N stretching vibration is observed in the range 1266–1382 cm⁻¹ [44]. In the presented study, the band at 1274–1280 cm⁻¹ and 1288–1291 cm⁻¹ in the experimental and calculated spectra are attributed to the stretching C-N vibration (Table 1).

3.2.5. C-Cl Vibration

According to literature data, the C-Cl stretching vibrations show strong bands in the region of 760–505 cm⁻¹ [43,44]. In our experiment, the computed and experimental C-Cl bands are observed at 812–815 cm⁻¹ and 828–829 cm⁻¹, respectively (Table 1).

3.3. ¹H and ¹³C NMR Study

Hybrids 1–4 were characterized by ¹H and ¹³C NMR spectra and the corresponding chemical shifts were simulated by the GIAO method using the optimised limitations from the 6-311+G(d.p) basis set [45,46]. The hybrids 1–4 contain more than 50 protons and 30 carbon atoms. Of these, the more important protons and carbon atoms were selected for analysis. The both experimental and calculated chemical shifts for the selected peak of compounds 1–4 are presented in

Table 2. The experimental and calculated chemical shifts in ¹H NMR and ¹³C NMR spectra for selected atoms of hybrids 1–4.

Atoms	1		2		3		4
	Exp	Calc	Exp	Calc	Exp	Calc	Exp	Calc
H-2q	9.07	9.19	9.06	9.11	9.04	9.13	9.04	9.09
H-4q	8.50	8.35	8.50	8.23	8.50	8.25	8.50	8.26
H-3q	7.73	7.36	7.73	7.31	7.71	7.27	7.72	7.27
H-28Abet	4.87	4.79	4.88	4.70	4.28	5.14	4.40	4.17
H-29Abet	4.71	5.09	4.71	5.01	4.71	4.99	4.72	4.87
H-29Bbet	4.62	4.72	4.62	4.60	4.61	4.52	4.62	4.92
H-28Bbet	4.29	3.94	4.28	3.82	3.88	3.28	4.00	3.76
H-3bet	3.21	3.50	4.49	4.76	4.81	4.78	4.81	4.74
H-19bet	2.44	2.48	2.43	2.49	2.47	2.34	2.45	2.38
CH3C(O)O	-	-	2.06	2.02	2.09	2.05	-	-
CHCC(O)O	-	-	-	-	-	-	2.91	2.57
C-5q	178.2	181.8	178.2	181.5	178.1	181.9	178.1	181.9
C-8q	177.7	181.6	177.8	180.6	177.9	181.3	177.6	181.4
C-7q	158.1	156.1	158.2	155.3	158.3	156.9	158.1	156.9
C-2q	154.6	155.3	154.7	153.9	154.8	153.1	154.5	153.7
C-20bet	150.1	154.2	150.1	153.2	150.1	152.6	153.3	153.5
C-8Aq	146.6	143.4	146.8	143.2	146.6	143.5	149.9	143.0
C-4q	134.8	132.4	134.9	131.6	134.8	131.7	134.9	131.7
C-6q	131.5	137.7	129.0	137.4	129.1	137.8	128.9	138.1
C-4Aq	128.2	125.6	128.2	124.7	128.2	124.8	128.1	125.2
C-3q	128.1	125.1	127.8	124.5	127.8	124.6	127.8	124.8
C-29bet	110.1	110.8	109.9	109.8	109.9	109.8	110.0	110.4
C28bet	79.0	80.1	80.8	81.2	73.9	64.6	74.6	72.1
C-3bet	73.9	76.8	74.0	79.1	80.9	99.2	92.6	98.7
C(O)O	-	-	170.9	180.3	171.0	180.2	169.4	162.6
CH3C(O)O	-	-	16.1	15.5	15.8	15.2	-	-
CHCC(O)O	-	-	-	-	-	-	74.8	75.6
CHCC(O)O	-	-	-	-	-	-	74.7	72.2

. Figures S1–S8 present the experimental and calculated ¹H and ¹³C NMR spectra for all derivatives 1–4.

The calculated ¹H NMR spectra of compounds 1–4 reproduce the experimental ones well (Figure 7a) and the correlation coefficients are in the range of 0.986–0.995.

Comparison of experimental and calculated ¹H NMR spectra of 1–4 shows the differences in values of proton signals in the range of 4.17–5.15 ppm. According to literature data [47,48,49], the signals of protons at the H-29 position of betulin are singlet, while signals of H-28 protons are dublet. The analysis of experimental ¹H NMR shows that H-29A and H-29B protons are observed at 4.28–4.87 ppm and 4.62–4.61 ppm, respectively, while two dublet signals of H-28 (A and B) protons are observed at 4.88–4.40 ppm and 3.88–4.29 ppm, respectively. In the calculated ¹H NMR spectra, the signal of the H-29A proton is at 4.87–5.09 ppm, while H-28A is at 4.17–5.14 ppm. For all compounds, the calculated signal of H-29A is at higher chemical shift than H-28A (Table 2).

The 5,8-quinolinedione moiety contains only three protons in the pyridine ring, which prevents the use of 2D NMR spectra to analyse the carbon atoms in heterocyclic moiety. The carbon atoms in 5,8-quinolinedione moiety were assigned signals on the basis of the comparison of the experimental and calculated spectra (Table 2). The calculated ¹³C NMR spectra reproduce the experimental ones well and their correlation coefficients are in the range 0.989–0.996 for 1–4 (Figure 7b). Identification of the chemical shift for carbon atoms at C-6 and C-7 of 5,8-quinolinedione moiety allows to determine the position of the substituent in 5,8-quinolinedione. In the experimental spectra of compounds 1–4, the signal of C-6 carbon in 5,8-quinolinemoiety is shifted toward lower chemical shifts than in calculated ones. The signals of C-7 carbon are observed at 158.1–158.3 ppm and 155.3–158.3 in the experimental and calculated spectra, respectively. The good reproducibility of these two signals in experimental and calculated spectra confirms that hybrids 1–4 contain betulin moiety at the C-7 position of 5,8-quinolinedione moiety.

3.4. Frontier Molecular Orbitals

Frontier molecular orbitals theory describing HOMO and LUMO orbitals is useful for determining the energy distribution and energetic behaviour of a molecule. The HOMO and LUMO energies denote the ability of molecules to donate and acquire an electron, respectively [50]. The energy gap (E_HOMO-E_LUMO) makes it possible to designate the kinetic stability and reactivity of a compound [13,46,50,51]. The HOMO and LUMO orbitals and its energy gap was calculated in the Gaussian09 package at the B3LYP level using the 6-311G+(d.p) basis set [31].

As seen in Figure 8, the isocontour of the HOMO orbital of all compounds are mainly delocalized on the isopropenyl group and E ring of betulin moiety. The energy of the HOMO orbital (E_HOMO) is within the range of (−6.605 to −6.839) eV meaning that the molecules are low reactive with the electrophilic molecule [50]. The LUMO orbital is delocalized on 5,8-quinolinedione moiety and its energy shows that compounds 1–4 are high reactive against the nucleophilic molecule [50]. The negative energy of HOMO and LUMO orbitals shows that the compounds are stable. For hybrids 1–4, the large HOMO–LUMO energy gap varies in the range of (−3.015 to −3.125) eV, and this can be related to some molecular hardness and less polarizability. According to literature data [52,53], the van der Waals interactions between two molecules depends on energy gap ΔE and its value could be correlated with the biological activity of the compound. The ΔE energy gap for compounds 1–4 has low value, suggesting possible interactions between hybrids and proteins responsible for their biological activity (

Table 3. Calculated HOMO and LUMO orbital energies and global reactivity descriptors for compounds 1–4.

Parameters	6-311G+(d.p)
	1	2	3	4
Fourier molecular orbital energies (eV)
EHOMO	−6.605	−6.839	−6.826	−6.756
ELUMO	−3.585	−3.783	−3.702	−3.741
ΔE = EHOMO − ELUMO	−3.020	−3.056	−3.125	−3.015
Global reactivity descriptors (eV)
Ionization potential (I)	6.605	6.839	6.826	6.756
Electron affinity (A)	3.585	3.783	3.702	3.741
Hardness (η)	1.510	1.528	1.562	1.508
Softness (s)	0.662	0.654	0.640	0.663
Chemical potential (μ)	−5.095	−5.311	−5.264	−5.249
Electronegativity (ϰ)	5.095	5.311	5.264	5.249
Electrophilicity index (ω)	8.598	9.230	8.868	9.137

, Figure 8).

The energies of the HOMO-LUMO orbital facilitate calculation of the global reactivity descriptors, such as: ionization potential (I), electron affinity (A), hardness (η), softness (s), chemical potential (µ), electronegativity (χ) and electrophilicity index (ω) [54]. Global reactivity descriptors are calculated according to equations: I = −EHOMO, A = −ELUMO, η = (I − A)/2; s = 1/η; μ = −(I + A)/2, χ = −μ, ω = μ²/2η and present in Table 3.

The calculated reactivity descriptors show high stability of hybrids 1–4, signifying their resistance toward deformation of the electron cloud of the chemical system under small perturbations and less polarizability.

3.5. Molecular Electrostatic Potential

The molecular electrostatic potential (MEP) provides knowledge about the arrangement of the positive and negative charge of the molecule. The colours on the MEP surface represent different values of the electrostatic potential. The most negative and positive regions are represented by red and blue, respectively [55]. The MEP maps of hybrids 1–4 were calculated using the DFT/B3LYP level using the 6-311G+(d.p) basis set and they are sketched for an order of −44.3774 kcal/mol to 44.3774 kcal/mol. The calculated MEP maps for hybrids 1–4 are presented in Figure 9a–d.

For all compounds 1–4, the nucleophilic regions (red colour) are localized in three main areas. The first area contains a nitrogen atom, carbonyl group at the C-8 position and an oxygen atom at C-7 position of 5,8-quinolinedione moiety. The second and third area includes the carbonyl group at the C-5 position, and the hydroxy (1) or carboxyl group (2–4) of betulin moiety, respectively (Figures S9–S12).

For all hybrids, in the first area, following two potential minima are observed: −(2.72–2.88) eV and −(2.07–2.18) eV. In second area near carbonyl group at C-5 position of 5,8-quinolinedione moiety the potential minimum is equal −(2.18–2.29) eV. For hybrid 1, in third area only one potential minimum (−2.99 eV) can be observed. For hybrid with acyl group (2–4) two potential minima: −(2.19–2.45) eV and −(1.63–1.74) eV at this areas are observed. The negative region could be involved in the formation of a hydrogen bond with some biological targets.

The electrophilic regions (blue colour) are localized near hydrogen atoms of 5,8-quinolinedione. Additionally, in hybrid 4, the positive region is situated on a hydrogen atom of propynoil group. The light green colour region containing the betulin scaffold is charge neutral.

3.6. Molecular Docking

In this study, we use the AutoDock Vina program (referred to as Vina) in in silico research. As a molecular target, we use two main SAR-CoV-2- polymerase, named Mpro and PLpro. We used dexamethasone, remdesivir, and chloroquine as reference ligands. The tested compounds ranked by AutoDock Vina are shown in

Table 4. Vina affinity scoring values of tested compounds.

Compound	ΔG [kcal/mol]
	Mpro	PLpro
1	−8.7	−8.1
2	−8.5	−7.9
3	−8.5	−7.7
4	−8.9	−7.6
Dexamethasone	−7.4	−6.5
Remdesivir	−7.5	−7.2
Chloroquine	−5.6	−5.6

. The lowest scores of binding energy (kcal/mol) of protein-ligand complexes correspond to a strong binding affinity, and the most probable ligand-protein system in vivo.

Results obtained in the Vina program indicate that all tested compounds show lower binding energy with Mpro and PLpro compared to the reference dexamethasone, remdesivir and chloroquine (Table 4).

As in other coronaviruses, the main protease plays an important role in viral maturation by processing many polyproteins that are translated from the viral RNA.

According to the crystallographic data, amino acids His 41, His 164, Met 49, Met 165, Thr 190, and Gly 143 play an important role in stabilizing the ligand-Mpro complexes [56,57]. Figure 10a and Figure 11 and

Table 5. Interaction of hybrids 4 and 1 with active sites of Mpro and PLpro, respectively.

Protein		Ligand		Interaction
Name	Residue	Name	Residue	Type	Distance [Å]
Mpro	Gly143	4	carbonyl oxygen at C-8q	conventional hydrogen bond	2.40
	Gly143		carbonyl oxygen at C-8q	carbon hydrogen bond	2.85
	Thr26		pyridine ring	π-donor	3.11
	Thr25		1,4-quinone ring	π-sigma	2.92
	Cys145		C-2	alkyl-alkyl	4.72
	Pro168		cyclopentane ring	alkyl-alkyl	5.13
	Cys145		C-24	alkyl-alkyl	4.43
	Cys145		C-25	alkyl-alkyl	5.20
	Met49		1,4-quinone ring	alkyl-alkyl	5.13
	Met165		C-29	alkyl-alkyl	5.10
	Leu167		C-29	alkyl-alkyl	5.48
	His41		C-2	π-orbitals	5.36
	His41		C-8q	π-orbitals	4.99
PLpro	Glu167	1	carbonyl oxygen at C-8q	carbon hydrogen bond	2.77
	Glu167		pyridine ring	π-anion	3.78
	Met206		C-7	alkyl-alkyl	5.24
	Met206		C-16	alkyl-alkyl	5.19
	Val202		C-27	alkyl-alkyl	5.06
	Tyr207		C-23	π-alkyl	4.74
	Tyr207		C-24	π-alkyl	4.58

present the possible interaction of the best scored compound 4 inside the binding pocket of Mpro.

Figure 11 presents the possible interaction of compound 5 inside the binding pocket of Mpro after 2D analysis using the Discovery Studio Visualizer.

The corresponding amino acids that are significantly involved in the hydrophobic interactions are as follows: Thr25 (π-sigma), Thr26 (π-donor), His41 (π-orbitals), Met49, Cys145, Met165, Leu167, Pro168 (alkyl-alkyl). Generally, strong hydrogen bond interaction (with donor-acceptor distance of 2.40 Å) between Gly143 and carbonyl oxygen atoms increases the stability of the ligand-receptor complex.

Molecular dynamics (MD) calculation provides a better understanding of structure-function relationships in motion and other conformational changes by the proteins. Therefore, an MD study was conducted in this paper to evaluate the dynamic behavior of the top-scoring complexes for a period of 10 ns. Analyzing the complexes ligand root-mean-square deviation (RMSD) and hydrogen bonding can provide insights into a structural conformation faced by the protein and ligand during simulation (Figure 12, Figure 13 and Figure 15).

Mpro has the most optimal RMSD profile with a medium value below 2 Å (average value 1.448 Å). Hybrid 4, during the first 7 ns of simulation, maintains a constant value, but then it quickly raises RMSD to a value of 0.9 Å, and stabilizes at an RMSD value of 0.5 Å (average value 0.506 Å).

The data shown in Figure 13 shows the results of analysis of hydrogen bond contacts in the studied complexes. The compound 4 forms two hydrogen bonds, and the binding with Gly146 residue occurs in the complex at the 28% time of measurement.

PLpro is responsible for processing non structural proteins Nsp1, Nsp2 and Nsp3, which are released after the cleavage of the N-terminus of the replicas polyprotein. In addition, the binding site contained S3/S4 pockets (residues Asp164, Val165, Arg166, Glu167, Met 208, Ala246, Pro247, Pro248, Tyr 264 and Gly266) [58].

Figure 10b and Figure 14 and Table 5 presents the possible interaction of the top docked pose of 1 inside the binding pocket of COVID-19 PLpro. The corresponding amino acids that are significantly involved in the hydrophobic interactions are as follows: Glu167 (carbon hydrogen bond, π-anion), Val202, Met206 (alkyl-alkyl), Tyr207 (π-alkyl).

The plots for PLpro Cα versus time for the simulation with compound 1 are shown in Figure 15. The visual analysis of the trajectory confirms the stability of the structure of the protein at around 6 ns. The compound 1 has the most optimal RMSD profile with a medium value below 1 Å (average value 0.727 Å). These results indicate, that compounds possessing the highest Dock Score (1 and 4) are good anti-SARS-CoV-2 drug candidates.

4. Conclusions

The various spectroscopic investigations and computational analysis enabled the examination of the 5,8-quinolinedione-betulin hybrids 1–4 as candidates for interaction with SARS-CoV-2.

The experimental and calculated FT-IR spectra demonstrated a good correlation for all compounds. The most significant difference in the FT-IR spectra of hybrids 1–4 occurred in the region of the carbonyl bands. The type of substituent impacted the frequency position and separation of the carbonyl bands ∆ν.

The computed HOMO and LUMO energy gap supported the possibility of hybrids to be chemically active, as they are high reactive against nucleophilic molecules. The molecular electrostatic potential (MEP) maps showed that the negative potential sites are present in the electronegative atoms, i.e., nitrogen and oxygen atoms.

The molecular docking study showed that hybrids 1–4 were locked onto the COVID-19 proteins with a lower negative dock energy than the reference compounds. The lower negative dock energies were obtained for the Mpro–hybrid 4 and PLpro–hybrid 1 complexes. Using the molecular dynamics simulation for these two complexes, the stability of the docking pose was determined. On the whole, the obtained results show that compounds possessing the highest Dock Score are good potential anti-SARS-CoV-2 drug candidates. However, further experimental studies are necessary to validate these findings.