Structure features of peptide-type SARS-CoV main protease inhibitors: Quantitative structure activity relationship study.

In the present work, an extensive QSAR (Quantitative Structure Activity Relationships) analysis of a series of peptide-type SARS-CoV main protease (MPro) inhibitors following the OECD guidelines has been accomplished. The analysis was aimed to identify salient and concealed structural features that govern the MPro inhibitory activity of peptide-type compounds. The QSAR analysis is based on a dataset of sixty-two peptide-type compounds which resulted in the generation of statistically robust and highly predictive multiple models. All the developed models were validated extensively and satisfy the threshold values for many statistical parameters (for e.g. R2 ​= ​0.80-0.82, Q2 loo ​= ​0.74-0.77, Q 2 LMO  ​= ​0.66-0.67). The developed QSAR models identified number of sp2 hybridized Oxygen atoms within seven bonds from aromatic Carbon atoms, the presence of Carbon and Nitrogen atoms at a topological distance of 3 and other interrelations of atom pairs as important pharmacophoric features. Hence, the present QSAR models have a good balance of Qualitative (Descriptive QSARs) and Quantitative (Predictive QSARs) approaches, therefore useful for future modifications of peptide-type compounds for anti- SARS-CoV activity.

Introduction

Coronaviruses (subfamily Coronavirinae, family Coronaviridae, order Nidovirales) have been classified into four genera: Alphacoronavirus, Betacoronavirus, Gammacoronavirus and Deltacoronavirus [1]. Of these, especially Betacoronavirus, have been found to cause respiratory, enteric, hepatic, and neurological diseases in many animals, and also in humans [2]. The two Betacoronaviruses (βCoVs) viz. severe acute respiratory syndrome CoV (SARS-CoV) and Middle East respiratory syndrome CoV (MERS-CoV) spread in 2003 and 2012, respectively [2,3]. SARS-CoV and MERS-CoV have fatality rate of 10% and 35%, respectively [2,3]. Unfortunately, till this date, there is no appropriate treatment for SARS-CoV and MERS-CoV [2,3]. The situation has worsened with the recent outbreak of more contagious novel coronavirus SARS-CoV-2, which is the causative agent for COVID-19 pandemic. This disease has a long-term socio-economic impact on many countries due to high infection and mortality rate. At present, this disease is responsible for more than 732,499 deaths and 20 million confirmed infected cases [4]. Therefore, there is an urgent need to curb this deadly disease.

The recent outbreak of COVID-19 pandemic, with no approved treatment for infections, has spread quickly in many countries. Even though, COVID-19 is a new disease caused by the novel coronavirus SARS-CoV-2, but a good number of studies suggest that it has significant similarity with SARS-CoV [3,[5], [6], [7]]. This similarity is reflected from the facts that SARS-CoV and SAR-CoV-2 have a 79% similarity at the genome level, in addition, a lot of proteins like glycosylated spike (S) protein (76% of sequence similarity and a highly conserved receptor-binding domain), papain-like protease (PLpro) (83% similarity with similar active sites), RNA-dependent RNA polymerase (RdRp), and coronavirus main protease (3CLpro) are essential for both. To add further, RdRp and 3CLpro protease of SARS-CoV-2 share over 95% of sequence similarity with those of SARS-CoV [3,[5], [6], [7]]. In short, viral proteins essential for SARS-CoV-2 entry into host cells and subsequent replication are highly similar to those associated with SARS-CoV. Consequently, research and development on SARS-CoV could be useful for the development of a therapeutic or preventive agent for COVID-19. To add further, designing and synthesis of a new drug and its congeners followed by their bio-screening need a lot of time. Therefore, high similarity of SARS-CoV-2 with SARS-CoV and considering the potential threat of this pandemic, it is reasonable that a therapeutic drug, which has been previously tested against SARS-CoV could be easily optimized to be effective against SARS-CoV-2.

To optimize a compound to become a lead or a drug, a researcher needs to follow an easy, efficient, economical, and eco-friendly approach (e-Chemistry approach). A feasible solution to achieve these goals is to use Computer-Aided Drug Designing (CADD). CADD is a contemporary approach with a good number of benefits like cheaper, result oriented, minimizes animal testing as well as trials and errors, less time-consuming, a few advantages to mention. In recent time, the thriving branches of CADD such as QSAR (Quantitative Structure Activity Relationships), Molecular docking, Pharmacophore modelling, etc. have contributed significantly in optimization of lead and drug candidates [8,9].

QSAR analysis involves finding mathematical correlation between the structural features of congeneric molecules with bio-activity. The main steps followed during a typical QSAR analysis are (1) Selection of an appropriate dataset (2) 3D-structure generation and their optimization using suitable technique (3) Molecular descriptor calculation followed by their pruning (4) Generation of QSAR model using a proper algorithm for feature (that is, molecular descriptor) selection and (5) Adequate validation of developed QSAR model [10].

QSAR is an established CADD branch, which has been used successfully to identify noticeable and hidden structural patterns/features having correlation with a desired activity/property (Qualitative QSAR). In addition, it is helpful to predict the activity/property before the actual synthesis and testing of a molecule (Quantitative QSAR). Therefore, QSAR analysis is performed routinely by the researchers to get important qualitative and quantitative idea about the congeneric molecules for optimization. An adequately validated and statistically robust QSAR analysis offers in-depth knowledge about the structural patterns that have good correlation with the desired activity/toxicity/property of a drug candidate and enhanced intuition for the mechanism of drug action [[11], [12], [13], [14], [15], [16], [17], [18]].

A good number of researchers have synthesized peptide-type compounds and tested against SAR-CoV [[19], [20], [21], [22]]. Recently, Zhang et al. [23] synthesized and tested α-Ketoamides (peptide-type compounds) as broad-spectrum inhibitors of coronavirus (SARS-CoV). Despite these efforts, optimization of peptide-type compounds is in pipeline to get a drug candidate. Henceforth, in the present work, we have developed Qualitative cum Quantitative SAR models for a series of sixty-two peptide-type compounds for anti- SARS-CoV activity. The results could be beneficial for future optimizations of peptide-type compounds with better activity profile.

Experimental methodology

Selection of dataset: The dataset selected for the present work comprises sixty-two peptide-type compounds having moderate to high activity against SARS-CoV [19,21,22]. The selected peptide-type compounds possess Ki ​= ​3 to 56,000 ​nM. The reported activity Ki values were converted to pKi (pKi ​= ​–logKi) before actual QSAR analysis. For the sake of ease and understanding the chemical space covered by the present dataset, most and least active three molecules have been depicted in Fig. 1 . The SMILES notation for all selected molecules along with their reported activity values Ki, and pKi are present in Table S1 in supplementary material.

Fig. 1

Variations in activity and chemical structure in the present dataset.

Structure optimization and molecular descriptor calculation

The structures were drawn using ChemSketch 12 Freeware (www.acdlabs.com) followed by their conversion to 3D-structures using OpenBabel 2.4. Then, the force field MMFF94 was used to optimize the structures. The 3D-structures were optimized using default settings available in TINKER. The molecules in each set were then aligned using Open3DAlign. In the next step, the optimized and aligned molecules were then used for calculation of molecular descriptors using PyDescriptor [24] and PaDEL [25].

Molecular descriptor pruning

PyDescriptor and PaDEL provided more than 30,000 molecular descriptors for each molecule in all sets. Therefore, molecular descriptor pruning was essential to remove redundant molecular descriptors. For this, molecular descriptors with high co-linearity (|R|>0.90) and nearly constant (>95%) were removed to avoid the inclusion of multi-collinear and spurious variables in GA-MLR (Genetic Algorithm–Multi-linear Regression) model, using objective feature selection in QSARINS ver. 2.2.2 [[26], [27], [28], [29]]. The resulting reduced molecular descriptor pool comprised of only 668 molecular descriptors only but large enough to cover 1D-to 3D-descriptor space.

QSAR model building and their validation:

The reduced molecular descriptor set comprises zero-, one-, two- and three-dimensional descriptors, charge descriptors and molecular properties, thereby covering broad descriptor space. Subjective feature selection (SFS) was executed to build the statistically acceptable GA-MLR based QSAR models using QSARINS ver. 2.2.2. The developed models were subjected to thorough statistical validation (internal and external validation) according to OECD principles; models with high internal and external predictive ability have been reported.

The general procedure for building the QSAR models is as follows:

The set was bifurcated randomly, using random splitting option in QSARINS, into a training and a prediction set of 52 (i.e. 80 % training set) and 12 (i.e. 20 % prediction set), respectively. Then, the training set was used for model development, and the prediction set for the external validation, that is, to judge the predictive ability on new chemicals.

QSARINS was used to build GA-MLR based QSAR models using default settings. The selected fitness function to maximize in GA was Q2, this ensured the double cross-validation as well [26,27,30]. During model development it was observed that there was growth in the value of Q2 up to six variables, but then, it had visible and significant reduction (see Fig. 2 ). Therefore, molecular descriptor selection was restricted to a set of six descriptors to avoid over-fitting and develop easy and informative QSAR models.

Fig. 2

Correlation of statistical parameters with number of molecular descriptors.

The values for molecular descriptors, which are present in QSAR models, are available in the supplementary information for each molecule.

One of the OECD guideline suggests to thoroughly validating a QSAR model, therefore all the models were subjected to internal and external validation, Y-scrambling along with model applicability domain (AD) analysis using QSARINS. The statistical quality and strength of a GA-MLR based QSAR model was judged on the basis of: (a) internal validation based on leave-one-out (LOO) and leave-many-out (LMO) procedure (i.e. cross-validation (CV)); (b) using External validation; (c) Y-randomization (or Y-scrambling) and (d) fulfilling of respective threshold value for the statistical parameters [31,32]: R 2 tr ​≥ ​0.6, Q 2 loo ​≥ ​0.5, Q 2 LMO ​≥ ​0.6, R 2 ​> ​Q 2, R 2 ex ​≥ ​0.6, RMSE  ​< ​RMSE , ΔK ​≥ ​0.05, CCC ​≥ ​0.80, Q 2-F n ​≥ ​0.60, r 2 m ​≥ ​0.6, (1-r 2/r o 2) ​< ​0.1, 0.9 ​≤ ​k ​≤ ​1.1 or (1-r 2/r’o 2) ​< ​0.1, 0.9 ​≤ ​k’ ​≤ ​1.1,| r o 2− r’o 2| ​< ​0.3 with RMSE and MAE close to zero. A QSAR model that did not satisfy above mentioned criteria was consequently excluded.

Thus, the complete procedure involving molecular descriptor calculation and their pruning, followed by subjective feature selection along with model building and validation was performed on all the four sets. It was observed that the MMFF94 optimized set resulted in development of statistically better QSAR models, which have been reported in the present work.

Result and discussions

Though, the dataset used in the present study contains only sixty-two compounds but the presence of positional isomers, heterocyclic rings, etc. significantly augment the chemical space covered by the peptide-type compounds. In our previous work on QSAR analysis related to small and moderate sized datasets [29], we have demonstrated that a QSAR model built using undivided whole dataset provides advantages like identification of maximum useful information, capturing of maximum relevant molecular descriptors, and benchmark for comparison and assessment of QSAR models constructed using divided datasets. In addition, this approach also helps to capture unrevealed structural features, which govern the bio-activity profile of congeneric molecules. Therefore, in the present work, models have been derived using divided and undivided datasets.

The derived QSAR models are as follow:

Model-1.1 (Undivided Set model):

pKi ​= ​1.781 (±1.357) ​+ ​1.073 (±0.317) ∗ sp2O_aroC_7B ​+ ​0.377 (±0.111) ∗ APC2D3_C_N −1.264 (±0.647) ∗ APC2D9_N_N ​+ ​−2.538 (±0.751) ∗ KRFPC3478 ​+ ​−0.743 (±0.233) ∗ fringNsp3C8B −0.61 (±0.246) ∗ APC2D6_C_S.

Model-1.2 (Divided Set model: Training: 80%, Prediction: 20%):

pKi ​= ​1.7 (±1.658) ​+ ​1.182 (±0.353) ∗ sp2O_aroC_7B ​+ ​0.363 (±0.115) ∗ APC2D3_C_N −1.351 (±0.642) ∗ APC2D9_N_N −2.422 (±0.731) ∗ KRFPC3478–0.772 (±0.269) ∗ fringNsp3C8B −0.169 (±0.085) ∗ ringC_sp3S_9B

The statistical parameters for developed models 1.1 and 1.2 have been presented in Table 2 . The symbols have their usual meaning, which are available in the supplementary material also.

Table 2

Statistical parameters for developed QSAR models 1.1 and 1.2

Statistical Parameter	Model-1.1	Model-1.2
Fitting
R2tr	0.801	0.824
R2adj.	0.78	0.8
R2tr- R2adj.	0.022	0.025
LOF	0.326	0.324
Kxx	0.241	0.271
ΔK	0.06	0.05
RMSEtr	0.461	0.433
MAEtr	0.375	0.328
RSStr	13.15	9.366
CCCtr	0.89	0.904
s	0.489	0.467
F	36.953	33.606
Internal validation
R2cv(Q2loo)	0.741	0.769
R2-R2cv	0.06	0.055
RMSEcv	0.526	0.496
MAEcv	0.427	0.378
PRESScv	17.14	12.312
CCCcv	0.857	0.874
Q2LMO	0.673	0.665
R2Yscr	0.097	0.124
Q2Yscr	−0.182	−0.22
External validation
RMSEex	–	0.527
MAEex	–	0.461
PRESSext	–	3.332
R2ex	–	0.758
Q2-F1	–	0.741
Q2-F2	–	0.74
Q2-F3	–	0.739
CCCex	–	0.841
R2-ExPy	0.743	0.758
R’o2	0.684	0.623
k’	0.996	0.975
1-(R2/R’o2)	0.08	0.178
r’2m	0.562	0.48
Ro2	0.741	0.753
k	0.997	1.019
1-(R2-ExPy/Ro2)	0.002	0.006
r2m	0.713	0.709

A good number of statistical parameters for model 1.1 and 1.2, which are related to fitting, internal and external validation and Y-scrambling, have been tabulated in Table 2. From Table 2, it is clear that R2 tr, CCCtr, CCCcv, R2 adj. and F satisfy the recommended threshold value, which shows that the QSAR models are statistically robust with adequate number of molecular descriptors in the models. The values for different cross-validation parameters such as R2 cv, RMSEcv, MAEcv, CCCcv, and Q2 LMO support the statistical robustness of the QSAR models. The external predictive ability of the models is established by the high values of R2 ex, Q2F1, Q2F2, Q2F3, and CCCex [12,14,15,[33], [34], [35]].

In short, the developed QSAR models fulfill the recommended threshold values for many internal and external validation parameters. In addition, for a better validation of derived models, the model applicability domain (AD) was assured by plotting Williams plots for models 1.1 and 1.2 (see Fig. 3 ). Therefore, these models are statistically robust and possess good external predictive ability. To add further, satisfaction of recommended threshold values for these parameters along with low correlation among the molecular descriptors point out that these models are not developed by chance (see supplementary information).

Fig. 3

(a) Graph of experimental vs Predicted pKi values for model 1.1 (b) Williams plot for model 1.1 (c) Graph of experimental vs Predicted pKi values for model 1.2 (d) Williams plot for model 1.2.

Interpretation of QSAR models

The models 1.1 and 1.2 have been built using the undivided and divided dataset, respectively. They comprise four common molecular descriptors. Therefore, the approach to develop QSAR models using undivided and divided dataset has been successful in identification of greater number of important molecular descriptors, which is useful to capture maximal information. Though we have compared the activities of the molecules of the dataset in terms of a single descriptor, we make it clear that the combined or converse effect of confounding factors/descriptors do have additional influence on the activity profile of the compounds.

A molecular descriptor with a positive coefficient in both the model is sp2O_aroC_7 ​B (number of sp2 hybridized Oxygen atoms within seven bonds from aromatic Carbon atoms). Since, in the present series of compounds, sp2 hybridized Oxygen atoms are always present as a part of carbonyl group (>CO), therefore it is rational to consider that this molecular descriptor also points out toward the presence of number of carbonyl groups in conjugation with aromatic Carbon atoms, that is, aromatic rings. Therefore, increasing the number of carbonyl groups within seven bonds from aromatic Carbon atoms could increase the anti-SARS activity of peptide-type of compounds. This observation is supported by the fact that molecule number 2 (pKi ​= ​5.658 ​M) possesses only two such carbonyl groups, whereas the molecule number 47 (pKi ​= ​8.523 ​M) and 48 (pKi ​= ​8.387 ​M) have three such carbonyl groups. Another such comparison is possible between molecule number 30 (pKi ​= ​5.77 ​M) and 31 (pKi ​= ​7.174 ​M). This descriptor has been depicted in Fig. 4 . The sp2-hybridized Oxygen and aromatic Carbon atoms have been shown using blue colour and the seven bonds are red coloured.

Fig. 4

Representation of sp2O_aroC_7 ​B using molecule number 47 and 2 as representatives only.

The atom-pair molecular descriptor APC2D3_C_N stands for the presence of Carbon and Nitrogen atoms at a topological distance of 3. This molecular descriptor has a positive coefficient in both the models. Therefore, higher value of this descriptor could lead to better activity profile for a molecule. A comparison of 31 (pKi ​= ​7.187 ​M) with 45 (pKi ​= ​4.854 ​M) as well as among 42 (pKi ​= ​7.658 ​M), 43 (pKi ​= ​6.097 ​M) and 45 (pKi ​= ​4.854 ​M) also support this observation. These molecules along with their APC2D3_C_N values have been depicted in Fig. 5 .

Fig. 5

Exemplification of molecular descriptor APC2D3_C_N.

An atom pair molecular descriptor with a negative coefficient in both the models is APC2D9_N_N, which stands for the number of Nitrogen atoms separated from each other by a topological distance of nine. It appears that increasing its value could cause diminish the anti-SARS activity of peptide-type compounds. A comparison of molecule number 51 (pKi ​= ​7.658 ​M) with 54 (pKi ​= ​6.658 ​M) and 55 (pKi ​= ​6.658 ​M) supports this observation (see Fig. 6 ). Therefore, such a combination of Nitrogen atoms should be avoided for better activity.

Fig. 6

Pictorial representation of molecular descriptor APC2D9_N_N.

KRFPC3478 is a finger print molecular descriptor, which represents the presence of Carbon atom at the position number 3 of an Indole ring. This descriptor with a negative coefficient in both models has negative contribution towards the anti-SARS activity of peptide-type compounds. An analysis of molecule 31–42 indicates that the presence of Indole moiety does not have negative contribution each time. Hence, it is rational to consider that the presence of Carbon atom at the position number 3 of an Indole ring has negative contribution. Therefore, it must be avoided to have better activity. The molecule 40 (pKi ​= ​5.174 ​M) and 41 (pKi ​= ​5.125 ​M) have relatively lower activity than other analogues bearing Indole ring (see Fig. 7 ).

Fig. 7

Representation of molecular descriptor KRFPC3478.

Another molecular descriptor with a negative coefficient in model 1.1 and 1.2 and hence with a negative correlation with activity is fringNsp3C8B. This descriptor represents frequency of occurrence of sp3-hybridized Carbon atoms exactly at eight bonds from the ring Nitrogen atoms. As the number of such Nitrogen atoms increases, the activity decreases. This observation is supported on comparing molecule 31 (pKi ​= ​7.187 ​M) and 46 (pKi ​= ​6.167 ​M). These two molecules are positional isomers of each other, but they have good difference in their activity. This could be attributed to higher frequency of occurrence of sp3-hybridized Carbon atoms exactly at eight bonds from the ring Nitrogen atoms in case of molecule 46 than the molecule 31. This descriptor has been shown in Fig. 8 using molecule number 31 and 46 as representatives. The sp3-hybridized Carbon atoms and ring nitrogen atoms have been highlighted using blue colour, while the eight bonds by red colour.

Fig. 8

Depiction of molecular descriptor fringNsp3C8B using 31 and 46 as representatives.

Though, the molecular descriptors APC2D6_C_S (presence of Carbon and Sulfur at a topological distance of six) and ringC_sp3S_9B (number of ring Carbon atoms within nine bonds from sp3-hybridized Sulfur atoms) convey interrelation of Carbon and Sulfur and their subsequent effect on activity, but both molecular descriptors provide different level and type of information in varying details. Both molecular descriptors have a negative coefficient. Therefore, increasing their value could reduce the anti-SARS activity profile. A comparison of molecule number 42 (pKi ​= ​7.658 ​M) with 43 (pKi ​= ​6.097 ​M) vindicates this observation. The molecule 42 has a lower number of ringC_sp3S_9B and a higher value of APC2D3_C_N than 43.

Conclusions

In conclusion, statistically robust QSAR models with good external predictive ability have been developed, which have successfully highlighted a good number of molecular features. The developed models 1.1 to 2.2 satisfy the threshold values for many statistical parameters that are necessary to establish the quality and usefulness of a QSAR model. Thus, the developed QSAR models have a good balance of Quantitative and Qualitative aspects. Therefore, the developed models could be useful for future optimization of the activity profile of the molecules used in the present dataset.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.