key: cord-0920679-2yfzmxhe
authors: Ničkčović, Vanja P.; Nikolić, Gordana R.; Nedeljković, Biserka M.; Mitić, Nebojša; Danić, Snežana Filipović; Mitić, Jadranka; Marčetić, Zoran; Sokolović, Dušan; Veselinović, Aleksandar M.
title: In silico approach for the development of novel antiviral compounds based on SARS-COV-2 protease inhibition
date: 2022-04-03
journal: Chem Zvesti
DOI: 10.1007/s11696-022-02170-8
sha: 420fe1f176c90f082584a3777b2421e6ae950fc9
doc_id: 920679
cord_uid: 2yfzmxhe

The COVID-19 pandemic emerged in 2019, bringing with it the need for greater stores of effective antiviral drugs. This paper deals with the conformation-independent, QSAR model, developed by employing the Monte Carlo optimization method, as well as molecular graphs and the SMILES notation-based descriptors for the purpose of modeling the SARS-CoV-3CLpro enzyme inhibition. The main purpose was developing a reproducible model involving easy interpretation, utilized for a quick prediction of the inhibitory activity of SAR-CoV-3CLpro. The following statistical parameters were present in the best-developed QSAR model: (training set) R(2) = 0.9314, Q(2) = 0.9271; (test set) R(2) = 0.9243, Q(2) = 0.8986. Molecular fragments, defined as SMILES notation descriptors, that have a positive and negative impact on 3CLpro inhibition were identified on the basis of the results obtained for structural indicators, and were applied to the computer-aided design of five new compounds with (4-methoxyphenyl)[2-(methylsulfanyl)-6,7-dihydro-1H-[1,4]dioxino[2,3-f]benzimidazol-1-yl]methanone as a template molecule. Molecular docking studies were used to examine the potential inhibition effect of designed molecules on SARS-CoV-3CLpro enzyme inhibition and obtained results have high correlation with the QSAR modeling results. In addition, the interactions between the designed molecules and amino acids from the 3CLpro active site were determined, and the energies they yield were calculated. SUPPLEMENTARY INFORMATION: The online version contains supplementary material available at 10.1007/s11696-022-02170-8.

The period from the end of 2019 and the beginning of 2020 will be remembered as the onset of the most hazardous pandemic in modern history, i.e. the COVID-19 outbreak. COVID-19 is caused by coronavirus 2 (SARS-CoV-2), or by the positive-sense single-stranded RNA virus Guo et al. 2020; Rothan and Byrareddy 2020) . The Center for Disease Control and Prevention (CDC) first established the existence of this virus in Wuhan, China in the final days of 2019, and from then on, the virus has spread to other countries, posing a considerable threat to global health (Lai et al. 2020; Pedersen and Ho 2020) . The droplets produced when those infected exhale, cough or sneeze represent the main means of SARS-CoV-2 distribution. the crucial factor in spreading the virus is the ease with which it is transmitted, since breathing in the proximity of someone who has COVID-19, or touching a virus contaminated surface before making contact with one's mouth, eyes, or nose, constitute the main routs of infection Singhal 2020; Sohrabi et al. 2020) . Conversely, the chief danger of contracting COVID-19 lies in the possibility of an asymptomatic person spreading the virus, since 10% of the infections originate from people who do not exhibit any symptoms whatsoever, in spite of the fact that COVID-19 is the most infectious once someone develops symptoms related to COVID-19 (Kampf et al. 2020; Xiao et al. 2020) . The chief symptoms attributed to COVID-19 are the potential loss of smell, cough, or taste, as well as the occurrence of fever, or having trouble breathing (Shang et al. 2020; Wu et al. 2020) . Fortunately, to date, proper care and treatment usually lead to the recovery of most people infected with COVID-19 Sanders et al. 2020; Vellingiri et al. 2020) . The WHO (World Health Organization) has determined two groups which pose a greater risk of suffering from severe illness and even mortality caused by the SARS-CoV-2 infection-adults over 60 years of age, and people whose immune systems are compromised by various health issues , Kampf et al. 2020 . Since SARS-CoV-2 is a recently discovered virus, a lot of its characteristics are still unknown.

One promising approach for COVID-19 treatment is blocking the replication of the virus by inhibiting viral protease (Lin et al. 2005; Hoffmann et al. 2020; Jin et al. 2020) . One of the key steps in viral replication, which also includes SARS-CoV-2, is the processing of replicase polypeptides 1a and 1ab into functional proteins, which is a process determined with 3C-like protease (3CLpro) (Chen et al. 2005; Deng et al. 2014; Ramajayam et al. 2010; Zhavoronkov et al. 2020) . Two proteases, papain-like protease (PLpro) and main protease or 3-chymotrypsin-like protease (Mpro/3CLpro), are essential components for replication of SARS-CoV-2, whit main role related to cleaving the two polyproteins, i.e., PP1A and PP1AB into several functional components. In this process, N-terminal domain is split to viral precursor protein at the three sites by PLpro, while the C-terminal domain is split to precursor protein at the 11 sites by 3CLpro (Akaji and Konno 2020) . Targeting these enzymes had been used for the treatment of other pathogenic coronaviruses (i.e., MERS-CoV) (Amin and Jha 2020) and since the sequence homology of SARS-CoV-2 3CLpro is 96% structurally closer to the SARS-CoV 3CLpro targeting SARS-CoV-2 3CLpro as a potential therapeutic target is quite a feasible approach for anti-CoVID-19 drug development (Amin et al. 2021; Zhang et al 2020) .

Chemoinformatics and computational chemistry have developed methodologies which have proven to be of great importance when designing potential therapeutics. When it comes to drug development research, the most significant contributions made by these in silico methods are the search for new leading compounds, or the optimization of therapeutic activities (or pharmacokinetic properties) of the series of chemical compounds whose biological activities have already been determined (Ekins et al. 2007; Tabeshpour et al. 2018 ). Among the most used methodologies are the quantitative structure-activity relationship (QSAR) and molecular docking. The mathematical equation linking the biological activities of the studied molecules with their chemical characteristics, defined as molecular descriptors, represents one of the most widely used practices for QSAR model representation (Cherkasov et al. 2014) . In current practice, the applied representation of molecule structure is used for the calculation of a considerable number of descriptors, calculated in such fashion that could be further utilized for the development of a relevant QSAR model (Liu and Long 2009; Pérez González et al. 2008) . A novel approach involves the use of conformation-independent optimal descriptors which are based on the constitutional and topological molecular features, and the descriptors based on the Simplified Molecular Input Line Entry System (SMILES) notation for QSAR model development, where the Monte Carlo optimization approach is used for this purpose. One of the recently suggested approaches to the manner of overcoming issues related to mechanistic interpretation is the application of SMILES descriptors, since these have a physical meaning and can be associated with molecular fragments (Toropova et al. 2016; Veselinović et al. 2015; Zivkovic et al. 2020) , which is why this type of QSAR modeling is conducted under the OECD (Organization for Economic Cooperation and Development) guidelines for QSAR model development.

Several in silico methods were administered for the purpose of revealing new compounds having a potential 3C-like protease (3CLpro) inhibition activity in the presented research. Since molecular descriptors based on the SMILES notation were used, along with local graph invariants, the developed QSAR models were conformation-independent. One of the chief aims of this study was to determine the structural requirements or molecular fragments responsible for the 3CLpro inhibition effect. Moreover, this research defined the structural attributes found in small molecules concerning ligand-receptor interactions, which could be used in designing and developing SARS-CoV infection treatment therapeutics. Since the SMILES notation is used for chemical representation by most of the currently utilized chemical databases, this type of modeling could be employed by other scientists for screening purposes, as the results of the developed QSAR models have displayed considerably good values for various validation metrics, both internal and external ones.

A group of 84 molecules encompassing a wide range of heterocyclic compound classes, with an established 3C-like protease inhibiting activity against the infectious bronchitis virus (SARS-CoV virus), were selected from the binding database, according to "The Scripps Research Institute Molecular Screening Center" (Gilson et al. 2015) , and used for QSAR model development. What needs to be specially emphasized is that the 3CLpro protease of SARS-CoV-2 has a sequence similarity greater than 95% when compared with that of the bronchitis virus (SARS-CoV virus), which was the main reason for the use of the surrogate enzyme for developing inhibitors against SAR-CoV-2 3CLpro . All the used molecules have an activity in IC 50 (nM) values, which was determined using the same bioassay protocol (QFRET-based dose-response biochemical high throughput screening assay) (Jacobs et al. 2013) . The IC 50 values were converted to pIC 50 (pIC 50 = − logIC 50 ) values as a studied activity. Afterwards, Open Babel was utilized to canonize the SMILES notations obtained from the abovestated database. Table S1 in the Supplementary Material cites the SMILES notations of all the used molecules, with their pIC50. Furthermore, the obtained data set was divided randomly into the training set (63 compounds, 75%) and the test set (21 compounds, 25%) for three random splits. In addition, the published method was used to verify the normality of the activity distribution .

The Monte Carlo optimization method was applied as the main algorithm for the conformation-independent QSAR models development. A hybrid approach was employed for the selection of adequate molecular descriptors as this approach represents the combination of the best features both from the molecular graph and from the descriptors based on the SMILES notation.

Molecular graph-based descriptors were derived from local graph invariants, and the most elementary theoretical graph concepts, such as paths and walks, were used in this research. Literature contains detailed mathematical definitions and applications in QSAR model development (Stoičkov et al. 2018 ). The following local graph invariants were applied for the purpose of defining the optimal descriptors used in conformation-independent QSAR models: Morgan extended connectivity indices of increasing orders (EC0), path numbers with the length of 2 and 3 (p 2 , p 3 ), valence shells in the range of 2 and 3 (s 2 , s 3 ), as well as the number of carbon atom neighbors (NumberOfCarbon), and the non-carbon atom neighbors (NumberofNonCarbon).

Medicinal chemists find descriptors with mechanistic interpretations very appealing, in addition to those that could be correlated with molecular fragments. Unfortunately, most molecular graph-based descriptors do not have this feature, and the SMILES notation was used as a highly convenient alternative to them. One of the key features of the SMILES notation-based descriptors is the possibility of relating them to the appropriate molecular fragment. Data concerning the detailed description of the SMILES notation, as well as its application for defining molecule structure can be found in literature (Toropova et al. 2016; Veselinović et al. 2015) . The descriptors based on the SMILES notation need a numerical value that is to be used in further QSAR model development, as is the case with all molecular descriptors. The numerical value used for this purpose is defined as correlation weight. Its mathematical representation is the sum of all the defined SMILES descriptors, according to Eq. 1.

Z, x, y, t, α, β and γ can have the values of 1 (yes) or 0 (no). These values represent coefficients based on which these descriptors are used in QSAR modeling, which is done by following a simple rule-if the coefficient is 1, the descriptor is used, and if it is 0, then it is discarded. In Eq. 1, the S k symbol is related to the local descriptors which are associated with one SMILES notation symbol (or two that cannot be separated)-SMILES atoms. The new optimal descriptors, represented with the SS k and SSS k symbols, respectively, are defined by the linear combinations of two and three SMILES atoms. With the exception of the local ones used in the QSAR modeling, the global SMILES notation descriptors could be used. When compared to the local ones, the global descriptors differ since they are related to the global features of a studied molecule, not only its structure. The following global SMILES notation-based descriptors were used in the presented research: NOSP, HALO, BOND, and ATOMPAIR. These are defined on the basis of the presence or absence of the following chemical elements: nitrogen, oxygen, sulfur and phosphorus (NOSP); fluorine, chlorine and bromine (HALO); double, triple, or stereochemical bonds (BOND), as well as the presence of seven chemical elements: F, Cl, Br, N, O, S, and P (ATOMPAIR). The presence/absence of eight chemical elements (fluorine, chlorine, bromine, iodine, nitrogen, oxygen, sulfur and phosphorus) and different types of chemical bonds (stereo chemical bond, double bond and triple bond) is defined as with the HARD-index a global SMILES notation-based descriptor and represented as a line of eleven symbols. A combination of both SMILES notation (both local and global) and local graph invariant descriptors is utilized in the presented research for the QSAR model deDCW, and is defined as the sum of all the used optimal descriptors' correlation weights, according to Eq. 2.

(1) DCW(T, Nepoch) = zCW(ATOMPAIR) + xCW(NOSP) + yCW(BOND) + tCW(HALO)

where S k is one SMILES symbol, defined as the SMILES atom, SS k and SSS k are the linear combinations of two or three neighbor SMILES atoms, EC0 k is the Morgan connectivity index of the zero order (a hydrogen suppressed graph was used in this research), PT2 k and PT2 k are paths in the length of 2 and 3, VS2 k and VS3 k are valence shell 2 and 3, and NNC k is Nearest Neighbors (Toropov et al. 2003) .

The calculations of all the above mentioned molecular descriptors and the QSAR model development were conducted with the use of the Monte Carlo optimization method and its algorithm, implemented in CORAL (CORrelation and Logic) (http:// www. insil ico. eu/ coral) software. The Monte Carlo optimization method depends on two essential parameters for the QSAR model development-the threshold (T), used for eliminating rare molecule features, and the number of epochs (N epoch ), which represents the number of iterative processes of the algorithm for the purpose of reaching the top correlation coefficient value. Within the QSAR model developing process, the Monte Carlo optimization method matches molecular descriptors with their numerical values (CW) until CWs are determined. In this manner, the Least Squares method can calculate the DCW (T,N epoch ) for the training and test set compounds, as defined in Eq. 3. A systematic review of the applied method is described in the literature (Toropova et al. 2016; Veselinović et al. 2015) . The search for the best combination of T and N epoch in the presented research was performed within the values of 0-10 for T, and 0-70 for N epoch .

The quality of the developed QSAR model needs to be evaluated by employing different validation methods. The process should provide information on the benefits of the developed model and whether it could be used for future predictions of the studied activities. The methodology from the published paper was used entirely for this purpose (Veselinović et al. 2015; Roy et al. 2008) , which involved the calculation of various statistical parameter values, such as the regular and cross-validated correlation coefficient, standard estimation error, mean absolute error (MAE), the Fischer ratio, root-mean-square error, R m 2 , and MAE-based metrics. What is more, the validation of the developed QSAR models in this research was achieved through data randomization (Y-scrambling test) and with the determination of concordance correlation coefficient (CCC), as well as the novel parameter, dubbed the index of ideality (Stoičkov et al. 2018; Toropov and Toropova 2017; Veselinović et al. 2018 ).

For molecular docking studies were applied the Molegro Virtual Docker (MVD) software. The studies molecules were drawn using the Marvin sketch (Marvin 6.1.0, 2013, ChemAxon), whereas the MMFF94 force field was utilized to gain their optimal 3D geometry. The protein databank provided the structure of the studied enzyme with the PDB id: 6lu7 (the crystal structure of COVID-19 main protease in a complex with an inhibitor N3) ). The purpose of MVD was obtaining an adequate geometrical orientation of the flexible ligand within the active site of the studied enzyme, surrounded by rigid amino acids, as well as identifying the hydrogen bonds and hydrophobic interactions between them. Finally, MVD was also used to calculate the relevant binding energies, also defined as "scoring" functions (Thomsen and Christensen 2006) . The use of these functions is assessing the studied molecules' inhibitory effect, and the following "scoring" functions were calculated in this research-Hbond, NoHbond, VdW, Steric, Pose energy, MolDock, and Rerank Score. The energy from hydrogen and no hydrogen bond interactions were calculated with HBond and NoHbond90, respectively; the energies from the Van der Walls and Steric interactions were calculated with VdW and Steric "scoring" functions; the overall energy of the best-calculated pose was calculated with the Pose energy. MolDock Score and Rerank Score were calculated as the final estimators of ligand and amino acids from the enzyme's active site interaction energies, and the whole molecular docking protocol was validated in accordance with the published methodology (Amin et al. , 2019 Jain et al. 2020; Zivkovic et al. 2020 ). All the crystallized water molecules were removed before the docking. The binding site was computed within spacing in such a way that it was well sampled with the grid resolution of 0.3 Å. The MolDock SE was used as a search algorithm, with the set number of runs being up to 100. The docking procedure parameters were: population size − 50; maximum number of iterations − 1.500; energy threshold − 100.00, and the maximum number of steps − 300. The maximum number of poses to generate was increased to 10 from a default value of 5. Discovery Studio Client v20.1.0.19 was used for showing two-dimensional representations of the interactions between the studied molecules and the amino acids 3CLpro active site.

To consider any QSAR model for prediction purposes, its applicability domain (AD) should be defined prior to use (Gadaleta et al. 2016; Gramatica 2007) . The methodology described in the literature was used to define the AD for all the developed QSAR models in this study (Toropova et al. 2016; Veselinović et al. 2015) . Based on the results obtained, all the molecules fall within the defined AD, and there were no outliers. The numerical values for all the metrics used for the determination of the developed QSAR model's quality are cited in Table 1 . On the basis of the results presented, the Monte Carlo optimization method was able to produce a QSAR model which exhibited good reproducibility, and a high predictability potential. When it comes to the 3CLpro inhibitory activity, the second split yielded the best QSAR model, developed with the T value of 1, whereas the best N epoch value was 20. Figure 1 presents the graphical representations of the best Monte Carlo optimization runs (the highest value for r 2 ) for the developed QSAR models (all the three splits). Equation 4 contains the mathematical representation of the best-developed model. The reproducibility concordance correlation coefficient (CCC) was used for the assessment of the developed QSAR models, and the obtained numerical values for the CCC indicate that the presented QSAR models show high reproducibility values. The metrics results based on the MAE indicate a GOOD model, thus classifying the developed QSAR model as valid. The numerical values of the index of ideality of correlation (IIC) were calculated for the purpose of making a final estimation regarding the quality of the developed QSAR models, and the obtained results exhibit a high predictive potential for the developed QSAR models. The numerical values calculated for Y-randomization, in which the Y values were scrambled in 1000 trials in ten separate runs, are presented in Table 2 , and these were used to evaluate the sturdiness of the developed QSAR models. The presented results indicate that the developed QSAR models were free from correlation by chance because the numerical value for C R 2 p was higher than 0.5.

A QSAR model was developed by Kumar and Roy for the same activity, with the use of 2D molecular descriptors and the multiple linear regression (MLR) technique (Kumar and Roy 2020) . When compared with this research, where the authors used 69 molecules (16 molecules with outlier behavior were omitted from the original database), all the molecules were used in the presented research as they all fall into the defined AD. The following validation parameters (4) pIC 50 = −4.9496(± 0.0040) + 0.0096 × DCW(1, 20) Table 1 The statistical quality of the developed QSAR models for 3CLpro inhibitory activity Table 1 , the QSAR models developed with the application of the Monte Carlo optimization method had Avg r m2 and Δr m2 values of 0.8966-0.6097 and 0.1885-0.0191, respectively. When comparing the two approaches, by applying solely the reported statistical parameters, in terms of predictability, the QSAR model developed with the application of the Monte Carlo optimization method is superior.

The determination of molecular fragments, defined as the SMILES notation optimal descriptors having a positive and negative impact on the studied activity, was one of the chief aims of this research (Halder 2018; Kumar et al. 2019; Manisha et al. 2019; Toropov et al. 2019) . The full list of molecular descriptors, both based on the molecular graph and the SMILES notation, are shown in Table S2 (Supplementary material), while the selected Table 3 . The calculation example is presented in Table S2 , where molecular graph-based descriptors were omitted with the aim of achieving an easier interpretation. The selected molecular fragments, defined as the SMILES notation descriptors which have a positive influence on the studied activity (a lower concentration of the compound is required to achieve IC 50 ) are as follows: "(… (…….", "(…N…(…" and "(…O…(…"-additional branching in the molecule, and additional branching in the molecule realized by adding aliphatic nitrogen and oxygen atom; " + + + + CL-N = = = ", " + + + + N--O = = = " and " + + + + O--B2 = = "-the presence of chlorine and nitrogen atoms in the molecule, the presence of oxygen and nitrogen atoms in the molecule, and the presence of oxygen and a double bond in the molecule; "BOND10000000", "C…/…….", and "C…/…C…"-double bond between two carbon atoms; "C…C…C…"-propyl group; "N…C…C…"-ethyl amine group; "O…1……." oxygen atom inside the ring; "O…C…….", "O…C…C…" and "C…O…C…"-methoxy, ethoxy and dimethoxy groups, respectively; "s…c……." and "s……….."-aromatic sulfur atom; "c…c……." and "c…c…c…"-aromatic carbon atoms. The same analysis may be performed for molecular fragments, defined as the SMILES notation descriptors with a negative impact on the studied activity. The selected molecular fragments are: "C…(…….", "N…C…(…" and "N…(……."-simple molecular branching on the carbon and nitrogen atom; "(…Br..(…" and "(…Cl..(…"-additional molecular branching with bromine and chlorine atoms; "C…2…(…" and "N…2……."-second ring inside the molecule, and second ring inside the molecule with a nitrogen atom; " = …1……."-double bond within the ring; "N……….." and "-aliphatic nitrogen; "N… = …….", "N… = …C…", "O… = ……." and "O… = …C…"-the double bond between the nitrogen or oxygen atom and the carbon atom; "n…n…….", "s…n……." and "s…n…n…"aromatic nitrogen atoms and the presence of both aromatic nitrogen and sulfur atoms. The graphical representation of the impact of molecular fragments on the selected molecule activity is displayed in Fig. 2 . The green color in Fig. 2 indicates a positive influence, while red indicates a negative impact on the studied activity.

Identified molecular fragments were used for the computer-aided design of novel 3CLpro inhibitors. The molecule with the highest inhibitory effect (compound 1 from the dataset) was used as a starting/template molecule for the purpose of producing molecules with a higher inhibition potential. Table 4 indicates the SMILES notation of all the designed molecules, along with their calculated values for the -pIC50, utilizing the best-developed QSAR model, while their chemical structures are presented in Fig. 3 .

Added molecular fragments were identified as structural alerts that have a positive effect on the studied activity in all the designed molecules (A1-A5). Molecules A1-A3 had the following molecular fragments added-the isopropyl and cyclopropyl group on the different parts of template molecule A. The groups are defined as the following SMILES notation descriptors, all promoters of the studied activity increase-"(…(……."; "1……….."; "C……….."; "C…C……."; "C…C…C…". Molecules A4 and A5 have added methyl ether group. They all have a positive impact on the studied activity and are defined with the following SMILES notation descriptors -"C…O…C…"; "O…C……."; "O…C…C…"; "O………..". The approach suggested in the literature was taken to assess the predictability of the developed QSAR model, and to further estimate inhibitory potential of designed molecules (Zivkovic et al. 2020) . Within this approach all designed molecules and template molecule A were subjected to molecular docking studies with SARS-CoV-3CLpro enzyme to assess their binding potential and numerical values for all calculated "scoring" functions are presented in Table 5 . One of assumption is that the more molecule is bounded to enzyme, indicated with higher binding energies between molecule and amino acids form active site and calculated in regards to appropriate "scoring" functions, the higher inhibition would be, resulting to higher values of pIC 50 , smaller concentrations are needed to achieve inhibition. Also, different ligand-amino acids interactions could be associated with various scoring functions, so everything should be taken into consideration when performing the assessment of the inhibitory potency. Literature contains the detailed definitions of other "scoring" functions (Thomsen and Christensen 2006, Zivkovic et al. 2020) . Calculated values for energies with the application of defined "scoring functions" are much higher in comparison to real binding energies and most likely they are unrealistic. However, since same approach was used for all molecules, obtained values could be compared to assess the binding preferences, meaning the higher the binding energy, the higher binding preferences, leading to higher inhibitory potential. According to the obtained results for MolDock and ReRank, molecule A5 has the potentially highest inhibitory activity and molecule A the lowest. This result is in correlation with the results from the QSAR modeling since the calculated values for -pIC 50 using the best QSAR model obtained with the Monte Carlo optimization method are the most preferable for molecule A5 and the least preferable for molecule A. All the interactions between the selected molecules and amino acids from the 3CLpro enzyme's active site have been identified, and Figures S1-S6 in the Supplementary Information section present the 2D representation of hydrogen bonds, in addition to hydrophobic, and hydrophilic interactions inside the 3CLpro binding pocket. The best-calculated poses for all studied molecules inside the active site of 3CLpro are cited in Fig. 4 .

It is of extreme importance to determine the studied compounds' drug-likeness when the aim is for them to become therapeutics. Compounds should possess physical and chemical parameters enabling them to be considered as potential medication. The application of Lipinski's "Rule of Five" represents one of the extensively used approaches. These rules determine the absorption/permeation of the compound, and if the compound possesses a molecular weight higher than 500, logP over 5, and the number of hydrogen bond donors and acceptors is higher than 5 and 10, respectively, it will most likely have poor absorption/permeation. In addition, the number of rotatable bonds, associated with conformational molecular flexibility, is used to assess the efficacy of binding to receptors/channels, as well as their bioavailability. The molecule should have a satisfactory number of rotatable bonds, i.e. 10 or fewer (Lipinski et al. 2001 ). In addition to Lipinski's "Rule of Five", Veber proposed having an additional set of rules for a molecule's drug-likeness assessment (Veber et al. 2002) , and they can be summarized as acceptable when the oral bioavailability is more likely to occur when the molecule has 10 or fewer rotatable bonds, a polar surface area equal to or less than 140 Å 2 , and 12 or fewer hydrogen donors and acceptors. For template molecule A and all the designed molecules, the physical and chemical parameters defined above are calculated using the Molinspiration Cheminformatics software (www. molin spira tion. com), and presented in Table 6 . In accordance with the obtained results, all the molecules were in accordance with Lipinski's "Rule of Five", and the Veber rules, which may be associated with their potentially good solubility and permeability through biological membranes, thus leading to satisfactory bioavailability, which is why they may be considered as potential therapeutics.

The present research demonstrates the development of new robust and reliable QSAR models for 3CLpro inhibition based on the Monte Carlo optimization method as the main model, and developed on the basis of optimal descriptors derived both from a local graph and the SMILES notation invariants. The assessment of the developed QSAR model's robustness and its predictive potential was achieved by applying various statistical techniques and obtained numerical values, used to validate the developed QSAR model, indicating that it has high applicability. The applied methodology was able to determine the molecular fragments, used as the SMILES notation fragments in QSAR modeling, which have a positive and negative impact on 3CLpro inhibition. They were used for the computer-aided design of novel 3CLpro inhibitors with a potentially higher activity when compared to template molecule A. To assess the developed QSAR models' predictability all designed molecules and template molecule A were subjected to molecular docking studies with 3CLpro. Their binding potential was assessed according to the comparison of numerical values for all calculated "scoring" functions, related to energies from the interactions between the molecules and the amino Fig. 4 The best-calculated poses for designed molecules (A1-A5) and template molecule A inside the active site of 3CLpro 50 show high inter-correlation. Also, all designed molecules have potentially good pharmacokinetic properties and possess drug-likeness. The methodology presented in this research paper can be applied in the search for novel therapeutics for the treatment of COVID-19 based on 3CLpro inhibition.

The online version contains supplementary material available at https:// doi. org/ 10. 1007/ s11696-022-02170-8.

Design and evaluation of anti-SARS-coronavirus agents based on molecular interactions with the viral protease

Fight against novel coronavirus: a perspective of medicinal chemists

Exploring pyrazolo[3,4-d]pyrimidine phosphodiesterase 1 (PDE1) inhibitors: a predictive approach combining comparative validated multiple molecular modelling techniques

Reliable structural information for rational design of benzoxazole type potential cholesteryl ester transfer protein (CETP) inhibitors through multiple validated modelling techniques

Protease targeted COVID-19 drug discovery: what we have learned from the past SARS-CoV inhibitors?

Enzymatic activity characterization of SARS coronavirus 3C-like protease by fluorescence resonance energy transfer technique 1

Emerging coronaviruses: genome structure, replication, and pathogenesis

QSAR modeling: Where have you been? Where are you going to

Coronaviruses resistant to a 3C-like protease inhibitor are attenuated for replication and pathogenesis, revealing a low genetic barrier but high fitness cost of resistance

In silico pharmacology for drug discovery: applications to targets and beyond

Applicability domain for QSAR models: where theory meets reality

BindingDB: a public database for medicinal chemistry, computational chemistry and systems pharmacology

Principles of QSAR models validation: internal and external

The origin, transmission and clinical therapies on coronavirus disease 2019 (COVID-19) outbreak-An update on the status

Finding the structural requirements of diverse HIV-1 protease inhibitors using multiple QSAR modelling for lead identification

SARS-CoV-2 cell entry depends on ACE2 and TMPRSS2 and is blocked by a clinically proven protease inhibitor

Discovery of non-covalent inhibitors of the SARS main proteinase 3CLpro

Chemical-informatics approach to COVID-19 drug discovery: Monte Carlo based QSAR, virtual screening and molecular docking study of some in-house molecules as papain-like protease (PLpro) inhibitors

Structure of Mpro from SARS-CoV-2 and discovery of its inhibitors

Persistence of coronaviruses on inanimate surfaces and its inactivation with biocidal agents

Development of a simple, interpretable and easily transferable QSAR model for quick screening antiviral databases in search of novel 3Clike protease (3CLpro) enzyme inhibitors against SARS-CoV diseases

In silico design of diacylglycerol acyltransferase-1 (DGAT1) inhibitors based on SMILES descriptors using Monte-Carlo method

Corona virus: A review of COVID-19

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) and coronavirus disease-2019 (COVID-19): the epidemic and the challenges

Anti-SARS coronavirus 3C-like protease effects of Isatis indigotica root and plant-derived phenolic compounds

Experimental and computational approaches to estimate solubility and permeability in drug discovery and development settings

Current mathematical methods used in QSAR/ QSPR studies

Research and development on therapeutic agents and vaccines for COVID-19 and related human coronavirus diseases

Development of prediction model for fructose-1,6-bisphosphatase inhibitors using the Monte Carlo method

Comparative QSARs for antimalarial endochins: importance of descriptor-thinning and noise reduction prior to feature selection

Further exploring rm2 metrics for validation of QSPR models

SARS-CoV-2: a storm is raging

Transmission routes of 2019-nCoV and controls in dental practice

Variables selection methods in QSAR: an overview

Synthesis and evaluation of pyrazolone compounds as SARS-coronavirus 3C-like protease inhibitors

The epidemiology and pathogenesis of coronavirus disease (COVID-19) outbreak

Exploring the impact of size of training sets for the development of predictive QSAR models

Pharmacologic treatments for coronavirus disease 2019 (COVID-19): a review

Structural basis of receptor recognition by SARS-CoV-2

A review of coronavirus disease-2019

World Health Organization declares global emergency: a review of the

Development of non-peptide ACE inhibitors as novel and potent cardiovascular therapeutics: an in silico modelling approach

Computer-aided drug design and drug pharmacokinetic prediction: a mini-review

MolDock: a new technique for high-accuracy molecular docking

The index of ideality of correlation: a criterion of predictive potential of QSPR/QSAR models?

Structure-toxicity relationships for aliphatic compounds based on correlation weighting of local graph invariants

Idealization of correlations between optimal simplified molecular inputline entry system-based descriptors and skin sensitization

The utilization of the Monte Carlo technique for rational drug discovery

Molecular properties that influence the oral bioavailability of drug candidates

COVID-19: a promising cure for the global panic

Application of smiles notation based optimal descriptors in drug discovery and design

Design and development of novel antibiotics based on FtsZ inhibition-in silico studies

The SARS-CoV-2 outbreak: what we know

Clinical findings in a group of patients infected with the 2019 novel coronavirus (SARS-Cov-2) outside of Wuhan, China: retrospective case series

Crystal structure of SARS-CoV-2 main protease provides a basis for design of improved α-ketoamide inhibitors

Potential COVID-2019 3C-like protease inhibitors designed using generative deep learning approaches

The application of the combination of Monte Carlo optimization method based QSAR modeling and molecular docking in drug design and development

Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations

We have no conflicts of interest to disclose.