COSMO-RS-Based Descriptors for the Machine Learning-Enabled Screening of Nucleotide Analogue Drugs against SARS-CoV-2.

Chemical similarity-based approaches employed to repurpose or develop new treatments for emerging diseases, such as COVID-19, correlates molecular structure-based descriptors of drugs with those of a physiological counterpart or clinical phenotype. We propose novel descriptors based on a COSMO-RS (short for conductor-like screening model for real solvents) σ-profiles for enhanced drug screening enabled by machine learning (ML). The descriptors' performance is hereby illustrated for nucleotide analogue drugs that inhibit the ribonucleic acid-dependent ribonucleic acid polymerase, key to viral transcription and genome replication. The COSMO-RS-based descriptors account for both chemical reactivity and structure, and are more effective for ML-based screening than fingerprints based on molecular structure and simple physical/chemical properties. The descriptors are evaluated using principal component analysis, an unsupervised ML technique. Our results correlate with the active monophosphate forms of the leading drug remdesivir and the prospective drug EIDD-2801 with nucleotides, followed by other promising drugs, and are superior to those from molecular structure-based descriptors and molecular docking. The COSMO-RS-based descriptors could help accelerate drug discovery for the treatment of emerging diseases.

On March 11, 2020, the World Health Organization declared the disease COVID-19, caused by the new virus severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), as a pandemic.[1] SARS-CoV-2 is a positive, nonsegmented, single-stranded ribonucleic acid (RNA) β-coronavirus that belongs to the large family of human coronaviruses, along with SARS-CoV-1 and Middle East respiratory syndrome coronavirus (MERS-CoV).[2,3] To date, a vaccine and effective therapies for the prevention or treatment of COVID-19 are not available.[4] The research community has vigorously launched or planned numerous COVID-19 research projects and clinical trials.[5] Considering the limited knowledge of the SARS-CoV-2 infection process, researchers are currently developing treatments based on the knowledge of other SARS viruses.[6] The COVID-19 drug discovery effort is to a large extent focused on the inhibition of the angiotensin converting enzyme II (ACE2), assisting in viral entry, and RNA-dependent RNA polymerase (RdRp), catalyzing viral transcription.[7−9] Currently, the leading anti-SARS-CoV-2 drug remdesivir and highly promising EIDD-2801 with active forms designed to inhibit the RdRp-catalyzed RNA replication by mimicking nucleotides are under intensive clinical studies.[10,11]

Drug development faces multifold challenges, such as extended time to bring new drugs to the market, high attrition rates, and changing regulatory requirements.[12,13] The repurposing of market-available or investigational drugs for the treatment of new or rare diseases is increasingly becoming attractive, because it reduces the risk, cost, and timeline for drug development.[13,14] Computational data-driven approaches for drug screening are particularly useful for the identification of drug candidates. The initial identification step is typically followed by a mechanistic assessment in preclinical models and ultimately by clinical trials.[13] Among the computational approaches, signature matching and molecular docking are based on chemical structure analysis, while genome association, pathway mapping, and retrospective clinical analysis involve genome data and electronic health records screening.[15] Chemical and structural similarity-based methods rely on the hypothesis that structural similarity implies correlation at the therapeutic effect or drug-target level.[16] Molecular docking and chemical similarity-based strategies have been employed to study the pharmacological activity and screen antiviral drugs against key enzymes of SARS-CoV-2.[17−19] However, in some cases, computational screening methods are not in full agreement with in vitro, preclinical, and clinical studies,[20] and even the clinical studies of remdesivir have different conclusions.[10,21] Due to the global urgency and large number of possible drugs, more accurate and enhanced computational methods are needed.

One of the potentially promising ways to improve predictability is through the use of artificial intelligence (AI) and machine learning (ML). Increasingly popular in drug development, ML techniques allow the automation of multidimensional data analysis,[22] and have shown great promise for new drug development and repositioning.[12,23,24] Based on a training set, ML approaches enable the screening of efficient drug candidates without having a deep knowledge of the underlying processes. The key to the efficiency of ML methods is the use of representative descriptors or molecular fingerprints to translate the molecular system to a machine-understandable language.

A number of molecular descriptors have been developed and implemented into the leading software packages, such as BIOVIA Discovery Studio 2020 and Schrödinger Suite, which provide good results for drug screening.[25,26] The majority of these descriptors are based on the molecular structure, such as the number of rings, hydrogen bond donors and acceptors, induced charges, etc. It is worth noting that such descriptors are important but not sufficient for the chemical similarity analysis of ligands and its extension to ML, as these yield mixed results regarding the ability to enhance the accuracy of ensemble-averaged quantity prediction.[27] Typically, small drug molecules interact with a protein active site localized on a small surface area in the presence of a complex environment. To capture the full picture of interactions, the ideal drug screening descriptors should combine both structural and chemical thermodynamics-based features. Methods based on ML have successfully been employed to replace the computationally expensive docking and binding energy calculations.[28,29]

In order to better represent the ligand and its interaction with the active site of a protein, the list of descriptors should be kept short, and strong correlations among them should be avoided to maintain low degeneracy. In addition, some descriptors might be chosen to avoid overlap dependence from the local structure of the binding site. In many cases, ML techniques can identify internal correlations among variables and reduce the number of descriptors. One such technique is the principal component analysis (PCA), a nonparametric statistical ML technique, which is a very important hypothesis-generating tool for approaching drug discovery by a systemic perspective.[30] The contributions of hydrophobic and aromatic character, and molecular shape to the chemical similarity of small molecule drugs have been analyzed using PCA to establish correlations of drug structures with transcriptional responses.[31] The PCA approach has been employed to reduce the dimensions of molecular descriptors and construct a chemical space to choose descriptors which are presumed best for the particular effect.[32,33] The molecular property distributions of drug-like, non-drug-like, and natural compounds from traditional Chinese medicines are analyzed by PCA.[34]

The translation of molecular properties into a format suitable for ML requires an appropriate combination of structural and chemical data. However, while the geometrical properties are well represented by molecular fingerprints, the characterization of chemical reactivity is a more challenging task. Typically, the number of hydrogen bond donors or acceptors, functional groups, heats of formation, polarizabilities, etc. are used to represent the ligand’s ability to interact with a protein’s active site. However, in many cases, the chemical reactivity information is incomplete and/or overlaps with the molecular structure, making it insufficient to thoroughly represent the interaction.[27] As a possible alternative to the predictive methods based on the group contribution approach, the conductor-like screening model for real solvents (COSMO-RS) type models, originally developed by Dr. A. Klamt,[35] has attracted substantial interest in the recent years. For the first time, the efficiency of COSMO-RS σ-profiles for the detection of new bioisosteric drug candidates has been studied by Thormann and Klamt.[36] Our idea is to further combine the ability of σ-profiles to represent the chemical potentials of different compounds in almost arbitrary phases with the strength of ML in order to enhance σ-profiles’ drug screening and design capabilities.

The COSMO-RS is a fast and efficient tool for molecular similarity screening and the study of physicochemical and physiological properties originating from solvation thermodynamics and computational quantum mechanics.[37] These methods rely on solvent-accessible surfaces with induced charges screening the electrostatic potential produced by nuclear charges and electronic density. The resulting surfaces are specific to each molecule and projected further to the one-dimensional (1D) function called the σ-profile, which represents the probability distribution of a molecular surface segment having a specific charge density.[35] The COSMO-RS σ-profiles capture in a 1D format the entire set of chemical properties calculated from quantum chemistry, including hydrogen bond donor and acceptor, and hydrophobic (lipophilic) interactions.[38] The σ-profiles have been employed to calculate phase diagrams, azeotropes, and the solubility of drugs and pesticides in water.[38] The calculation of drug absorption, distribution, metabolism, and excretion properties has been demonstrated for high-throughput screening by using the molecular fragment-based approach COSMOfrag.[37] Despite this potential, COSMO-RS σ-profiles have rarely been employed for drug similarity screening, and not at all for SARS-CoV-2 inhibitor screening to date.

In this study, we propose a novel set of drug screening descriptors based on COSMO-RS σ-profiles, augmented by dipole moment and induced charge of the phosphorus atom, to evaluate the chemical similarity of the drugs with nucleotides, as RNA replication transcription initiation activators. The COSMO-RS-based descriptors are more suitable for ML because these effectively account for the chemistry, as expressed by chemical thermodynamics and interaction preferences, and are not strongly correlated. These advantages are highlighted in comparison with standard molecular structure-based descriptors.

We illustrate the proposed approach, using a series of 18 active compounds of market-available drugs built to mimic the chemical structures of adenosine, guanosine, uridine, and cytidine monophosphates (AMP, GMP, UMP, and CMP) that bind to messenger RNA (mRNA), when catalyzed in the central region of RdRp. The drug series includes the pyrimidine and purine nucleoside monophosphate analogues in the form that inhibits mRNA replication from pubchem.ncbi.nlm.nih.gov and recent publications (Supporting Information Table S1). The 64 most stable conformations of these drugs and nucleotides are determined using BIOVIA Discovery Studio 2020.[25] The geometries of these conformations, generated using the Discovery Studio 2020 Visualizer, are optimized using Becke’s three-parameter hybrid functional with the local term of Lee, Yang, and Parr (B3LYP),[39] and double-numerical basis set with polarization functions on H atoms (DNP)[40] in DMol3. The COSMO-RS calculations,[35,41] parametrized to optimally reproduce thermodynamic properties,[42] are conducted by using DMol3, as implemented in BIOVIA Discovery Studio 2020.[25]

To identify the most important hydrogen bond donor and acceptor regions in the σ-profile, the geometry optimization of remdesivir monophosphate interacting with RdRp (PDB ID: 7BV2) is conducted using quantum mechanics/molecular mechanics (QM/MM), as implemented in DMol3 in BIOVIA Discovery Studio 2020. The B3LYP/DNP[39,40] approach and CHARMm36 force field with Momany-Rone charges[43] are applied for the QM (F86101, ALA688, U20, U10, A11) and MM system (remaining), respectively. The QM/MM interactions are treated using electronic embedding with dispersed boundary charges. The molecular docking analysis is conducted using the site-features docking algorithm LibDock[44] implemented in BIOVIA Discovery Studio 2020.

In Figure a, we present the QM/MM optimized model system of the leading drug remdesivir in its monophosphate form integrated into the primer strand mRNA interacting with RdRp (Protein Data Bank entry 7BV2).[7] The hydrogen bonding network with the target U10 is shown in dashed green lines pointing up from the remdesivir monophosphate molecule and with the pattern U20 pointing down from the remdesivir monophosphate molecule. The remdesivir monophosphate-RdRp interactions also include attractive van der Waals, π–cation, π–π stacking, π–alkyl, and covalent as well as repulsive electrostatic interactions with RdRp functional groups, as shown in Figure b. The COSMO-RS σ-profile of remdesivir monophosphate (Figure c) relates to the hydrogen bond donor (large negative values) and acceptor (large positive values) as well as hydrophobic interactions (small negative and positive values). The target U10, pattern U20, and Ala688 hydrogen bonds are highlighted. The hydrogen bond densities (Table ) indicate the donor and acceptor types as negative and positive values, respectively, with absolute values corresponding to the hydrogen bonding strength. The 1D σ-profiles are projected to the scalar values by weighted integration into four donor/acceptor regions. The weights and region positions are selected based on the positions of hydrogen bonds between the drug and RdRp nucleosides (Figure c).

Figure 1

(a) Noncovalent bonding of remdesivir (ball and stick) to mRNA (helix) in RdRp (green ribbon) optimized using QM/MM. Green lines are the hydrogen bonds with U10 (target, two upper bonds) and U20 (pattern, two lower bonds). (b) 2D noncovalent bonding network of remdesivir bonded to mRNA. (c) COSMO-RS σ-profile of remdesivir and its projections W1–W4, with the hydrogen bonds highlighted.

Table 1

Charge Density (e/Å2) for the Hydrogen Bonds Highlighted in Figure c

Nrmsd-ALA688	Nrmsd-T:U10	Ormsd-P:U20	Hrmsd-T:U10	Hrmsd-T:U20
0.0092	0.0130	0.0062	–0.0140	–0.0159

The calculated σ-profile projections (W1–W4) augmented by the COSMO-RS solvent accessible surface, COSMO-RS solvation free energy, phosphorus atom induced charge, and dipole moment constitute descriptor set I. The COSMO-RS solvation free energy is calculated by a statistical thermodynamic treatment of the pairwise interactions of the solvent surface segments.[42] The induced Mulliken charge of the phosphorus atom is selected as a precursor of the chemical activity of the phosphate group toward covalent bonding to mRNA.

The performance of the descriptor set I, made up of eight descriptors, is compared to descriptor set II, made up of all ten drug screening descriptors (nine for molecular structure and one for hydrophobicity, listed in Table ) implemented in BIOVIA Discovery Studio 2020. To select the most important combinations of descriptors for the series of 22 nucleotides and active compounds of nucleotide drug analogues, we employed the PCA implemented in MATLAB[45] to the descriptor sets I and II, without any preprocessing. The PCA linear combination coefficients listed in the Table give the loadings of each of the four principal components (PC). Of these, the most important PC1–PC3 for each descriptor set are presented in the score plot in Figure .

Table 2

Principal Component (PC) Analysis Results for Descriptor Sets I and II

parameter	PC1	PC2	PC3	PC4
Descriptor Set I
COSMO-RS Solvent accessible surface (Å2)	–0.286	0.242	0.567	–0.400
Charge on P (|e|)	–0.265	0.453	–0.138	–0.154
COSMO-RS solvation energy (kcal/mol)	0.522	0.232	–0.048	0.031
Dipole moment (Debye)	–0.393	–0.281	–0.111	–0.544
COSMO-RS σ-profile projection W1 (Å2)	–0.488	0.212	–0.200	0.253
COSMO-RS σ-profile projection W2 (Å2)	–0.035	–0.396	0.694	0.256
COSMO-RS σ-profile projection W3 (Å2)	–0.075	0.595	0.323	0.294
COSMO-RS σ-profile projection W4 (Å2)	–0.420	–0.220	–0.137	0.550
Latency	2.887	1.751	1.190	0.902
Descriptor Set II
AlogP	–0.361	0.232	0.333	–0.010
Molecular weight (g/mol)	0.272	0.233	0.526	0.242
Number of hydrogen bond donors	0.486	0.021	–0.008	0.359
Number of hydrogen bond acceptors	0.446	0.113	0.070	–0.164
Number of rotatable bonds	–0.199	0.238	–0.382	0.803
Number of rings	0.095	0.459	0.381	0.024
Number of aromatic rings	–0.048	0.540	–0.295	–0.246
Molecular fractional polar surface area (Å2)	0.508	–0.102	–0.369	–0.127
Molecular surface area (Å2)	–0.048	0.540	–0.295	–0.246
Molecular polar surface area (Å2)	0.470	0.135	–0.080	0.055
Latency	0.430	0.280	0.133	0.062

Figure 2

(a) Principal components (PC) 1–3 for descriptor set I (Table ) with labels of the COSMO-RS-derived projections W1–W4 of the σ-profile, COSMO-RS solvent accessible surface (SAS), ECOSMO-RS, P atom induced charge (QP), and dipole moment (DM). (b) PC 1–3 for descriptor set II (Table ). (green squares, remdesivir; green diamonds, EIDD-2801; green circles, AMP; yellow squares, GMP; yello triangles, IDX-184; yello circles, mizobirine; blue squares, ribavirin, penciclovir, acyclovir; red circles, remaining drugs and nucleotides.)

From the score plot (Figure a), it is clear that the reference group of conformers (marked by green color) represented by AMP and the active monophosphate forms of the drugs remdesivir and EIDD-2801 known to exhibit some level of efficiency in the inhibition of SARS-CoV-2 mRNA replication[10,11] form a cluster within the calculated PCA. The similar but sparser distribution could be seen for descriptor set II (Figure b). This shows that the similarity of these compounds is clear and could be captured with both sets. However, on a closer look, one can see that for descriptors set II the conformers for each component are closer together than for the proposed new descriptor set I. This could be explained by the structural effects and degeneracy of descriptor set II, as the molecular weight, number of hydrogen bonds, number of donors and acceptors, number of rotatable bonds, number of rings, and number of aromatic rings (listed in Table ) are (nearly) identical for different conformers of the same compound. Another important feature of the proposed descriptor set I is that the other drugs, IDX-184 and mizobirine in monophosphate forms (colored in yellow, Figure a), which have been identified to potentially inhibit the SARS-CoV-2 mRNA replication,[46,47] are better correlated with the reference conformers than those of descriptor set II (Figure b). In a yield reduction assay, mizoribine has inhibited the replication of SARS-CoV-1 more strongly than ribavirin, reducing the virus titers to one-tenth or less.[46] The anti-hepatitis C virus replication inhibitor drug IRX-184 has shown promising results as SARS-CoV-2 RdRp inhibitor based on molecular docking that is yet to be confirmed experimentally.[47] Moreover, the extent of correlation among the conformers of the same drug is enhanced when descriptor set I is used (Figure a). The monophosphate forms of the drugs ribavirin, penciclovir, and acyclovir, shown in blue squares, are also anticipated to have some RdRp replication inhibition activity. The herpesvirus RNA and deoxyribonucleic acid (DNA) replication inhibitors ribavirin and penciclovir have been shown to reduce the SARS-CoV-2 viral infection against a clinical isolate in vitro at high concentrations.[20] Acyclovir treatment has been associated with significantly extending the life and improving the circulatory and pulmonary function in patients with ventilator-associated pneumonia and high levels of herpesvirus replication,[48] making it promising for the treatment of COVID-19 symptoms. Currently, proven drug efficiency results are scarce, making further descriptor set evaluation challenging. This analysis could be extended in the future when more results become available in the literature.

The PC loadings listed in Table represent the normalized coefficients of linear combinations for transformation properties represented by descriptors to the new PCs. The PC loadings also indicate the weight of the descriptors, while the PC latency indicates its strength or order of importance of PCA components. For descriptor set I, the COSMO-RS solvation energy and σ-profile projections W1 and W4 have the largest weights in PC1 (with the highest latency), followed by the dipole moment, accessible surface, and phosphorus atom charge. The largest weights in PC2 are calculated for the σ-profile projections W2 and W3, and the phosphorus atom charge. In PC3, the largest weights are for σ-profile projection W2 and the solvent accessible surface. For descriptor set II, the molecular polar surface area and fractional polar surface area as well as the number of hydrogen bond donors and acceptors and hydrophobicity are the most important for PC1, based on the weights. In PC2, the molecular surface area and number of aromatic rings have the largest weights. The large weights of the molecular surface area, hydrogen bonding, and hydrophobicity descriptors of set II effectively justify the use of the more advanced descriptors of the proposed COSMO-RS-based set I. The polar surface area descriptors of set II have close physical meaning to the COSMO-RS σ-profiles, as both describe the “polarization” of the molecular surface. However, the advantage of the proposed descriptors is that the COSMO-RS σ-profiles are based on more accurate DFT calculations, while polar surface area descriptors are calculated from molecular mechanics force fields and cannot fully account for the correlation of atoms within a molecule. This advantage is particularly important to separate the different conformers that are treated as distinct molecules with calculated atomic charges and COSMO-RS σ-profiles rather than with the same set of force field parameters using descriptor set II.

The superior performance of proposed descriptor set I exhibited by the previously noted enhanced clustering and correlation with the reference conformers is attributed to the statistical thermodynamics the COSMO-RS chemical potential accounts for. This effectively highlights the importance of the COSMO-RS σ-profiles that are calculated by integration over the molecular surface, account for hydrogen bonding and hydrophobic interactions, and are parametrized to reproduce chemical thermodynamic properties.[42]

In order to enhance comparison and compensate for the lack of clinical trial results, we also perform LibDock analysis for the 22 pyrimidine and purine nucleotides and analogue drugs. The site-features docking algorithm LibDock performs as well as docking programs based on genetic/growing and Monte Carlo driven algorithms.[49] Of course, a good docking score does not necessarily mean that the ligand is a good inhibitor, but only a good candidate for potential drug selection. In Figure , we present the drug docking results, where 11 drugs and 3 nucleotides (shown as Group 1 in Figure a) dock around an active site of RdRp and have high LibDock scores (Figure b). The results also show some agreement with the PCA plot (Figure ). Penciclovir, ribavirin, and the nucleotide GMP have the highest LibDock scores, followed by acyclovir and the nucleotide CMP. Among these, only the nucleotide GMP is correlated well with the reference conformers including AMP, based on the PCA in Figure a, but the LibDock score of AMP is low. The drugs ribavirin, penciclovir, and acyclovir, anticipated to have some inhibition activity based on the LibDock score, are weakly correlated with the reference conformers group, based on the PCA (Figure a). Clearly, the difference between PCA and docking results could be attributed to the chemical thermodynamics properties accounted for by descriptor set I. The ML results include additional information not only about complex processes but also about the interactions in the real system.

Figure 3

(a) Docking of the drugs and nucleotides on the mRNA of RdRp, showing the localization of the drugs of Group 1. (b) LibDock score of the 22 pyrimidine (Pyr) and purine nucleotides and drug analogues.

A novel set of descriptors based on COSMO-RS σ-profiles and chemical thermodynamics is proposed and evaluated using PCA for the initial screening of a series of nucleotides and nucleotide-analog RdRp replication inhibitor drugs to help accelerate the discovery of COVID-19 treatments. The proposed descriptor set that accounts for chemical thermodynamics in addition to molecular structure and features descriptors with low degeneracy is evaluated with respect to a commonly used descriptor set based on molecular structure alone and a molecular docking algorithm. The PCA results show that the novel σ-profile-based descriptor set I clearly correlates the leading COVID-19 drugs remdesivir and EIDD-2801 in monophosphate forms and highlights weaker correlations with drugs that have been reported to exhibit anti-SARS-CoV-2 activity. The comparison shows that the proposed descriptor set is superior to both the commonly used descriptor set II and molecular docking scores.

The presented approach for conducting drug screening by using COSMO-RS-based descriptors and ML can be further augmented by taking into account (pre)clinical trial result-based descriptors. In addition, the proposed set of descriptors can be further enhanced based on the target, and combinations of antiviral drugs can be evaluated and optimized to continue improving COVID-19 treatments.