Identification of 1,2,3-triazole-phthalimide derivatives as potential drugs against COVID-19: a virtual screening, docking and molecular dynamic study.

In this work we aimed to perform an in silico predictive screening, docking and molecular dynamic study to identify 1,2,3-triazole-phthalimide derivatives as drug candidates against SARS-CoV-2. The in silico prediction of pharmacokinetic and toxicological properties of hundred one 1,2,3-triazole-phtalimide derivatives, obtained from SciFinder® library, were investigated. Compounds that did not show good gastrointestinal absorption, violated the Lipinski's rules, proved to be positive for the AMES test, and showed to be hepatotoxic or immunotoxic in our ADMET analysis, were filtered out of our study. The hit compounds were further subjected to molecular docking on SARS-CoV-2 target proteins. The ADMET analysis revealed that 43 derivatives violated the Lipinski's rules and 51 other compounds showed to be positive for the toxicity test. Seven 1,2,3-triazole-phthalimide derivatives (A7, A8, B05, E35, E38, E39, and E40) were selected for molecular docking and MFCC-ab initio analysis. The results of molecular docking pointed the derivative E40 as a promising compound interacting with multiple target proteins of SARS-CoV-2. The complex E40-Mpro was found to have minimum binding energy of -10.26 kcal/mol and a general energy balance, calculated by the quantum mechanical analysis, of -8.63 eV. MD simulation and MMGBSA calculations confirmed that the derivatives E38 and E40 have high binding energies of -63.47 ± 3 and -63.31 ± 7 kcal/mol against SARS-CoV-2 main protease. In addition, the derivative E40 exhibited excellent interaction values and inhibitory potential against SAR-Cov-2 main protease and viral nucleocapsid proteins, suggesting this derivative as a potent antiviral for the treatment and/or prophylaxis of COVID-19.Communicated by Ramaswamy H. Sarma.

Introduction

Nowadays respiratory diseases caused by virus are among the main causes of morbidity and mortality due to their easy spreading and the high pathogenicity of viral agents (Piñana et al., 2020). At the end of 2019 a new Coronavirus, called SARS-CoV-2, appeared in Wuhan, China and became one of the most aggressive pandemics in history (Meo et al., 2020), exceeding 3 million infected individuals and causing thousands of deaths worldwide by the beginning of May 2020. The clinical manifestations of Coronavirus disease 2019 (COVID-19) are characterized by fever, fatigue, cough, which can rapidly progress into Acute Respiratory Discomfort Syndrome (ARDS) and death (Huang et al., 2020). Although vaccinal candidates against COVID-19 are already under development stage, a safe and effective vaccine for large-scale vaccination is far from becoming a reality. Hydroxychloroquine or chloroquine alone or in combination with a second-generation of macrolide, have been preconized for the treatment of COVID-19 despite no conclusive evidence of their benefits (Mehra et al., 2020). In this regard, the search for a safe and effective treatment against the severe form of COVID-19 is still necessary (Zhang et al., 2020).

Among the arsenal of bioactive compounds, the 1,2,3-triazole group has been attracting interest of scientific community due to their wide-ranging and versatile medicinal applications such as antivirals (Sun et al., 2020), antibacterial (Gonzaga et al., 2013), antifungal (Thanh et al., 2019) and antiparasitic (Leite et al., 2018) agents. Furthermore, 1,2,3-triazoles can be easily synthesized by low-cost copper catalysed reactions between alkynes and azides, namely click chemistry. Besides being used to synthesize triazoles derivatives, click chemistry is also widely applied for the fusion of two or more chemical groups in a single hybrid molecule with improved biological activity (Ouyang et al., 2018). As triazole derivatives, the phthalimide group has been shown to be useful in the formation of hybrid compounds with antiviral activity against HIV (Al-Masoudi et al., 2016), cytomegalovirus and varicella-zoster (Mandić et al., 2020). Importantly, both triazole and phthalimide derivatives also have relevant anti-inflammatory activity (Assis et al., 2019). Therefore, 1,2,3-triazole-phthalimide hybrids may not only be able to prevent viral replication but also treat the severe form of COVID-19, characterized by an intense inflammatory process called “cytokine storm” (Zhang et al., 2020).

One of the strategies for the discovery of antiviral candidates against SARS-CoV-2 is the association of in silico screening-based approaches with the molecular docking of selected molecules against SARS-CoV-2 target proteins, to identify compounds having good drug-likeness properties. Because experimental approaches for drug discovery are highly costly and time-consuming, screening-based methods offer a fast and low-cost way for prospecting new potential drugs (Hall & Ji, 2020). The SARS-CoV-2 viral genome encodes for 29 proteins—many of which could serve as potential targets for anti-viral drugs. Among these proteins, spike protein (S), C-like proteases and nucleocapsid proteins have gained growing interest due to their relevant roles in the host cell recognition, viral replication, processing and formation of nucleocapsid structure. Spike glycoproteins, found on the surface of SARS-CoV-2, are a class I fusion proteins which bind to human Angiotensin Converting Enzyme 2 (ACE2) receptor mediating both the fusion of the viral envelope with the host-cell membrane and the virus entry into host cells. The chymotrypsin-like cysteine protease, also known as main protease (Mpro) or 3-C like protease (3CLPro), is also an interesting target for drugs against COVID-19. This protease cleaves the C-terminal sequence of the polyprotein PP1A and PP1AB into 16 functional non-structural proteins (nsps), which play an important role in viral replication (Chen et al., 2020). Besides proteases and spike proteins, the nucleocapsid protein of SARS-CoV-2 has been investigated as a target protein for drug development. The nucleocapsid protein comprises two structured domains, the N-terminal domain (NTD) and the C-terminal domain (CTD), interspersed by three intrinsically disordered regions (IDRs) known as the N-terminal arm (N-arm), the central linker region (LKR), and the C-terminal tail (C-tail) arm (Zhou et al., 2008, 2020). The main function of these proteins is to form complexes with ribonucleoproteins. However, they also regulate the replication, transcription of viral RNA and, in the host cells, inhibit protein translation, leading to disruption of host cell metabolism. Hence, targeting this protein can lead to blockbuster therapeutic agents for COVID-19 treatment (Bhatia, Narang & Rawal, 2020; Zeng et al., 2008). In this study, we have performed docking study of selected 1,2,3-triazole-phthalimide derivatives against COVID-19 targets.

Materials and methods

Selection and preparation of ligands

The access to the chemical structures of 1,2,3-triazole-phthalimide hybrid compounds was performed using the SciFinder® database (https://scifinder.cas.org). For this study, only compounds that were obtained by click chemistry and have biological activity already reported in the literature were chosen for the analysis. The chemical structure of compounds was built using ACD/ChemSketch 2.1 and optimized using the PM3 semi-empirical method (Stewart, 1989) implemented in Mercury 4.3.1 (Macrae et al., 2008). The protonation states of the ligands were verified by MarvinSketch 20.19.0 software assuming a physiological pH (Butcher et al., 2015). The canonical smiles were obtained for ADMET in silico analysis (Karwath & Raedt, 2006).

In silico analysis of physico-chemical, pharmacokinetic and toxicological profile (ADMET) of ligands

The analysis of the pharmacokinetic and druglike properties of 1,2,3-triazole-phthalimide derivatives was carried out using the Web server SwissADME, http://www.swissadme.ch/ (Swiss Institute of Bioinformatics). Only the compounds that have good gastrointestinal adsorption and did not violate the following parameters of Lipinski’s Rule of Five: (i) clogP ≤ 5; (ii) Molecular weight (MW) ≤ 500 g/mol; (iii) Number of hydrogen bond acceptors (HBA) (sum of N and O atoms) ≤ 10, (iv) Number of hydrogen bond donors (HBD) (sum of OH and NH groups) ≤ 5, were selected (Lipinski, 2004). The number of rotatable bonds (nRotb) ≤ 10, were also considered for the further toxicity and molecular docking analysis. The prediction of toxicological potential of these compounds was performed using the pKCSM on-line platforms (Pires et al., 2015). The AMES toxicity, hepatotoxicity, and immunotoxicity parameters were also used to filter out compounds.

Molecular docking of 1,2,3-triazole-phthalimide hybrids on SARS-CoV-2 target proteins

The molecular docking of selected compounds on SARS-CoV-2 target proteins was carried out using iGemDock 2.1 and AutoDock 4.2 software. The main protease (Mpro), spike (SP) and viral nucleocapsid (NCP) proteins were chosen according to the stereochemical parameters provided by Ramachandran analysis, using the PROCHECK server (Laskowski et al., 2006). Molecular docking was performed based on the confrontation between optimized ligands structure and biomacromolecules from SARS-CoV-2. The crystal structure of viral proteins was obtained from the Protein Data Bank (PDB) under the codes: spike: 6VSB, in its open (6VYB) and closed states (6VXX) (Walls et al., 2020; Wrapp et al., 2020); nucleocapsid protein, 6VYO (Yadav et al., 2020) and Mpro, 5R80 (Razzaghi-Asl et al., 2020). The validation of target protein-ligand complex structures was performed using the co-crystallized standard ligand of target proteins to ensure the virtual screening process. For this, Mpro was redocked with its co-crystalized ligand RGZ: methyl 4-sulfamoylbenzoate. For each ligand-protein complex 10-docking poses were generated using Lamarckian Genetic Algorithm (López-Camacho et al., 2015). The pose with the lowest binding energy was selected as the final docking result. The interactions of the target-ligand complex were analysed and rendered using Discovery Studio v.20.1.0.19295.

Molecular fractionation with conjugate caps (MFCC—ab initio)

To calculate the interaction energy between receptors and ligands, the complexes obtained by molecular docking were used as an input for the molecular fractionation method with conjugated caps (MFCC). For this purpose, a convergence radius of 15 Å was defined (having its origin in the geometric center of the ligands). All the residues with at least one atom inside the spherical volume were taken into account. All atomic positions were kept fixed, with the exception of hydrogen atoms, which were optimized using consistent valence force field (CVFF). Subsequently, simulations within the Density Functional Theory formalism using the Local Density Approximation (DFT/LDA) were carried out using the DMOL3 code.

To obtain the interaction energy between each eligible amino acid and the ligands, a fractionation was applied, respecting Eq. (1): where represents the interaction energy between the ligand L and the amino acid Ri, where cap ans is arranged including the predecessor amino acid and posteriorly respectively of The term is the total energy of the system consisting of the ligand and the residue with its caps; the is the total energy of the waste and the isolated caps the total energy of the system formed by the ligand and only the caps; and finally, the energy of the system composed by the caps. The total energy of the ligand interaction with the entire receptor is obtained by adding the individual energy of each eligible amino acid, as in Eq. (2). where N is the number of eligible amino acids.

Molecular dynamics and binding energy calculations

Molecular dynamics and simulation (MDS) studies were carried out in order to determine the stability and convergence of E38 and E40 molecules complexed with Mpro. To set up the simulation, initially, the systems were built for complexes of E38-Mpro and E40-Mpro, respectively, in a system builder. For this purpose, Desmond 2018-4 was used to set up the initial parameters within an explicit SPC water model and placed in the orthorhombic box 4.0 × 4.0 × 4.0 Å. All the protease-ligand complexes were neutralized with NaCl by adding 0.15 M Na+ ions. ASL module was used to select the specific residues of ligand and protein molecule. The prepared systems were relaxed using Desmond default protocol of relaxation (Khan et al., 2020). MDS run of 20 ns was set up at constant temperature and constant pressure (NPT) for the final production run. The NPT ensemble was set up using Nosé–Hoover chain coupling scheme (Dayer et al., 2017) at temperature 300 K for final production and throughout the dynamics with relaxation time 1 ps. RESPA integrator was used to calculate the binding interactions for a time step 2 fs (Kaiser, 2005). All the other parameters were associated in the settings as described elsewhere (Chowdhury, 2020). After the final production run, the simulation trajectories of main protease complexed with E38 and E40 molecules were analyzed for final outcome of RMSD, RMSF and number of hydrogen bonds formed from the simulation. Binding energies of the complexes were calculated using MM-GBSA (Genheden & Ryde, 2015) for every 1 ns trajectory up to 20 ns and the average binding energies with standard deviations were measured for accurate binding approximation and stability described elsewhere (Genheden & Ryde, 2015).

Results and discussion

In silico pharmacokinetic and ADMET assays

In silico screening of compounds against new coronavirus targets is considered a useful strategy for the detection of molecules with great probability of presenting good pharmacokinetic properties and in vitro and in vivo low toxicity (Yang et al., 2020). Thus, we used in silico approaches to predict the pharmacokinetics and the toxicological potential of 101 phthalimide-1,2,3-triazole derivatives, aiming to identify leader compounds against SARS-CoV-2 (Figure 1). These compounds were categorized according to the radicals attached to the 1,2,3-triazole-phthalimide nucleus into groups A-H, as showed in supplementary material (Tables S1–S8): (A) derivatives containing benzyl substituents (Tehrani et al., 2019); (B) derivatives containing phenyl or benzothiazole groups (Da Silva et al., 2019); (C) derivatives containing carbohydrate groups and an additional phthalimide (Assis et al., 2012); (D) compounds with N-phenylacetamides substituents (Phatak et al., 2019); (E) presence of benzimidazole group between phtalimide-triazole nucleus and three radicals variations (Singh et al., 2020); (F) thalidomide derivative analogues (Ronnebaum & Luzzio, 2016); (G) derivatives of phthalimide-benzamide-1,2,3-triazole (Sadat-Ebrahimi et al., 2020); and (H) derivatives with an additional phthalimide (López-González et al., 2016).

Figure 1.

The chemical structures of 1,2,3-triazole-phthalimide bioactive derivatives obtained from SciFinder database. R1: 4-methylbenzyl 2-methylbenzyl, 4-methoxybenzyl, 2-fluorobenzyl, 3-fluorobenzyl, 4-fluorobenzyl, 2-chlorobenzyl, 3-chlorobenzyl, 2,3-dichlorobenzyl, 3,4-dichlorobenzyl, 2-bromobenzyl, 4-bromobenzyl, 2-nitrobenzyl; R2: 2-(methylsulfanyl)-1,3-benzothiazole or phenyl; R3: carbohydrate or phthalimide; R4: 4-methyl, 4-chloro, H, 4-bromo, 3-nitro, 4-methoxy, 2-chloro, 3-chloro, 4-fluoro, 2-methoxy, 3-methoxy, 2-methyl; R5: Methyl or H; R6: methylbutyl, methyl and 2-methylpropyl; R7: 4-fluorophenyl, p-tolyl, H, 1-pentafluorophenyl, 4-(trifluoromethyl)phenyl, 1-benzyl, 1-phenyl, 4-methoxyphenyl, 2,4-difluorophenyl, 2,6-difluorophenyl, 1-perfluorophenyl, 2-fluorophenyl, 4-pyridinyl; R8: H, 4-methoxybenzyl, 2-fluorobenzyl, 4-fluorobenzyl, 4-chlorobenzyl, 2,6-dichlorobenzyl, 3-bromobenzyl, 4-bromobenzyl, 2-nitrobenzyl, 3-nitrobenzyl, 4-nitrobenzyl.

For hit compounds selection, molecules from the above-mentioned series, that violated the Lipinski’s rule of five (RO5) and therefore, those which have poor druglike properties, were filtered out. The lipophilicity and solubility are the key molecular properties for drug absorption. Lipinski's rule of five is a rule of thumb that describes the druggability of a given molecule. Our pharmacokinetic analysis (supplementary material, Tables S9 and S10) revealed that 44 out 101 compounds violated at least one of the Lipinski’s rules. Out of these, 11 derivatives were predicted to have low gastrointestinal absorption (GIA). Although the compounds B03 and B04 did not violate the RO5, they showed a low GIA profile and were also rejected as hit compounds. Thus, 57 compounds were further submitted to ADMET analysis (Table 1).

Table 1.

Analysis of physicochemical and pharmacokinetic profiles of selected phthalimide-1,2,3-triazole derivatives.

	Physicochemical properties					Absorption		Metabolism					Excretion
	HBA	HBD	CLogP	MW	NRB	BBB	PGP	1A2	2C19	2C9	2D6	3A4	TC
A01	4	0	2.47	332.36	4	Yes	No	Yes	Yes	Yes	No	No	−0.07
A02	4	0	2.43	332.36	4	Yes	No	Yes	Yes	Yes	No	No	−0.06
A03	5	0	2.10	348.36	5	No	No	Yes	Yes	Yes	No	Yes	−0.04
A04	5	0	2.42	336.32	4	Yes	No	Yes	Yes	Yes	No	No	0.15
A05	5	0	2.42	336.32	4	Yes	No	Yes	Yes	Yes	No	No	0.12
A06	5	0	2.44	336.32	4	Yes	No	Yes	Yes	Yes	No	No	0.10
A07	4	0	2.63	352.77	4	Yes	No	Yes	Yes	Yes	No	No	0.34
A08	4	0	2.64	352.77	4	Yes	No	Yes	Yes	Yes	No	No	0.28
A09	4	0	3.15	387.22	4	Yes	No	Yes	Yes	Yes	No	No	0.46
A10	4	0	3.17	387.22	4	Yes	No	Yes	Yes	Yes	No	No	0.40
A11	4	0	2.71	397.23	4	Yes	No	Yes	Yes	Yes	No	No	0.32
A12	4	0	2.75	397.23	4	Yes	No	Yes	Yes	Yes	No	No	0.25
A13	6	0	1.49	363.33	5	No	No	No	Yes	No	No	No	−0.07
B01	5	0	3.03	407.47	5	No	No	Yes	Yes	Yes	No	Yes	0.41
B02	5	0	3.19	421.50	6	No	No	Yes	Yes	Yes	No	Yes	0.49
B05	4	0	2.04	304.30	3	Yes	No	Yes	Yes	No	No	No	0.11
B06	4	0	2.23	318.33	4	Yes	No	Yes	Yes	Yes	No	No	0.11
B07	4	0	2.54	332.36	5	Yes	No	Yes	Yes	Yes	No	No	0.19
B08	4	0	2.87	346.38	6	Yes	No	Yes	Yes	Yes	No	Yes	0.15
C01	10	0	1.03	484.46	11	No	No	No	No	No	No	Yes	1.37
C04	6	0	1.46	401.37	5	No	Yes	No	Yes	Yes	No	Yes	−0.23
C05	6	0	1.78	415.40	6	No	Yes	No	Yes	Yes	No	Yes	−0.14
C06	6	0	2.09	429.43	7	No	Yes	No	Yes	Yes	No	Yes	−0.13
D01	5	1	1.75	375.38	6	No	Yes	No	No	Yes	No	Yes	−0.07
D02	5	1	1.92	395.80	6	No	No	No	Yes	Yes	No	Yes	0.05
D03	5	1	1.41	361.35	6	No	Yes	No	No	No	No	No	−0.07
D04	5	1	2.06	440.25	6	No	No	No	Yes	Yes	No	Yes	0.03
D06	6	1	1.45	391.38	7	No	Yes	No	No	Yes	No	Yes	−0.05
D07	5	1	2.00	395.80	6	No	No	No	Yes	Yes	No	Yes	0.32
D08	5	1	1.94	395.80	6	No	No	No	Yes	Yes	No	Yes	0.11
D09	6	1	1.72	379.34	6	No	No	No	No	No	No	No	−0.24
D10	6	1	1.46	391.38	7	No	Yes	No	No	Yes	No	Yes	−0.11
D11	6	1	1.42	391.38	7	No	Yes	No	No	Yes	No	Yes	0.02
D12	5	1	1.75	375.38	6	No	Yes	No	No	Yes	No	Yes	−0.06
E03	5	1	3.37	430.50	6	No	No	No	Yes	Yes	Yes	Yes	0.51
E21	6	0	3.31	491.54	7	No	Yes	No	Yes	Yes	No	Yes	0.73
E25	6	0	3.47	466.47	5	No	Yes	No	Yes	Yes	No	No	0.54
E26	5	0	3.53	462.50	5	No	Yes	No	Yes	Yes	No	No	0.72
E27	5	0	3.22	448.48	5	No	Yes	No	Yes	Yes	No	No	0.71
E29	6	0	3.58	466.47	5	No	Yes	No	Yes	Yes	No	No	0.56
E31	5	0	3.29	462.50	5	No	Yes	No	Yes	Yes	No	Yes	0.81
E33	5	1	2.17	374.40	4	No	No	No	Yes	Yes	No	Yes	0.84
E35	5	0	3.90	476.53	5	No	Yes	No	Yes	Yes	No	Yes	0.72
E37	5	1	2.54	388.42	4	No	No	No	Yes	Yes	No	Yes	0.85
E38	5	0	3.62	462.50	5	No	Yes	No	Yes	Yes	No	No	0.72
E39	5	0	3.68	476.53	6	No	Yes	No	Yes	Yes	No	Yes	0.81
E40	6	0	2.89	463.49	5	No	Yes	No	Yes	Yes	No	Yes	0.76
E45	6	0	3.99	494.52	6	No	Yes	No	Yes	Yes	No	No	0.44
F01	6	1	0.33	325.28	2	No	No	No	No	No	No	No	−0.17
F02	6	1	0.24	339.31	3	No	No	No	No	No	No	No	−0.19
G01	5	1	2.87	437.45	7	No	Yes	No	Yes	Yes	No	Yes	0.05
G02	6	1	2.90	467.48	8	No	Yes	No	Yes	Yes	Yes	Yes	0.04
G03	6	1	3.19	455.44	7	No	Yes	No	Yes	Yes	No	Yes	−0.15
G04	6	1	3.21	455.44	7	No	Yes	No	Yes	Yes	No	Yes	−0.20
G05	5	1	3.36	471.90	7	No	Yes	No	Yes	Yes	No	Yes	−0.20
H01	6	0	1.39	387.35	4	No	Yes	No	Yes	No	No	No	−0.25
H02	6	0	1.33	387.35	4	No	Yes	No	Yes	No	No	No	−0.26

HBA: Hydrogen-Bond Acceptors, HBD: Hydrogen-Bond Donors, CLogP: Octanol-water partition coefficient, MW: Molecular Weight, NRB: Number Rotatable Bonds, BBB: Blood-Brain Barrier, PGP: P-glycoprotein Subtract, TC: Total Clearance (Log mL/min/kg).

The high oral bioavailability, preconized by Lipinski’s rules, is one of main desirable characteristics of a druggable compound and an important factor for the optimization of bioactive molecules as therapeutic agents. Oral route is one of the preferred pathways of drug administration due to its unique advantages, including sustained and controllable drug delivery, easy administration and high patient compliance (Ekins et al., 2010; Homayun et al., 2019). Passive intestinal absorption (associated with low MW), reduced molecular flexibility (measured by NRB), low TPSA or total hydrogen bond counts (HBA and HBD) are important predictors of good oral bioavailability (Brito, 2011). A compound having CLogP ranging 0.5 to 3.5 and molecular mass < 500 kDa is considered a druggable candidate (Tetko et al., 2016).

The presence of more than two heterocycles groups in the same molecule can result in compounds with high molecular weight and LogP values. These characteristics can be observed in most of the compounds belonging to the group E (1,2,3-triazole-benzoimidazole-phthalimide derivatives) and in compounds G07 and G08 (1,2,3-triazol-benzamide-phthalimide), which are therefore excluded from our study. Independently of molecular weight, compounds with high probability of good oral bioavailability have no more than 5 hydrogen bond donors and no more than 10 hydrogen bond acceptors. According to these latter criteria, our predictive analysis showed that the C03 and C02, D05, G09, G10, and G11 have low permeability throughout biological barriers.

The blood–brain barrier (BBB) protects the central nervous system (CNS) preventing certain substances (mostly harmful) from entering the brain tissues. The BBB limits the passage of most of the external compounds to maintain CNS at a steady state (Durán-Iturbide et al., 2020; Miao et al., 2019). In addition, P-glycoprotein (P-gp) actively transports a wide variety of compounds out of cells. This protein is highly associated with the ADMET properties, playing a major role in the multidrug resistance (MDR) phenomenon (Li et al., 2014). Our ADMET analysis showed that 15 out 57 compounds are predicted to cross the BBB, and 27 compounds were identified as a potential substrate of P-gp. Total Clearance (TC) is an important parameter associated with both the half-life and bioavailability of a xenobiotic compound. This parameter has a direct impact in determining the dose regimen (how often) and dose size (how much) of a given drug. Thus, its prediction helps us to determine the feasibility of clinical dosing and provides a framework for the starting dose in in vivo studies (Durán-Iturbide et al., 2020). In this work, the TC values of 1,2,3-triazoles-phtalimide derivatives varied from −0.26 to 0.85 log (mL/min/kg), for H2 and E37 derivatives, respectively. Higher values of TC suggest that the xenobiotic is removed rapidly from the body, whereas a low clearance value indicates slower removal.

Prediction of toxicological potential of 1,2,3-triazole-phtalimide compounds

Poor pharmacokinetics and adverse side effects, due to compound toxicity, are pointed as the main causes of late-stage failures in drug development (Li et al., 2019). In this regard, in silico predicted toxicological potential of the above-mentioned compounds was investigated. Our results showed that 15 out 57 compounds were positive for AMES and 46 compounds presented potential to be hepatotoxic. Of these, 10 compounds were both positive for AMES and hepatotoxic. The analysis of 8 selected compounds, that were either negative for AMES test or presented no hepatotoxic potential, revealed a predicted LD50 for acute toxicity in rats ranging from 1,888 to 2,980 mol/kg. The minimum value for oral chronic toxicity (ORCT) was 0.160 and the maximum value was 2,046 mg/kg/day. Besides presenting positive score for AMES test and hepatotoxicity the compound A09 also showed potential to cause adverse effects on the skin (Table 2).

Table 2.

In silico toxicological prediction of 1,2,3-triazole-phthalimide derivatives that showed better pharmacokinetic profile.

Comp.	AMES toxicity	MTD	hERGI	hERG II	ORAT (LD50)	ORCT	Hep	IT	SS
A01	No	−0.43	No	Yes	1.984	1.317	Yes	No	No
A02	No	−0.40	No	Yes	1.994	1.306	Yes	No	No
A03	Yes	−0.30	No	Yes	2.112	1.373	Yes	No	No
A04	No	−0.13	No	Yes	1.991	1.823	Yes	No	No
A05	Yes	−0.22	No	Yes	2.062	1.404	Yes	No	No
A06	No	−0.31	No	Yes	2.050	1.384	Yes	No	No
A07	No	−0.39	No	Yes	2.044	1.236	No	No	No
A08	No	−0.40	No	Yes	2.035	1.269	No	No	No
A09	No	−0.40	No	Yes	2.178	1.225	Yes	No	Yes
A10	No	−0.40	No	Yes	2.156	1.257	Yes	No	No
A11	No	−0.39	No	Yes	2.053	1.226	Yes	No	No
A12	No	−0.42	No	Yes	2.040	1.236	Yes	No	No
A13	Yes	−0.38	No	Yes	2.980	1.595	Yes	Yes	No
B01	No	0.38	No	Yes	2.309	0.850	Yes	No	No
B02	No	0.40	No	Yes	2.244	0.924	Yes	No	No
B05	No	0.48	No	Yes	2.541	1.133	No	No	No
B06	No	0.47	No	Yes	2.429	1.298	Yes	No	No
B07	No	0.51	No	Yes	2.460	1.411	Yes	No	No
B08	No	0.55	No	Yes	1.888	1.171	Yes	No	No
C01	No	0.63	No	No	2.439	1.309	Yes	No	No
C04	No	−0.14	No	Yes	2.272	2.020	Yes	No	No
C05	No	−0.13	No	Yes	2.243	2.046	Yes	No	No
C06	No	−0.09	No	Yes	2.221	1.992	Yes	No	No
D01	No	−0.21	No	Yes	2.153	1.674	Yes	No	No
D02	No	−0.18	No	Yes	2.148	1.631	Yes	No	No
D03	No	−0.25	No	Yes	2.109	1.696	Yes	No	No
D04	No	−0.19	No	Yes	2.158	1.624	Yes	No	No
D06	No	−0.14	No	Yes	2.080	1.654	Yes	No	No
D07	No	−0.25	No	Yes	2.151	1.625	Yes	No	No
D08	No	−0.22	No	Yes	2.171	1.693	Yes	No	No
D09	No	−0.12	No	Yes	2.065	1.605	Yes	No	No
D10	No	−0.25	No	Yes	2.093	1.648	Yes	No	No
D11	No	−0.21	No	Yes	2.111	1.691	Yes	No	No
D12	No	−0.28	No	Yes	2.155	1.668	Yes	No	No
E03	No	0.24	No	Yes	2.253	2.644	Yes	No	No
E21	No	0.23	No	Yes	2.480	1.690	Yes	No	No
E25	No	0.12	No	Yes	2.413	0.673	Yes	No	No
E26	Yes	0.08	No	Yes	2.416	0.803	No	No	No
E27	Yes	0.117	No	Yes	2.417	0.721	No	No	No
E29	Yes	0.121	No	Yes	2.416	0.555	Yes	No	No
E31	Yes	−0.028	No	Yes	2.435	0.822	No	No	No
E33	Yes	0.042	No	Yes	2.191	1.663	No	No	No
E35	No	0.200	No	Yes	2.463	1.074	No	No	No
E37	Yes	0.095	No	Yes	2.115	1.868	No	No	No
E38	No	0.217	No	Yes	2.465	0.993	No	No	No
E39	No	0.108	No	Yes	2.459	1.115	No	No	No
E40	No	0.225	No	Yes	2.441	1.590	No	No	No
E45	No	0.276	No	Yes	2.484	0.917	Yes	No	No
F01	No	−0.151	No	No	2.482	1.831	Yes	No	No
F02	No	0.095	No	No	2.364	1.965	Yes	No	No
G01	Yes	0.277	No	Yes	2.647	0.584	Yes	No	No
G02	Yes	0.236	No	Yes	2.853	0.578	Yes	No	No
G03	Yes	0.378	No	Yes	2.776	0.506	Yes	No	No
G04	Yes	0.292	No	Yes	2.816	0.532	Yes	No	No
G05	Yes	0.166	No	Yes	2.710	0.395	Yes	No	No
H01	No	0.015	No	Yes	2.110	1.889	Yes	No	No
H02	No	−0.018	No	Yes	2.216	1.985	Yes	No	No

MTD: Maximum Tolerated Dose—log(mg/kg/day), hERH I/II: Human Ether-a-go-go-related Gene, ORAT: Acute Oral Toxicity in Rats (LD50)—mol/kg, ORCT: oral rat chronic toxicity log(mg/kg/day), Hep: Hepatotoxicity, IT: Immunotoxicity; SS: Skin Sensitization.

In silico predictive toxicological evaluation is easy to perform, cost-effective, high-throughout alternative to in vitro assays and, even replaces time-consuming, expensive in vivo experiments (Myatt et al., 2018). Derivatives E26, E27, E31, E33, and E37 showed a positive predictive result for the AMES test, a method that assesses the mutagenic potential of chemical compounds (Kauffmann et al., 2020). The absence of a methyl at the carbon 5 of phthalimide ring and the presence of 1H-benzo[d]imidazole group may be related to the toxic effect of these compounds. In addition, changes in benzimidazole ligands can modulate their binding to deoxyribonucleic acid and the selectivity of these compounds towards the host genetic material (Gümüş et al., 2009).

Besides the mutagenic potential, another important variable taken in account in our predictive analysis was the evaluation of the hepatotoxic potential of selected phthalimide-1,2,3-triazole derivatives. Although the compounds A01, A02, A04, A06, A09, A10, A11 and A12, that have phthalimide-1,2,3-triazol-phenyl as the main structural group, were negative for AMES test, they showed important hepatotoxic potential. The presence of the methyl-benzyl group linked to 1,2,3-triazole-phthalimide nucleus, as well as the insertion of binders such as fluorine and bromine, may be related to the deleterious potential of these compounds. Interestingly, the compounds from the same series that did not show a predictive hepatotoxic effect, contained chlorine at the carbon 2 or 3 of phenyl group. However, this same result was not observed when 2,3-dichloro or 3,4-dichloro was present. Like the derivatives from the A series, the molecules E25, E45, G05, and G07 that also presented predictive hepatotoxic potential, have substituted halogenated groups in their structures. In order to better understand the structure-toxicity relationship of the compounds containing halogen atoms, further evaluation of other chemical-structural parameters is required (Kortagere et al., 2008).

The derivatives B01, B02, B06, B07, and B08 also showed positive predictive results for hepatotoxicity. The difference in the number of carbons in the aliphatic chain between phthalimide and triazole, also present in the molecules C04, C05 and C06, H01 and H02, could explain, at least in part, the predicted hepatotoxic potential found for compounds belonging to this series. Consistently, the number of carbons present in the aliphatic chain, that joins pharmacophoric groups in hybrid derivatives of 1,2,3-triazole-phthalimide, has been associated with the occurrence in vivo hepatotoxicity (Da Silva et al., 2019). The presence of the benzothiazole group in compounds of the E series did not appear to be determinant for the predicted hepatotoxicity of these molecules, since other triazole-phthalimide derivatives containing the same group, like E26, E27, E38, E31 and E33 did not present hepatotoxic potential. Heterocyclic molecules containing the benzothiazole group have shown to be more promising drug candidates to treat infectious diseases, reinforcing the pharmacological potential of benzothiazole-containing molecules (Papadopoulou et al., 2013). It has been reported that molecules containing thalidomide group in its structure, as those of series F, presented high hepatotoxic potential (Kamiya et al., 2020). In this regard, only the derivatives A7, A8, B5, E35, E38, E39, and E40 (Figure 2) showed good pharmacokinetic and toxicological profiles, and were therefore chosen for molecular docking.

Figure 2.

The chemical structures of derivatives A7, A8, B5, E35, E38, E39, and E40 chosen for molecular docking studies.

Molecular docking of selected 1,2,3-triazole-phtalimide compounds

Computational simulations of ligand-protein docking are an important component of drug discovery, being mainly used for the virtual screening of hit compounds from large databases and to identify and evaluate the effects of chemical changes during leader optimization (Bordogna et al., 2011). In this regard, the accuracy of modeled structures should be taken into account. Our stereochemical analysis of SAR-CoV-2 target proteins, obtained from PDB (Figures S1–S5), showed that the percentage of the sum of amino acid residues in the most favoured and additional allowed regions was ≥ 95%. This result confirmed the high quality of our selected target models (Laskowski et al., 2006). In order to validate our molecular docking, the co-crystal standard ligand (RGZ, methyl 4-sulfamoylbenozate) of Mpro complex (5R80) was re-docked against this protein and the root-mean-square deviation (RMSD) was calculated for predicting the stability of the protein and protein-ligand complexes. Our results showed that Autodock program was able to generate a successful pose (RMSD = 1.9 Å) (Figure S6). It is assumed that scored poses with an RMSD of less 2.0 Å are considered to be successful (Bordogna et al., 2011; Razzaghi-Asl et al., 2020).

Figure 6.

Energy for each nucleocapsid phosphoprotein amino acids in a radius of 15 Å from E40 ligand.

The docking analysis using iGemDock software revealed that the derivatives A7, A8, B5, E35, E38, E39, and E40 showed favourable interactions energies for all the tested proteins from SARS-CoV-2 (Table 3). E38 presented higher binding energy valued for spike protein in its closed state (Spike-C) (−116.6 kcal/mol) and spike protein in prefusion conformation (−107.1 kcal/mol). The higher binding energy value for main protease (Mpro) was found for its interaction with E38 (−117.2 kcal/mol). Among the tested compounds, E39 presented highest affinity to RNA binding domain of nucleocapsid phosphoprotein with binding energy of −108.8 kcal/mol. For all ligand-protein complexes, the interaction energies were predominantly due to Van der Waals, followed by hydrogen bonds. No electrostatic interactions were observed for all interactions of compounds with SARS-CoV-2 proteins.

Table 3.

iGemDock predictive profile of binding energies for A07, A08, B05, B35, E38, E39 and E40 triazole-phthalimide derivatives against SARS-CoV-2 targets.

		Mpro	Spike protein			NCP
Comp.	Energy (kcal/mol)	5R80	6VSB	6VXX	6VYB	6VYO
A07	E	−103.8	−94.83	−92	−89.47	−94.5
	VDW	−95.41	−71.12	−79.44	−74.98	−80.5
	HB	−8.4	−23.71	−12.56	−14.5	−14
	Elec	0	0	0	0	0
A08	E	−103.4	−93.35	−91.74	−101.5	−94.5
	VDW	−95.32	−83.81	−84.52	−86.46	−80.5
	HB	−8.1	−9.54	−7.22	−15.1	−14
	Elec	0	0	0	0	0
B05	E	−89.9	−99.48	−107.5	−84.81	−88.18
	VDW	−79.51	−66.62	−76.07	−67.26	−75.4
	HB	−10.4	−32.86	−31.48	−17.55	−12.7
	Elec	0	0	0	0	0
E35	E	−107.5	−98.8	−102.2	−109.1	−102.5
	VDW	−95.5	−78.2	−92.48	−91.94	−92.9
	HB	−12	−26.59	−9.80	−17.24	−9.6
	Elec	0	0	0	0	0
E38	E	−117.2	−107.1	−116.6	−93.2	−105.3
	VDW	−109.1	−91.53	−94.08	−80.36	−87.3
	HB	−8.1	−15.61	−22.61	−12.84	−18
	Elec	0	0	0	0	0
E39	E	−103.3	−99.32	−108.4	−98.55	−108.8
	VDW	−99.3	−83.36	−91.46	−90.85	−91.8
	HB	−4	−15.96	−16.94	−7.7	−17
	Elec	0	0	0	0	0
E40	E	−106.3	−105.1	−94.82	−96.75	−101.9
	VDW	−99.36	−82.52	−91.02	−87.53	−89.4
	HB	−7	−22.65	−3.8	−9.21	−12.5
	Elec	0	0	0	0	0

E: energy; VDW: van der Waals; HB: hydrogen-bond; Elec: electrostatic interactions.

The use of molecular docking has been largely applied as a strategy for the selection of possible leader candidates against the new coronavirus (Karypidou et al., 2018). Since the SARS-CoV outbreak in 2002, in silico molecular docking and structural studies have demonstrated a relationship between the SARS-CoV-2 spike receptor-binding domain and the angiotensin-2 receptor (ACE2) in the host cells. This association is strictly linked to the infection and transmission of new corona virus in humans (Utomo & Meiyanto, 2020). The wide distribution of the ACE2 receptor in human tissues and organs could be related to the appearance of symptoms in the severe phase of the disease, mainly in the pulmonary system, where it is largely expressed at the surface of type II alveolar cells (Sanchis-Gomar et al., 2020). In our analysis we docked the selected compounds on different conformational states of SARS-CoV-2 spike protein. Overall, our results showed that all the selected compounds bound with high affinity to spike proteins independently of its conformational state (open, close, or in prefusion conformation). However, slight differences in binding energy values were observed. The lowest binding energy value found for the interaction of E38, and E35 derivatives with different conformational states of spike protein, revealed the binding plasticity of 1,2,3-triazole-phthalmide derivatives towards this protein. Liu et al. (2004) described that the inhibition of peptides in the heptad repeat regions of the SARS-CoV spike protein results in the interference of the fusogenic mechanism. In this regard, it is possible that the selected Top-7 derivatives may interfere with the SARS-CoV-2 binding, entrance and replication processes in the host cells.

The predictive analysis of the inhibition constant values for best binding pose of protein-ligand complexes was performed using AutoDock software (Tables S11–S17). Our results showed that A7, E35, E38 and E40 presented higher binding affinity for at least one of the viral target SARS-CoV-2 proteins, with inhibition constant values (Ki) in the nanomolar range. All of these compounds showed high affinity for SARS-CoV-2 main protease. Comparing all ligand-protein interactions tested, the lowest values of inhibition constant (Ki) were found for the interaction of the compound E38 with Mpro with Ki = 23.85 nM and binding energy of −10.4 kcal/mol for the best pose. The second-best value of Ki was obtained for the E40-Mpro complex with Ki = 30.04 nM and a binding energy of −10.26 kcal/mol. Compound E40 also bound with high affinity to nucleocapsid protein, with Ki = 91.09 nM and binding energy of −9.6 kcal/mol (Figure 3).

Figure 3.

the graphical representation of binding energy of E38 and E40 derivatives on SARS-CoV-2 proteins. Spike (in prefusion conformation) Spike-O: spike protein in its open state; Spike-C: spike protein in its closed state; NCP: nucleocapsid protein; Mpro: main protease.

In addition to the spike protein, coronavirus proteases are also essential for viral transmission and virulence by processing viral proteins involved in viral replication (Báez-Santos et al., 2015). Compounds that act as inhibitors of these proteins can be considered as promising antiviral agents by reducing the severity of the infection (Hall & Ji, 2020). All the top-7 compounds bound with high affinity to SARS-CoV-2 main protease, being E38 and E40 the compounds that presented the highest affinity (i.e. the lowest binding energy) towards this target. The analogous performance of these compounds toward Mpro can be attributed to chemical-structural similarities of the two derivatives, which are only distinguished from each other by the presence of the phenyl and pyridin-4-yl groups in the N1 of triazole ring. Molecules containing 1,2,3-triazole group have already been described as having anti-coronavirus activity, by interacting with viral proteases, corroborating the relevance of these proteins as the main target of these compounds (Karypidou et al., 2018). The main amino acids involved in E40-SARS-CoV-2 binding were MET165, MET49, HIS41, CYS145, GLU166 and GLN189 (Mpro—Figure 4) and ALA134, GLY69, VAL133, TRP132, ILE130, ALA125, ILE131, LYS65, PHE66, ARG68 and TYR123 (NCP—Figure 5). The most prevalent interaction was Van der Waals and hydrogen bonds. The electrostatic binding energy values were insignificant when compared to the contributions of van der Waals and hydrogen bonds.

Figure 4.

The molecular docking of derivative E40 with Mpro. (A) Representative diagram of docking result of E40 with Mpro. (B) Binding interaction of E40 with amino acid residues of Mpro.

Figure 5.

The molecular docking of derivative E40 with NCP. (A) Representative diagram of docking result of E40 with NCP. (B) Binding interaction of E40 with amino acid residues of NCP.

In this study we also docked the Top-7 1,2,3-triazole-phthalimide derivatives on viral nucleocapsid protein. NCPs have been considered strategic targets to fight acute respiratory infections caused by SARS-CoV-2 due to their multiple functions in the viral replication cycle (Kang et al., 2020). All the compounds tested bound with high efficiency on NCP. Comparing the best-pose of selected molecules (1/10) for the NCP, the E40 presented the most favorable binding energy (−9.6 kcal/mol) and was the only one that achieved an inhibition constant in the nanomolar range for this target. The interaction of these compounds with nucleocapsid proteins may interfere with essential functions of the viral cycle, such as the formation of the helical ribonucleoprotein complex, during the packaging of the viral genome, and the control of basic cellular functions in the host cell, like the deregulation of the cell cycle (Cong et al., 2020; Surjit et al., 2006). Although further in vitro analysis are needed to better understand the biological effects of the interaction of E40-SARS-CoV-2 nucleocapsid proteins complex, the high predicted binding affinity of this compound to these proteins suggests that E40 may prevent nucleocapsid proteins from participating in the replication and release of viral particles (Zeng et al., 2008).

A characteristic of RNA viruses is the high rate of genetic mutation, which can result in the evolution of new viral strains and the inefficacy of treatment. Therefore, from a public-health perspective, the ability of our top-7 selected 1,2,3-triazole-phthalimide derivatives to bind to different viral proteins, may circumvent the drug resistance caused by single-target mutations or rare simultaneous mutations of several targets in different positions (Xu et al., 2020). Although some compounds were excluded from our study, it is important to take in mind that these molecules can be redesigned for more rational chemical-structural modifications, decreasing its toxic potential and/or improving their biological activity against target molecules.

Because our molecular docking analysis has pointed the compound E40 as the most promising compound, the best complexes obtained between this compound and the targets were subjected to quantum energy calculations by the MFCC technique. The results of the interactions between E40 and amino acids of the RNA binding domain of nucleocapsid phosphoprotein (PDB-ID 6VYO) are described in Figure 6 and Table 4.

Table 4.

Interaction energies between RNA binding domain of nucleocapsid phosphoprotein residues and E40.

Amino acids	Energy
ASP63	5.00
LEU64	6.10
LYS65	−7.00
PHE66	−4.32
PRO67	0.00
ARG68	−5.89
GLY69	−8.89
GLN70	−4.43
GLY71	−8.12
VAL72	−0.52
ILE84	−0.33
GLY85	−33.74
TYR86	−31.00
TYR87	−24.50
LEU121	−28.65
PRO122	−24.55
TYR123	−21.63
GLY124	−14.20
ALA125	−1.22
ASN126	−1.31
GLY129	−2.01
ILE130	−2.14
ILE131	1.87
TRP132	2.01
VAL133	2.02
ALA134	4.31
THR135	3.22
GLU136	4.87
ALA138	1.51
THR166	6.90
Total energy	−1.87E + 2 kcal/mol−8.03 eV

Overall, 30 amino acids of NCP showed some kind of interaction with E40, 11 showing repulsion and 19 attraction with the ligand. The PRO67 residue presented energy value of 0.00. Consequently, it cannot be confirmed whether the interaction occurring with the ligand is repulsive or attractive. The general energy balance points to a stronger attraction (−8.03 eV) between the connection site of protein and the E40 compound. The residues GLY85, TYR86, TYR87, LEU121, PRO122, TYR123 and GLY124 presented very low energies (negative values), suggesting very strong attractions by these amino acids, which are probably the anchor residues. Figure 7 and Table 5 show MFCC results between E40 and the open state of SARS-CoV-2 spike ectodomain structure (PDB-ID 6VYB).

Figure 7.

Energy for each Spike-O amino acids in a radius of 15 Å from E40 ligand.

Table 5.

Interaction energies between Spike-O amino acid residues and E40.

Amino acids	Energy
LYS1038	2.87
PHE906	−3.51
ASN907	−4.85
GLY908	0.31
ILE909	0.11
GLY910	0.20
VAL911	0.02
THR912	0.17
GLN1036	−2.14
SER1037	−22.17
LYS1038	−31.89
VAL1040	−15.72
ASP1041	−11.59
CYS1043	−9.25
GLY1046	−5.78
TYR1047	−5.23
HIS1048	−4.17
LEU1049	−10.83
THR1066	−4.22
TYR1067	−4.82
VAL1068	−10.11
GLU1092	−7.05
GLY1093	−4.04
GLN1106	4.58
ARG1107	2.89
ASN1108	−5.13
GLY885	−5.01
TRP886	−3.54
THR887	−0.17
GLY889	−5.12
ALA890	−3.21
GLN901	−3.71
ALA903	−2.63
TYR904	−1.72
ARG905	−2.47
ASN907	−18.12
GLY908	−2.14
ILE909	5.21
GLY910	5.87
GLU1031	3.17
GLY1035	3.16
GLN1036	−2.51
SER1037	−8.10
LYS1038	17.54
HIS1048	−5.81
Total energy	−1.81E +02 kcal/mol−7.77 eV

A total of 45 residues of spike-O showed energetic interactions with the drug, 13 of which are interacting repulsively and 32 attractively. The residues VAL1040, ASP1041, LEU1049, VAL1068, ASN907, and LYS1038 demonstrated a strong attraction for the ligand, possibly acting as anchoring residues. Figure 8 and Table 6 describe the results of the MFCC for interactions between spike glycoprotein in closed state (PDB-ID 6VXX) and E40.

Figure 8.

Energy for each Spike-C amino acid in a radius of 15 Å from E40 ligand.

Table 6.

Interaction energies between Spike-C amino acid residues and E40.

Amino acid	Energy
ARG357	3.12
SER359	−2.18
ASN360	−4.97
THR393	0.52
ASN394	0.00
GLU516	0.44
LEU518	0.17
ALA520	−0.07
PRO521	−2.13
ALA522	−18.18
THR523	−25.10
VAL524	−14.56
PRO561	−11.20
PHE562	−8.56
GLN563	−5.77
GLN564	−5.01
LYS41	−4.87
ILE128	−8.78
TYR170	−2.99
TYR200	−3.25
PHE201	−9.11
LYS202	−7.21
ILE203	−2.56
TYR204	−1.11
SER205	2.10
GLU224	−4.10
PRO225	−4.55
LEU226	−2.67
VAL227	−1.25
ASP228	−1.02
LEU229	−3.68
PRO230	−3.53
Total energy	−1.52E +02 kcal/mol−6.54 eV

There were 32 amino acids of spike protein interacting with E40, 5 with positive energies (repulsion), 26 with negative energies (attraction) and one with 0.0 energy (ASN394). The amino acids that stood out for having a strong attraction were: ALA522, THR523, VAL524 and PRO561. Although the structures obtained at inputs spike-C and spike-O represented different conformational states of the same structural protein of SARS-CoV-2 (spike glycoprotein), the compound E40 was strongly attached to both conformations but in different places. The results obtained for spike glycoprotein with a single receptor-binding domain up (6VSB) are shown in Figure 9 and Table 7.

Figure 9.

Energy for each spike glycoprotein (with a single receptor-binding domain up) amino acids in a radius of 15 Å from E40 ligand.

Table 7.

Interaction energies between spike glycoprotein (with a single receptor-binding domain up) residues and E40.

Amino acids	Energy
GLN321	3.01
PRO322	−3.02
SER325	−5.13
VAL539	0.43
ASN540	0.01
PHE541	0.05
ASN542	0.18
GLY545	−0.78
THR547	−1.89
GLY548	−21.58
THR549	−25.31
GLY550	−17.44
MET740	−12.37
CYS743	−9.81
GLY744	−5.11
ASP745	−5.12
SER746	−4.01
THR747	−9.78
GLU748	−3.55
CYS749	−4.02
SER750	−10.24
ASN751	−7.03
VAL976	−3.89
LEU977	−1.85
ASN978	2.54
ASP979	−4.12
ILE980	−5.08
LEU981	−2.50
SER982	−0.84
PRO986	−1.22
ALA989	−4.10
ILE993	−2.56
Total energy	−1.66E +02 kcal/mol−7.14 eV

A total of 32 residues of spike glycoprotein (with a single receptor-binding domain up) residues, interacted with the E40 ligand, 6 of which presenting repulsion and 26 attraction interactions, with emphasis on the amino acids GLY548, THR549, GLY550, MET740 and SER750, which showed strong attraction for the compound. Figure 10 and Table 8 showed the MFCC results for E40 interactions with the main protease in complex.

Figure 10.

Energy for each main protease amino acids in a radius of 15 Å from E40 ligand.

Table 8.

Interaction energies between main protease amino acid residues and E40.

Amino acids	Energy
ALA173	−0.25
ALA191	0.69
ARG188	0.00
ARG40	0.41
ASN142	0.02
ASN28	−0.50
ASN51	0.29
ASP187	−8.47
ASP48	−3.15
CYS145	−50.42
CYS44	−43.68
CYS85	−42.31
GLN189	−7.90
GLN192	−5.54
GLU166	−8.19
GLY138	−1.28
GLY143	0.26
GLY146	0.01
GLY170	0.03
GLY174	0.31
HIS163	−1.89
HIS164	−3.78
HIS172	0.23
HIS41	−7.06
ILE43	−0.22
LEU141	−3.94
LEU167	−0.23
LEU27	−0.11
LEU50	0.00
MET162	−4.10
MET165	2.00
MET49	−5.80
PHE140	−0.24
PHE185	−0.64
PRO168	−3.45
PRO39	0.00
PRO52	−0.50
SER144	−1.71
SER147	0.47
SER46	−0.23
THR169	0.29
THR175	−0.43
THR190	−3.40
THR25	−0.06
THR26	0.56
THR45	0.82
VAL171	0.50
VAL186	1.28
VAL42	0.67
Total energy	−2.01E +02 kcal/mol−8.63 eV

There are 49 amino acids of spike glycoprotein (with a single receptor-binding domain up) residues interacting with E40, 17 with slightly repulsive interactions, 3 with 0.0 energy and 29 with negative energy values. The CYS145, CYS44 and CYS85 amino acids residues, showed surprising attraction to E40, with energy values of −42.31 kcal/mol, −43.68 kcal/mol and −50.42 kcal/mol, respectively. These cysteine amino acids have already been identified as anchor residues in other studies on SARS-CoV-2 main protease (Komatsu et al., 2020). Taking in account the general energy balance, this protein also showed the lowest energy value in this study (−8.63 eV).

In a viral pandemic with rapid dispersion of cases, drug replacement strategy is promising to fight the disease, because substances that are commercially available in the pharmaceutical industry have already undergone several toxicity and safety tests. In chemistry, when atoms are held together, the energy needed to separate them is called binding energy. In covalent bonds, where the electrons are shared between atoms, the energy needed to break these bonds can vary, but energies less than −1.0 eV (Tkatchenko & Scheffler, 2009) suggest covalent bond among atoms. E40 showed lower energy values in all analyses. Interestingly, but not surprisingly, E40 showed affinity with all the three spike conformers (6VYB, 6VXX and 6VSB). Although, in this study, the ligand interacted with all Spike conformers, the anchorage residues were different in each interaction. This raises some concerns: would the E40 function as an inhibitor of the Spike-ACE2 interaction in any of the three conformers? Is there a possibility of allosteric inhibition by the E40? Although our results suggest a drug repositioning for COVID-19, further in vivo and in vitro studies are needed to answer these questions.

MD simulation and MM-GBSA binding energy calculations

The final convergence and the stability of E38 and E40 bound with the main protease was assessed by MD simulation. The analysis of RMSD and RMSF plots showed that after 20 ns of convergence, E3, bound to the main protease, displayed a significant stability (Figure 11A and 11B). The interactions of this protein with E38 displayed much less RMSD differences (∼0.25 Å) showing its stable conformation (Figure 11A). The average of H-bonds formed between E38 and the main protease showed to be 1.0 throughout the simulation period for the stable conformation of both the main protease and E38 molecule (Figure 11C). The structural superimposition of the main protease-E38 complex before and after simulation displayed less significant changes in the overall conformation (Figure 11D). In addition, accommodation of E38 at the binding cavity of this protein showed small-altered orientation after simulation, which seems to be due to a change in the visible secondary structure at the E38-M binding site (Figure 12).

Figure 11.

The molecular dynamic (MD) simulations of main protease bound to E38. (a) RMSD plot of Cα-backbone (blue) and ligand (red) displayed stable convergence till 20 ns. (b) RMSF plot depicting fluctuations of each residual position of main protease amino acids. (c) Number of H-bonds formed during the entire course of simulation. (d) Superimposed initial and final frame structures of main protease complexed with E38 (surface view) molecule before and after simulation.

Figure 12:

Surface view of main protease before MD simulation (orange) displaying little modifications in the E38-M binding site (yellow ring) after simulation (cyan teal). Lower panel depicted the conformational differences between the structures of E38 in the binding site provided orientation that is more accurate for higher binding.

After 20 ns of convergence, the E40-M complex displayed significant stability as demonstrated by RMSD and RMSF plots (Figure 13A and 13B). The interactions of this protein with E40 displayed much less RMSD differences (∼0.50 Å) conferring its stable conformation (Figure 13A). The average of H-bonds formed between E40 and the main protease showed to be 1.5 throughout the simulation period, for the stable conformation of both the main protease and the E40 molecule (Figure 13C). Structural superimposition of the main protease-E40 complex before and after simulation displayed less significant changes in the overall conformation (Figure 13D). In addition, accommodation of E40 at the binding cavity of Mpro showed small-altered orientation after the simulation which appeared to be due to the change in the visible secondary structure at the E40-MPro binding site (Figure 14).

Figure 13.

Molecular dynamic (MD) simulations of main protease bound to E40, (a) RMSD plot of Cα-backbone (blue) and ligand (red) displayed stable convergence till 20 ns. (b) RMSF plot depicting fluctuations of each residual position of amino acids of main protease. (c) Number of H-bonds formed during the entire course of simulation. (d) Superimposed initial and final frame structures of main protease complexed with E40 (surface view) molecule before and after simulation.

Figure 14.

Surface view of main protease before MD simulation (green) displaying little modifications in the E40-M binding site (yellow ring) after simulation (cyan teal). Lower panel represented the conformational differences between the structures of E40 in the binding site provided orientation that is more accurate for higher binding.

The free energies of binding using MMGBSA were calculated for the E38 and E40- Mpro complexes of SARS-CoV-2. The properties of MMGBSA calculations are displayed in Table 9. The free energy avarage of binding for every 1 ns trajectory up to 20 ns, for the main protease- E38 complex displayed a dG = −63.47 kcal/mol, with standard deviation of 3.0 kcal/mol, coulombic force of −8.70, solvent accessibility of 20.20 and high ligand binding affinity of −1.81. On the other hand, the E40-main protease complex displayed a dG = −63.31 with standard deviation of 7.0 kcal/mol and high binding efficiency, solvent accessibility and significant coulombic as well covalent energies. All the properties calculated in the MD simulation point the significant role of E38 and E40 compounds as a potent inhibitor of the main protease of SARS-CoV-2.

Table 9.

Properties calculated after MMGBSA for E38 and E40 with main protease.

Molecule	dG bind (kcal/ mol)	dG coulombic force (kcal/ mol)	dG covalent energy (kcal/ mol)	dG solvent accessibility (Å3)	Ligand efficiency (kcal/ mol)
E38	−63.47 ± 3	−8.70	2.11	20.20	−1.81
E40	−63.31 ± 7	−15.22	1.098	25.07	−1.79

Conclusion

Taken together, our results showed that the seven selected 1,2,3-triazole-phthalimide derivatives showed an excellent predicted pharmacokinetic and ADME properties, which are directly related to the solubility, permeability and toxicological profiles of such compounds. The in silico docking results showed that these selected compounds are potential multi-target ligands of essential SARS-CoV-2 spike, protease and nucleocapsid proteins suggesting the usefulness of 1,2,3-triazole-phthalimide derivatives for the development of drugs addressed to these proteins. Because the E40 presented predicted inhibition constant values in the nanomolar concentration range for both main protease and nucleocapsid proteins, this compound was considered the most promising as multi-target agent against SARS-CoV-2. Our MD simulation corroborated E40 and E38 as potent inhibitors of Mpro of SARS-CoV-2 However, we cannot rule out the possibility that other compounds, such as E38 can be included in further investigations alone or in combination with E40.

Click here for additional data file.