# Discovery and Validation of Traditional Chinese and Western Medicine Combination Antirheumatoid Arthritis Drugs Based on Machine Learning (Random Forest Model)

Authors: Jijia Sun; Qinghua Ni; Fengyan Jiang; Baocheng Liu; Jianying Wang; Lei Zhang; Jihan Huang
Journal: BioMed Research International
 (2023)

---

## Abstract

The combination of traditional Chinese medicine (TCM) and Western medicine is a promising method for treating rheumatoid arthritis (RA). Combining the two fully exploits the advantages of Western and TCM to treat RA and has the potential to greatly improve the therapeutic effect on RA. In this study, we developed a combination drug training set by using 16 characteristic variables based on the characteristics of small molecules of TCM ingredients and Food and Drug Administration-certified combination drug data downloaded from the DrugCombDB database. Furthermore, we compared the prediction and classification abilities of five models: thek-nearest neighbors, naive Bayes, support vector machine, random forest, and AdaBoost algorithms. The random forest model was selected as the classification and prediction model for Western and TCM and Western combination drugs. We collected data for 41 small molecules of TCM ingredients from the Traditional Chinese Medicine Systems Pharmacology database and 10 small molecule drugs commonly used in anti-RA treatment from the DrugBank database. Combinations of Western and TCM for anti-RA treatment were screened. Finally, the CellTiter-Glo method was used to determine the synergy of these combinations, and the 15 most predicted drug combinations were carried out experimental verification. Myricetin, rhein, nobiletin, and fisetin had high synergy with celecoxib, and rhein had high synergy with hydroxychloroquine. The preliminary findings of this study can be further applied for practical clinical anti-RA combined treatment strategies and serve as a reference for clinical treatment of RA with integrated Western and TCM.


---

## Body

## 1. Introduction

Modern medical research has demonstrated that many diseases often present complex pathogenic mechanisms, and consequently, the past methods for treating diseases with drugs are no longer adequate. In recent years, drug combination has been reported to offer clear advantages in treating tumors, AIDS, cardiovascular diseases, osteoarthritis, and other complex diseases; it has been widely implemented in clinical practice [1–4], with treatment methods based on combined drugs demonstrating advantageous synergistic effects, or a “1+1≥2” effect. The therapeutic effects of two or more drug combinations are greater than those achieved when each drug is used alone. Therefore, developing and evaluating potentially synergetic combination drugs, whether for development of new drugs or clinical treatment optimization, are essential. However, most existing combination drugs have been developed on the basis of limited clinical experience or experimental strategies. Verified combination drugs with synergistic therapeutic effects currently comprise a small portion of listed drugs. Therefore, making full use of knowledge of existing drug combinations, discovering new drug combinations based on well-accepted and understood drugs, and establishing a model for evaluating the synergistic effects of drug combinations are key concerns that must be addressed.Rheumatoid arthritis (RA) is a chronic, systemic, inflammatory disease of the synovial and other joints. RA is a chronic progressive autoimmune disease [5], and its underlying pathogenesis is currently unclear. However, studies have reported that the incidence of RA increases with age due to infection and inflammatory mediators, and approximately 0.3%–1.0% of the global population is affected by the disease every year [6, 7]. At present, the drugs used in clinical treatment of RA are mainly nonsteroidal anti-inflammatory drugs (NSAIDs), disease modifying antirheumatic drugs (DMARDs), and biological macromolecular therapy. Because RA is complex, several studies have explored different combination therapy strategies for its treatment [8–10]. However, some therapeutic drugs, such as methotrexate, may produce unfavorable side effects [11, 12]. Therefore, if sole use of anti-RA combination therapy has the potential to produce unfavorable side effects, long-term use in patients with RA may do greater harm. Traditional Chinese medicine (TCM) has long been used to treat RA, and several studies have found that combined anti-RA TCM and a compound prescription have not only anti-inflammatory, analgesic, immune regulation, and multilevel and multilink therapeutic effects but also the advantages of high safety, few adverse reactions, and low cost. Such a combination has gained attention in studies of RA treatment and has increasingly attracted international attention [13, 14].Therefore, we propose that the combination of Western medicine and TCM may be a promising method for improving the efficacy of anti-RA treatments and reducing side effects. Through their combination, the advantages of both Western medicine and TCM can be fully exploited to improve the therapeutic effects of RA treatment. However, because numerous active ingredients in TCM and hundreds of listed Western medicines have been reported to be effective in anti-RA treatments, the potential drug combinations are in the thousands or even tens of thousands. This multitude of combinations cannot feasibly be evaluated or verified through clinical practice and trials; therefore, the most practical method for identifying potential combinations would be to use drug data for existing combination drugs to model, evaluate, and screen potential anti-RA combinations of Western medicine and TCM and then experimentally validate these models. Therefore, in this study, we developed a classification and prediction model for combining Western and TCM drugs based on characteristics obtained through small molecule research on TCM and existing data on combined drugs. This model was developed using machine learning modeling followed by manual screening of anti-RA combination drugs. Through experimental verification, we identified effective combinations of Western and TCM drugs for treating RA; these findings can serve as a reference for clinical treatment of RA by using such combined treatment. The flow chart illustrating the process of drug discovery and verification for anti-RA integrative medicine through machine learning modeling is displayed in Figure1.Figure 1
Flow chart of drug discovery and verification for anti-RA integrative medicine through machine learning modeling.


## 2. Materials and Methods

### 2.1. Data Collection and Collation

The DrugCombDB database (http://drugcombdb.denglab.org/) is a comprehensive combinatorial drug database [15] containing integrated drug combinations from various data sources. To ensure the reliability of the research data, we first downloaded the Food and Drug Administration- (FDA-) provided known drug combination dataset from the database (up to May 31, 2019). The molecular formula, Chemical Abstracts Service (CAS) number, canonical simplified molecular-input line-entry system (SMILES), and other chemical informatics data corresponding to each combination of drugs and their action targets were collected from the DrugBank database (https://go.drugbank.com/).Furthermore, we searched “Rheumatoid Arthritis” on the Traditional Chinese Medicine Systems Pharmacology (TCMSP) database (https://tcmspw.com/) to retrieve target Chinese medicine ingredient data reported to be relevant to treating RA. We obtained the main small molecules chemical information from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/). In addition, several common anti-RA drugs were selected from the DrugBank database along with those discussed in the literature and those recommended by clinicians. All small molecule data for the TCM ingredients (collected from the TCMSP) and for known anti-RA Western drugs were combined in pairs to form the final prediction sample set.


### 2.2. Construction of Combined Drug Characteristic Variables

#### 2.2.1. Drug Combination Features Based on Molecular Fingerprint Similarity

Canonical SMILES of all drugs were collected from the DrugBank and PubChem databases, and the 1024-dimensional molecular fingerprints for each of the small molecule drugs were calculated using DRAGON software (version 7.0). The fingerprint similarity (SFP) between two small molecule drugs was then calculated using Tanimoto coefficient-based similarity:
(1)SFPdA,dB=A∙BA2+B2−A∙B,where A is the 1024-dimensional molecular fingerprint feature vector corresponding to the small molecule drug dA, A is the vector length, B is the 1024-dimensional molecular fingerprint feature vector corresponding to the small molecule drug dB, B is the vector length, and A∙B is the vector inner product.


#### 2.2.2. Drug Combination Characteristics Based on Small Molecule ADME Similarity

To predict the adsorption, distribution, metabolism, and excretion (ADME) characteristics of each drug, Canonical SMILES for all small molecule drugs were imported into the SwissADME online system (http://www.swissadme.ch/) [16]. We selected 12 indicators as characteristic indexes for our ADME similarity calculations: molecular weight, heavy atoms, aromatic heavy atoms, fraction Csp3, rotatable bonds, H-bond acceptors, H-bond donors, molecular refractivity, topological polar surface area, Silicos-IT Log P, and Silicos-IT LogSw. All ADME indexes were first normalized, and the Euclidean distance was then used to calculate the similarity value SADME of ADME:
(2)SADME=∑i=1nxiA−xiB21/2,where xiA is the ADME index corresponding to the drug small molecule dA and xiB is the ADME index corresponding to the drug small molecule dB (n=12).


#### 2.2.3. Drug Combination Characteristics Based on Sequence Similarity of Small Molecular Targets

The action targets of all drugs were obtained from DrugBank and the TCMSP, and the protr package (version 1.6.2) based on R language was used to obtain the sequence of action protein targets of all small molecule drugs. The sequence similarity between the two target proteins,SSeqdA,dB, was calculated using the Smith–Waterman algorithm. The maximum, minimum, median, and mean values of the sequence similarity SSeqdA,dB between all target proteins of the dAanddB were included in the characteristic variables as SSeqmax,SSeqmin,SSeqmed,andSSeqmean, respectively.


#### 2.2.4. Drug Combination Characteristics Based on GO Functional Similarity of Small Molecular Targets

To calculate the gene ontology (GO) functional similaritySGOdA,dB of the two drugs dAanddB, we used the GOSemSim package (version 2.14.2). Similarly, the minimum, maximum, median, and mean values of the GO functional similarity SGOdA,dB of the two drugs dAanddB were included in the characteristic variables as SGOmin,SGOmax,SGOmed,andSGOmean, respectively.


#### 2.2.5. Drug Combination Characteristics Based on Pathway Similarity of Small Molecule Targets

We used the clusterProfiler package (version 3.14.3) to perform a Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway enrichment analysis on all drug small molecule targets. We then screened the enrichment pathways of all targets based onp value < 0.05 and q value < 0.05. A vector, P=p1,p2,⋯,pn, consisting of all pathways was constructed according to the enriched pathways, pkk=1,2,⋯,n (where n is the number of pathways). If a target was enriched in a pathway pk in the target set TA of the drug small molecule dA, then the value was set as pkA=1 at the corresponding position of the drug small molecule A pathway vector PA, and if not, it was set as 0:
(3)pkA=1,TA∩pk≠∅,0,others.The target pathway vectorPB of the drug small molecule dB can be similarly obtained using a string similarity formula to calculate the pathway similarity value (SPathway) of two drugs:
(4)SPathwaydA,dB=PA∙PBPA×PB,where PA and PB are path vector lengths and PA·PB is the vector inner product.


#### 2.2.6. Drug Combination Characteristics Based on Distance of Small Molecule Targets in Human PPI Network

The distance and network proximity of small molecule targets of different combinations of drugs on the human protein–protein interaction (PPI) network were calculated to measure the synergy characteristics of combination drugs within the PPI network:(5)SPPImindA,dB=11+e−minda,b,SPPImaxdA,dB=11+e−maxda,b,SPPImeddA,dB=11+e−medda,b,SPPImeandA,dB=11+e−dAB,SnpdA,dB=SPPImeandA,dB−SPPImeandA,dA+SPPImeandB,dB2,whereda,b is the shortest path distance between two protein targets, aandb, within the PPI network (calculated using the igraph package [version 1.2.6]), dAB=1/A×B∑a∈A,b∈Bda,b, A is the number of targets for the drug small molecule dA, and B is the number of targets for the drug small molecule dB. The human PPI background network used in the calculation was derived from the literature [17]; the PPI network contained a total of 16677 proteins and 243603 interactions.The Tanimoto coefficient-based similarity, Euclidean distance, and cosine similarity calculations used in this paper were calculated using the philentropy package (version 0.4.0) based on R language (version 4.0.2).


### 2.3. Construction and Selection of Combination Drug Prediction Model

To comprehensively determine and select an appropriate machine learning model for modeling and prediction in this study, we selected different machine learning classifiers for use with the training set and the verification set and compared their classification accuracies. In this study, we selected and compared five classical machine learning models: thek-nearest neighbors (k-NN) [18], naive Bayes (NB) [19], support vector machine (SVM) [20], random forest (RF) [21], and AdaBoost classifiers [22]. The caret package (version 6.0-86), the klaR package (version 0.6-15), the svmRadial model in the kernlab package (version 0.9-29), the randomForest package (version 4.6-14), and the ada package (version 2.0-5) were used for the k-NN, naive Bayes, support vector machine, random forest, and AdaBoost classifiers, respectively.Furthermore, we used established sample training sets to compare the classification prediction effects of the five machine learning models. First, the sample set was divided into prediction training sets and internal validation sets at a 7 : 3 ratio. Then, the prediction training set was input into the five machine learning models, and the machine learning models were trained using three sets of 10-fold cross-validation. The trained machine learning models were evaluated using the internal validation sets, and the model with the best classification prediction effect was selected as the final drug prediction model for the Western and TCM drug combinations.


### 2.4. Evaluation of Machine Learning Models

In this study, we defined precision rate as the proportion of positive samples correctly predicted by the model out of all positive predictions, recall rate as the proportion of positive samples correctly predicted by the model out of all predictions, accuracy as the proportion of positive and negative samples to all samples correctly predicted by the model, andF1-score as the overall score of the model when precision rate and recall rate were given equal weight (calculated using Formulas (6)–(9), respectively). A higher score indicated a better classification effect.
(6)Precision=TPTP+FP,(7)Recall=TPTP+FN,(8)Accuracy=TP+TNTP+FN+FP+TN,(9)F1‐score=2×Precision×RecallPrecision+Recall,where true positive (TP) indicates a sample is positive and was predicted to be positive, false negative (FN) indicates a sample is positive but was predicted to be negative, false positive (FP) indicates a sample is negative but was predicted to be positive, and true negative (TN) indicates a sample is negative and was predicted to be negative.In the plot of the receiver operating characteristic (ROC) curve, the horizontal axis represents the false positive rate (FPR) and the vertical axis represents the true positive rate (TPR). The area under the ROC curve is the area surrounded by the ROC and its lower coordinate axis and is represented as AUROC in this paper. In addition, because the precision and recall rates of the model are not fixed, when the classification threshold of the model is adjusted, the precision and recall rates continue to change. The precision–recall curve represents the relationship between the precision rate of the ordinate and the recall rate of the abscissa with the change of the classification thresholdp value; the area beneath this curve is represented as the area under the precision–recall curve (AUPRC).


### 2.5.In Vitro Cell Experiment Verification

#### 2.5.1. Cell Culture and Reagents

RAW264.7 cells were purchased from the American Type Culture Collection. The culture reagents included Dulbecco’s modified Eagle medium (DMEM; Gibco, C11995500CP), Roswell Park Memorial Institute (RPMI) 1640 (Gibco, C11875500BT), fetal bovine serum (Bio IND, 04-002-1A), antimycotic (Lifetechnologies, 15240-112), phosphate-buffered saline (PBS), pH 7.4 (Gibco, 10010-023), trypsin-EDTA (0.05%; Lifetechnologies, 25300-054), bovine serum albumin (Lifetechnologies, 15561012), and CellTiter-Glo (CTG; PROMEGA, G7572). Consumables and instruments included a cell culture plate (Corning), cell culture flask (Corning), microplate tester (BioTek, HM-1), and conventional instruments, such as a CO2 incubator (Thermo 3111) and biosafety cabinet (Heal Force, HFSAFe-1200LC).


#### 2.5.2. Experimental Procedure and Methods

First, the different compounds were diluted with DMEM medium to the desired concentrations of the experimental design, namely, 0, 0.032, 0.16, 0.8, 4, 20, and 100μM. The cells were then counted and inoculated into 96-well culture plates at 2000 cells/100 μL/well; 90 μL of culture medium (including serum) was added to each well for overnight culture, and 10 μL of drug solution was added. After a 48 h culture, CTG was used to detect the proliferation of RAW264.7 cells; CTG (100 μL) was added to each well, the content of which was incubated at room temperature in a dark environment for 10 min. The chemiluminescence value of each well was measured at 500 ms using enzyme calibration.


#### 2.5.3. Combination Index Measurement

In this study, according to the principle of the median effect—specifically by employing the Chou–Talalay model [1]—the synergistic effect between each set of drugs was analyzed, and the effects of integrated Western and TCM drugs on RA cell lines were evaluated using the combination index (CI). The CI indicates the degree of drug synergy, addition, or antagonism for any given drug combination:
(10)CI=DADAx+DBDBx,where when CI=1, the two drugs have an additive effect, whenCI<1, the two drugs have a synergistic effect, and when CI>1, the two drugs are antagonistic. DAandDB represent the concentrations of the two drugs, dAanddB, when administered in combination, and DAxandDBx represent the concentrations of the two drugs, dAanddB, when they have been administered alone and the combined drug inhibition rate has been reached.
(11)Dx=Dmfafu1/m.In Formula (11), faandfu are the inhibition and survival rates, respectively, when the two drugs are combined, fa=1−fu, Dm is the median dose, and m is the curve coefficient (calculated according to half-maximal inhibitory concentration (IC50) theory).In this study, CompuSyn software (http://www.combosyn.com/) was used to calculate the CI values of the drug combinations to evaluate whether each screened combination of Western and TCM drugs had synergistic therapeutic effects.


### 2.6. Statistical Analysis

Data were expressed as themean±standarddeviation, and a t-test and Mann–Whitney nonparametric tests were performed using IBM SPSS Statistics (version 26.0).


## 2.1. Data Collection and Collation

The DrugCombDB database (http://drugcombdb.denglab.org/) is a comprehensive combinatorial drug database [15] containing integrated drug combinations from various data sources. To ensure the reliability of the research data, we first downloaded the Food and Drug Administration- (FDA-) provided known drug combination dataset from the database (up to May 31, 2019). The molecular formula, Chemical Abstracts Service (CAS) number, canonical simplified molecular-input line-entry system (SMILES), and other chemical informatics data corresponding to each combination of drugs and their action targets were collected from the DrugBank database (https://go.drugbank.com/).Furthermore, we searched “Rheumatoid Arthritis” on the Traditional Chinese Medicine Systems Pharmacology (TCMSP) database (https://tcmspw.com/) to retrieve target Chinese medicine ingredient data reported to be relevant to treating RA. We obtained the main small molecules chemical information from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/). In addition, several common anti-RA drugs were selected from the DrugBank database along with those discussed in the literature and those recommended by clinicians. All small molecule data for the TCM ingredients (collected from the TCMSP) and for known anti-RA Western drugs were combined in pairs to form the final prediction sample set.


## 2.2. Construction of Combined Drug Characteristic Variables

### 2.2.1. Drug Combination Features Based on Molecular Fingerprint Similarity

Canonical SMILES of all drugs were collected from the DrugBank and PubChem databases, and the 1024-dimensional molecular fingerprints for each of the small molecule drugs were calculated using DRAGON software (version 7.0). The fingerprint similarity (SFP) between two small molecule drugs was then calculated using Tanimoto coefficient-based similarity:
(1)SFPdA,dB=A∙BA2+B2−A∙B,where A is the 1024-dimensional molecular fingerprint feature vector corresponding to the small molecule drug dA, A is the vector length, B is the 1024-dimensional molecular fingerprint feature vector corresponding to the small molecule drug dB, B is the vector length, and A∙B is the vector inner product.


### 2.2.2. Drug Combination Characteristics Based on Small Molecule ADME Similarity

To predict the adsorption, distribution, metabolism, and excretion (ADME) characteristics of each drug, Canonical SMILES for all small molecule drugs were imported into the SwissADME online system (http://www.swissadme.ch/) [16]. We selected 12 indicators as characteristic indexes for our ADME similarity calculations: molecular weight, heavy atoms, aromatic heavy atoms, fraction Csp3, rotatable bonds, H-bond acceptors, H-bond donors, molecular refractivity, topological polar surface area, Silicos-IT Log P, and Silicos-IT LogSw. All ADME indexes were first normalized, and the Euclidean distance was then used to calculate the similarity value SADME of ADME:
(2)SADME=∑i=1nxiA−xiB21/2,where xiA is the ADME index corresponding to the drug small molecule dA and xiB is the ADME index corresponding to the drug small molecule dB (n=12).


### 2.2.3. Drug Combination Characteristics Based on Sequence Similarity of Small Molecular Targets

The action targets of all drugs were obtained from DrugBank and the TCMSP, and the protr package (version 1.6.2) based on R language was used to obtain the sequence of action protein targets of all small molecule drugs. The sequence similarity between the two target proteins,SSeqdA,dB, was calculated using the Smith–Waterman algorithm. The maximum, minimum, median, and mean values of the sequence similarity SSeqdA,dB between all target proteins of the dAanddB were included in the characteristic variables as SSeqmax,SSeqmin,SSeqmed,andSSeqmean, respectively.


### 2.2.4. Drug Combination Characteristics Based on GO Functional Similarity of Small Molecular Targets

To calculate the gene ontology (GO) functional similaritySGOdA,dB of the two drugs dAanddB, we used the GOSemSim package (version 2.14.2). Similarly, the minimum, maximum, median, and mean values of the GO functional similarity SGOdA,dB of the two drugs dAanddB were included in the characteristic variables as SGOmin,SGOmax,SGOmed,andSGOmean, respectively.


### 2.2.5. Drug Combination Characteristics Based on Pathway Similarity of Small Molecule Targets

We used the clusterProfiler package (version 3.14.3) to perform a Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway enrichment analysis on all drug small molecule targets. We then screened the enrichment pathways of all targets based onp value < 0.05 and q value < 0.05. A vector, P=p1,p2,⋯,pn, consisting of all pathways was constructed according to the enriched pathways, pkk=1,2,⋯,n (where n is the number of pathways). If a target was enriched in a pathway pk in the target set TA of the drug small molecule dA, then the value was set as pkA=1 at the corresponding position of the drug small molecule A pathway vector PA, and if not, it was set as 0:
(3)pkA=1,TA∩pk≠∅,0,others.The target pathway vectorPB of the drug small molecule dB can be similarly obtained using a string similarity formula to calculate the pathway similarity value (SPathway) of two drugs:
(4)SPathwaydA,dB=PA∙PBPA×PB,where PA and PB are path vector lengths and PA·PB is the vector inner product.


### 2.2.6. Drug Combination Characteristics Based on Distance of Small Molecule Targets in Human PPI Network

The distance and network proximity of small molecule targets of different combinations of drugs on the human protein–protein interaction (PPI) network were calculated to measure the synergy characteristics of combination drugs within the PPI network:(5)SPPImindA,dB=11+e−minda,b,SPPImaxdA,dB=11+e−maxda,b,SPPImeddA,dB=11+e−medda,b,SPPImeandA,dB=11+e−dAB,SnpdA,dB=SPPImeandA,dB−SPPImeandA,dA+SPPImeandB,dB2,whereda,b is the shortest path distance between two protein targets, aandb, within the PPI network (calculated using the igraph package [version 1.2.6]), dAB=1/A×B∑a∈A,b∈Bda,b, A is the number of targets for the drug small molecule dA, and B is the number of targets for the drug small molecule dB. The human PPI background network used in the calculation was derived from the literature [17]; the PPI network contained a total of 16677 proteins and 243603 interactions.The Tanimoto coefficient-based similarity, Euclidean distance, and cosine similarity calculations used in this paper were calculated using the philentropy package (version 0.4.0) based on R language (version 4.0.2).


## 2.2.1. Drug Combination Features Based on Molecular Fingerprint Similarity

Canonical SMILES of all drugs were collected from the DrugBank and PubChem databases, and the 1024-dimensional molecular fingerprints for each of the small molecule drugs were calculated using DRAGON software (version 7.0). The fingerprint similarity (SFP) between two small molecule drugs was then calculated using Tanimoto coefficient-based similarity:
(1)SFPdA,dB=A∙BA2+B2−A∙B,where A is the 1024-dimensional molecular fingerprint feature vector corresponding to the small molecule drug dA, A is the vector length, B is the 1024-dimensional molecular fingerprint feature vector corresponding to the small molecule drug dB, B is the vector length, and A∙B is the vector inner product.


## 2.2.2. Drug Combination Characteristics Based on Small Molecule ADME Similarity

To predict the adsorption, distribution, metabolism, and excretion (ADME) characteristics of each drug, Canonical SMILES for all small molecule drugs were imported into the SwissADME online system (http://www.swissadme.ch/) [16]. We selected 12 indicators as characteristic indexes for our ADME similarity calculations: molecular weight, heavy atoms, aromatic heavy atoms, fraction Csp3, rotatable bonds, H-bond acceptors, H-bond donors, molecular refractivity, topological polar surface area, Silicos-IT Log P, and Silicos-IT LogSw. All ADME indexes were first normalized, and the Euclidean distance was then used to calculate the similarity value SADME of ADME:
(2)SADME=∑i=1nxiA−xiB21/2,where xiA is the ADME index corresponding to the drug small molecule dA and xiB is the ADME index corresponding to the drug small molecule dB (n=12).


## 2.2.3. Drug Combination Characteristics Based on Sequence Similarity of Small Molecular Targets

The action targets of all drugs were obtained from DrugBank and the TCMSP, and the protr package (version 1.6.2) based on R language was used to obtain the sequence of action protein targets of all small molecule drugs. The sequence similarity between the two target proteins,SSeqdA,dB, was calculated using the Smith–Waterman algorithm. The maximum, minimum, median, and mean values of the sequence similarity SSeqdA,dB between all target proteins of the dAanddB were included in the characteristic variables as SSeqmax,SSeqmin,SSeqmed,andSSeqmean, respectively.


## 2.2.4. Drug Combination Characteristics Based on GO Functional Similarity of Small Molecular Targets

To calculate the gene ontology (GO) functional similaritySGOdA,dB of the two drugs dAanddB, we used the GOSemSim package (version 2.14.2). Similarly, the minimum, maximum, median, and mean values of the GO functional similarity SGOdA,dB of the two drugs dAanddB were included in the characteristic variables as SGOmin,SGOmax,SGOmed,andSGOmean, respectively.


## 2.2.5. Drug Combination Characteristics Based on Pathway Similarity of Small Molecule Targets

We used the clusterProfiler package (version 3.14.3) to perform a Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway enrichment analysis on all drug small molecule targets. We then screened the enrichment pathways of all targets based onp value < 0.05 and q value < 0.05. A vector, P=p1,p2,⋯,pn, consisting of all pathways was constructed according to the enriched pathways, pkk=1,2,⋯,n (where n is the number of pathways). If a target was enriched in a pathway pk in the target set TA of the drug small molecule dA, then the value was set as pkA=1 at the corresponding position of the drug small molecule A pathway vector PA, and if not, it was set as 0:
(3)pkA=1,TA∩pk≠∅,0,others.The target pathway vectorPB of the drug small molecule dB can be similarly obtained using a string similarity formula to calculate the pathway similarity value (SPathway) of two drugs:
(4)SPathwaydA,dB=PA∙PBPA×PB,where PA and PB are path vector lengths and PA·PB is the vector inner product.


## 2.2.6. Drug Combination Characteristics Based on Distance of Small Molecule Targets in Human PPI Network

The distance and network proximity of small molecule targets of different combinations of drugs on the human protein–protein interaction (PPI) network were calculated to measure the synergy characteristics of combination drugs within the PPI network:(5)SPPImindA,dB=11+e−minda,b,SPPImaxdA,dB=11+e−maxda,b,SPPImeddA,dB=11+e−medda,b,SPPImeandA,dB=11+e−dAB,SnpdA,dB=SPPImeandA,dB−SPPImeandA,dA+SPPImeandB,dB2,whereda,b is the shortest path distance between two protein targets, aandb, within the PPI network (calculated using the igraph package [version 1.2.6]), dAB=1/A×B∑a∈A,b∈Bda,b, A is the number of targets for the drug small molecule dA, and B is the number of targets for the drug small molecule dB. The human PPI background network used in the calculation was derived from the literature [17]; the PPI network contained a total of 16677 proteins and 243603 interactions.The Tanimoto coefficient-based similarity, Euclidean distance, and cosine similarity calculations used in this paper were calculated using the philentropy package (version 0.4.0) based on R language (version 4.0.2).


## 2.3. Construction and Selection of Combination Drug Prediction Model

To comprehensively determine and select an appropriate machine learning model for modeling and prediction in this study, we selected different machine learning classifiers for use with the training set and the verification set and compared their classification accuracies. In this study, we selected and compared five classical machine learning models: thek-nearest neighbors (k-NN) [18], naive Bayes (NB) [19], support vector machine (SVM) [20], random forest (RF) [21], and AdaBoost classifiers [22]. The caret package (version 6.0-86), the klaR package (version 0.6-15), the svmRadial model in the kernlab package (version 0.9-29), the randomForest package (version 4.6-14), and the ada package (version 2.0-5) were used for the k-NN, naive Bayes, support vector machine, random forest, and AdaBoost classifiers, respectively.Furthermore, we used established sample training sets to compare the classification prediction effects of the five machine learning models. First, the sample set was divided into prediction training sets and internal validation sets at a 7 : 3 ratio. Then, the prediction training set was input into the five machine learning models, and the machine learning models were trained using three sets of 10-fold cross-validation. The trained machine learning models were evaluated using the internal validation sets, and the model with the best classification prediction effect was selected as the final drug prediction model for the Western and TCM drug combinations.


## 2.4. Evaluation of Machine Learning Models

In this study, we defined precision rate as the proportion of positive samples correctly predicted by the model out of all positive predictions, recall rate as the proportion of positive samples correctly predicted by the model out of all predictions, accuracy as the proportion of positive and negative samples to all samples correctly predicted by the model, andF1-score as the overall score of the model when precision rate and recall rate were given equal weight (calculated using Formulas (6)–(9), respectively). A higher score indicated a better classification effect.
(6)Precision=TPTP+FP,(7)Recall=TPTP+FN,(8)Accuracy=TP+TNTP+FN+FP+TN,(9)F1‐score=2×Precision×RecallPrecision+Recall,where true positive (TP) indicates a sample is positive and was predicted to be positive, false negative (FN) indicates a sample is positive but was predicted to be negative, false positive (FP) indicates a sample is negative but was predicted to be positive, and true negative (TN) indicates a sample is negative and was predicted to be negative.In the plot of the receiver operating characteristic (ROC) curve, the horizontal axis represents the false positive rate (FPR) and the vertical axis represents the true positive rate (TPR). The area under the ROC curve is the area surrounded by the ROC and its lower coordinate axis and is represented as AUROC in this paper. In addition, because the precision and recall rates of the model are not fixed, when the classification threshold of the model is adjusted, the precision and recall rates continue to change. The precision–recall curve represents the relationship between the precision rate of the ordinate and the recall rate of the abscissa with the change of the classification thresholdp value; the area beneath this curve is represented as the area under the precision–recall curve (AUPRC).


## 2.5.In Vitro Cell Experiment Verification

### 2.5.1. Cell Culture and Reagents

RAW264.7 cells were purchased from the American Type Culture Collection. The culture reagents included Dulbecco’s modified Eagle medium (DMEM; Gibco, C11995500CP), Roswell Park Memorial Institute (RPMI) 1640 (Gibco, C11875500BT), fetal bovine serum (Bio IND, 04-002-1A), antimycotic (Lifetechnologies, 15240-112), phosphate-buffered saline (PBS), pH 7.4 (Gibco, 10010-023), trypsin-EDTA (0.05%; Lifetechnologies, 25300-054), bovine serum albumin (Lifetechnologies, 15561012), and CellTiter-Glo (CTG; PROMEGA, G7572). Consumables and instruments included a cell culture plate (Corning), cell culture flask (Corning), microplate tester (BioTek, HM-1), and conventional instruments, such as a CO2 incubator (Thermo 3111) and biosafety cabinet (Heal Force, HFSAFe-1200LC).


### 2.5.2. Experimental Procedure and Methods

First, the different compounds were diluted with DMEM medium to the desired concentrations of the experimental design, namely, 0, 0.032, 0.16, 0.8, 4, 20, and 100μM. The cells were then counted and inoculated into 96-well culture plates at 2000 cells/100 μL/well; 90 μL of culture medium (including serum) was added to each well for overnight culture, and 10 μL of drug solution was added. After a 48 h culture, CTG was used to detect the proliferation of RAW264.7 cells; CTG (100 μL) was added to each well, the content of which was incubated at room temperature in a dark environment for 10 min. The chemiluminescence value of each well was measured at 500 ms using enzyme calibration.


### 2.5.3. Combination Index Measurement

In this study, according to the principle of the median effect—specifically by employing the Chou–Talalay model [1]—the synergistic effect between each set of drugs was analyzed, and the effects of integrated Western and TCM drugs on RA cell lines were evaluated using the combination index (CI). The CI indicates the degree of drug synergy, addition, or antagonism for any given drug combination:
(10)CI=DADAx+DBDBx,where when CI=1, the two drugs have an additive effect, whenCI<1, the two drugs have a synergistic effect, and when CI>1, the two drugs are antagonistic. DAandDB represent the concentrations of the two drugs, dAanddB, when administered in combination, and DAxandDBx represent the concentrations of the two drugs, dAanddB, when they have been administered alone and the combined drug inhibition rate has been reached.
(11)Dx=Dmfafu1/m.In Formula (11), faandfu are the inhibition and survival rates, respectively, when the two drugs are combined, fa=1−fu, Dm is the median dose, and m is the curve coefficient (calculated according to half-maximal inhibitory concentration (IC50) theory).In this study, CompuSyn software (http://www.combosyn.com/) was used to calculate the CI values of the drug combinations to evaluate whether each screened combination of Western and TCM drugs had synergistic therapeutic effects.


## 2.5.1. Cell Culture and Reagents

RAW264.7 cells were purchased from the American Type Culture Collection. The culture reagents included Dulbecco’s modified Eagle medium (DMEM; Gibco, C11995500CP), Roswell Park Memorial Institute (RPMI) 1640 (Gibco, C11875500BT), fetal bovine serum (Bio IND, 04-002-1A), antimycotic (Lifetechnologies, 15240-112), phosphate-buffered saline (PBS), pH 7.4 (Gibco, 10010-023), trypsin-EDTA (0.05%; Lifetechnologies, 25300-054), bovine serum albumin (Lifetechnologies, 15561012), and CellTiter-Glo (CTG; PROMEGA, G7572). Consumables and instruments included a cell culture plate (Corning), cell culture flask (Corning), microplate tester (BioTek, HM-1), and conventional instruments, such as a CO2 incubator (Thermo 3111) and biosafety cabinet (Heal Force, HFSAFe-1200LC).


## 2.5.2. Experimental Procedure and Methods

First, the different compounds were diluted with DMEM medium to the desired concentrations of the experimental design, namely, 0, 0.032, 0.16, 0.8, 4, 20, and 100μM. The cells were then counted and inoculated into 96-well culture plates at 2000 cells/100 μL/well; 90 μL of culture medium (including serum) was added to each well for overnight culture, and 10 μL of drug solution was added. After a 48 h culture, CTG was used to detect the proliferation of RAW264.7 cells; CTG (100 μL) was added to each well, the content of which was incubated at room temperature in a dark environment for 10 min. The chemiluminescence value of each well was measured at 500 ms using enzyme calibration.


## 2.5.3. Combination Index Measurement

In this study, according to the principle of the median effect—specifically by employing the Chou–Talalay model [1]—the synergistic effect between each set of drugs was analyzed, and the effects of integrated Western and TCM drugs on RA cell lines were evaluated using the combination index (CI). The CI indicates the degree of drug synergy, addition, or antagonism for any given drug combination:
(10)CI=DADAx+DBDBx,where when CI=1, the two drugs have an additive effect, whenCI<1, the two drugs have a synergistic effect, and when CI>1, the two drugs are antagonistic. DAandDB represent the concentrations of the two drugs, dAanddB, when administered in combination, and DAxandDBx represent the concentrations of the two drugs, dAanddB, when they have been administered alone and the combined drug inhibition rate has been reached.
(11)Dx=Dmfafu1/m.In Formula (11), faandfu are the inhibition and survival rates, respectively, when the two drugs are combined, fa=1−fu, Dm is the median dose, and m is the curve coefficient (calculated according to half-maximal inhibitory concentration (IC50) theory).In this study, CompuSyn software (http://www.combosyn.com/) was used to calculate the CI values of the drug combinations to evaluate whether each screened combination of Western and TCM drugs had synergistic therapeutic effects.


## 2.6. Statistical Analysis

Data were expressed as themean±standarddeviation, and a t-test and Mann–Whitney nonparametric tests were performed using IBM SPSS Statistics (version 26.0).


## 3. Results

### 3.1. Collection and Collation of Combined Drug Data

We obtained 946 pairs (comprising 816 drugs) of combinatorial drug data reported by the FDA from the DrugCombDB database. We further downloaded all existing small molecule drug data from the DrugBank database. After name matching and sorting, some small molecule drug data without key information, such as targets, were removed. In total, we included 488 drugs and their corresponding 546 pairs of combinations. The data obtained for each small molecule drug mainly included DrugBank ID, canonical SMILES, and known targets (Supplementary1; pretreated FDA-certified combination drug dataset).


### 3.2. Construction of Combination Drug Training Set

First, using the method described in Section2.2, we obtained the molecular fingerprint similarity value (SFP), the ADME index similarity value (SADME), and the minimum, maximum, median, and mean of the similarity (SSeqmin,SSeqmax,SSeqmed,andSSeqmean, respectively) of the drug target sequence between the two small molecule drugs. The minimum, maximum, median, and mean values of the GO functional similarity of the drug targets wereSGOmin,SGOmax,SGOmed,andSGOmean, respectively. The KEGG pathway similarity value (SPathway); the minimum, maximum, median, and mean of the shortest path of different drug targets on the human PPI network (SPPImin,SPPImax,SPPImed,andSPPImean, respectively); and the proximity (Snp) of the different drug targets were also obtained. Thus, we obtained a total of 16 characteristic variables to predict combinations of Western and TCM drugs.In this study, we evaluated 16 characteristic variables between two drugs for 546 pairs of combined drugs and removed those for which characteristic variables could not be calculated or were missing. We obtained a total of 404 positive samples of combined drugs (comprising 488 drugs) reported by the FDA. Due to a lack of effective noncombination drug data, we combined the 488 drugs in pairs, removed the existing combination drug samples and those missing data, and obtained a total of 45603 noncombination drugs as negative samples (Supplementary2; positive and negative sample training dataset for all drug combinations).Finally, to construct the training set for the TCM-combined drug samples, because the number of noncombined drugs should be much larger than that of combined drugs and to avoid large deviation in positive and negative samples, we set the ratio of positive to negative samples as 1 : 5; the number of positive samples was 404, and the number of negative samples randomly extracted from the 45603 noncombined drug samples was404×5=2020. Through this, a training sample set was formed containing 2424 samples, and the training set was repeated 1000 times to randomly generate 1000 prediction training sets of combined drugs, as displayed in Figure 2.Figure 2
Training set of prediction samples comprised 1000 training sets of prediction samples randomly selected from 2020 samples of noncombination drugs and 404 samples of combination drugs.


### 3.3. Analysis and Screening of Characteristic Variables in Combination Drug Training Set

At-test and a Mann–Whitney nonparametric test were used to determine whether significant differences existed between the 16 characteristic variables in the combined drug group and those in the noncombined drug group (p<0.05). The results revealed that SFP,SADME,SSeqmin,SGOmax,SGOmed,SGOmean,andSPPImax did not significantly differ between the combined drug and noncombined drug groups (p>0.05). Significant differences were found in the SSeqmax,SSeqmed,SSeqmean,SPathway,SPPImin,SPPImed,SPPImean,andSnp between the two statistical tests (p<0.05). Statistical differences were found in SGOmin in the nonparametric tests (p<0.05). All statistical tests were performed in IBM SPSS Statistics (version 26.0), and their data are presented in Table 1. The violin diagrams of the distribution of characteristic variables in different groups were completed using ggplot2 (version 3.3.5) and are displayed in Figure 3.Table 1
Results oft-test and Mann–Whitney nonparametric test for 16 characteristic variables between different groups.
FeatureNoncombination groupCombination groupt-testMann–Whitney testSFP0.783±0.0870.783±0.0910.9140.664SADME3.347±2.9943.160±2.6970.1670.199SSeqmin0.014±0.0450.020±0.0700.1090.400SSeqmax0.066±0.1890.173±0.342∗∗∗∗SSeqmed0.021±0.0570.030±0.084∗∗∗SSeqmean0.026±0.0620.049±0.105∗∗∗∗SGOmin0.236±0.1750.221±0.1930.114∗∗SGOmax0.523±0.2390.550±0.2880.0610.092SGOmed0.374±0.1720.369±0.2050.6490.532SGOmean0.376±0.1660.374±0.1940.8250.644SPathway0.007±0.0290.014±0.039∗∗∗∗SPPImin0.867±0.0920.821±0.142∗∗∗∗SPPImax0.931±0.1160.934±0.0940.5000.889SPPImed0.917±0.0430.910±0.050∗∗∗∗SPPImean0.100±0.0640.121±0.091∗∗∗∗Snp−0.487±0.215−0.436±0.216∗∗∗∗∗p<0.05and∗∗p<0.01.Figure 3
Violin diagram of value distribution of 16 characteristic variables in training sample set.Anrfe function using the caret package was used to screen the 16 characteristic variables. A random forest model was used as the selection function, which was functions=rfFuncs, and the method was cross-validation, which was method=“cv.” The others parameters were set to default. The training set was repeated 1000 times to extract feature variables from 1000 randomly generated training set samples, and each extracted feature variable set was recorded. The frequency of each feature variable in the 1000 feature selections was counted; the results are displayed in Figure 4.Figure 4
Occurrence of 16 feature variables from 1000 feature extractions.As illustrated in Figure4, eight characteristic variables, namely, SPathway,SGOmin,SGOmax,SSeqmin,SSeqmax,SSeqmean,SPPImean,andSnp, had appearance frequencies of over 70% after 1000 characteristic screenings. SPathwayandSGOmin were identified as key characteristic variables for predicting combined drugs. The other eight indicators had appearance frequencies below 70%; for example, SFP had a rate of less than 10%. Thus, we removed the group of indicators with frequencies lower than 70% and included the top eight indicators as feature variables for subsequent machine learning.


### 3.4. Comparison and Selection Results of Machine Learning Models

The five included machine learning models were subjected to three replicates of 10-fold cross-validation under 1000 random samples to evaluate the predictive classification performance of different models. The results are displayed in Table2, and the ROC and PR curves of the five models are presented in Figure 5.Table 2
Prediction evaluation results of five machine learning classification models under 1000 random sample training sets.
ModelAccuracyPrecisionRecallF1 scoreAUROCAUPRk-NN0.827±0.0120.848±0.0120.965±0.0130.903±0.0080.651±0.0240.309±0.033NB0.813±0.0140.845±0.0120.949±0.0130.894±0.0090.617±0.0260.269±0.031SVM0.834±0.0110.838±0.0120.992±0.0060.909±0.0070.657±0.0270.330±0.038RF0.849±0.0120.861±0.0120.977±0.0080.915±0.0070.739±0.0240.460±0.039AdaBoost0.833±0.0110.834±0.0120.998±0.0040.909±0.0070.667±0.0270.317±0.038Figure 5
ROC and PR curves of five machine learning models for classification prediction.
(a)(b)Through evaluation and comparison of the predictive ability of the different models, we found that the classification index scores of the random forest model were generally better than those of the other prediction classifiers. Therefore, the random forest model was adopted as the combined drug prediction classifier for this study. The parameter was represented asntree=300.


### 3.5. Prediction Results of Integrative Medicine

We entered “Rheumatoid Arthritis” into the TCMSP database as a search term and identified 41 small molecule types in TCM ingredients that have been used to treat RA. We collected the target treatment information for each small molecule and verified it using the Uniprot database (https://www.uniprot.org/). We collected the molecular formula, PubChem CID, CAS number, and canonical SMILES data of each small molecule from PubChem (Table 3).Table 3
Small molecules of 41 anti-RA TCM ingredients collected from TCMSP database.
TCMSP IDNameMolecular formulaPubChem CIDCASMOL006505(-)-EpicatechinC15H14O672276490-46-0MOL006821(-)-Epigallocatechin gallateC22H18O1165064989-51-5MOL008680AcetaldehydeC2H4O17775-07-0MOL000475AnetholeC10H12O637563104-46-1MOL000008ApigeninC15H10O55280443520-36-5MOL002773Beta-caroteneC40H5652804897235-40-7MOL000358Beta-SitosterolC29H50O22228483-46-5MOL003973CaffeineC8H10N4O2251958-08-2MOL008842Chenodeoxycholic acidC24H40O410133474-25-9MOL000390DaidzeinC15H10O45281708486-66-8MOL000254EugenolC10H12O2331497-53-0MOL013179FisetinC15H10O65281614528-48-3MOL000392FormononetinC16H12O45280378485-72-3MOL000481GenisteinC15H10O55280961446-72-0MOL002467GingerolC17H26O444279323513-14-6MOL000666HexanalC6H12O618466-25-1MOL005916IrisolidoneC17H14O652817812345-17-7MOL000422KaempferolC15H10O65280863520-18-3MOL000305Lauric acidC12H24O23893143-07-7MOL000006LuteolinC15H10O65280445491-70-3MOL007990MilitarinC34H46O1717163858139-23-4MOL002008MyricetinC15H10O85281672529-44-2MOL003493NaphthaleneC10H893191-20-3MOL003403NicotineC10H14N28959454-11-5MOL005828NobiletinC21H22O872344478-01-3MOL006214ProgesteroneC21H30O2599457-83-0MOL012297PuerarinC21H20O952818073681-99-0MOL000098QuercetinC15H10O75280343117-39-5MOL012744ResveratrolC14H12O3445154501-36-0MOL002268RheinC15H8O610168478-43-3MOL011865Rosmarinic acidC18H16O8528179220283-92-5MOL006356SorbitolC6H14O6578050-70-4MOL007154Tanshinone IIAC19H18O3164676568-72-9MOL003186TripterineC29H38O412272434157-83-0MOL003187TriptolideC20H24O610798538748-32-2MOL004932Uralsaponin AC42H62O16128229103000-77-7MOL000511Ursolic acidC30H48O36494577-52-1MOL000635VanillinC8H8O31183121-33-5MOL000173WogoninC16H12O55281703632-85-9MOL009357Yakuchinone AC20H24O313314578954-23-1MOL009358Yakuchinone BC20H22O3644036581840-57-5The information for drugs commonly used in RA treatment was also obtained from DrugBank. In this study, acetaminophen, celecoxib, choline salicylate, diclofenac, hydroxychloroquine, indometacin, leflunomide, methotrexate, naproxen, sulfasalazine, and 10 other common drugs were selected, as presented in Table4.Table 4
Information from DrugBank database for 10 common anti-RA small molecule drugs.
IDNameMolecular formulaPubChem CIDCASDB00316AcetaminophenC8H9NO21983103-90-2DB00328IndometacinC19H16ClNO4371553-86-1DB00482CelecoxibC17H14F3N3O2S2662169590-42-5DB00563MethotrexateC20H22N8O512694159-05-2DB00586DiclofenacC14H11Cl2NO2303315307-86-5DB00788NaproxenC14H14O315639122204-53-1DB00795SulfasalazineC18H14N4O5S5339599-79-1DB01097LeflunomideC12H9F3N2O2389975706-12-6DB01611HydroxychloroquineC18H26ClN3O3652118-42-3DB14006Choline salicylateC12H19NO4546863502016-36-6The aforementioned 41 small molecule TCM (collected from TCMSP) and 10 common anti-RA Western medicine drugs were combined in pairs to obtain potential drug combinations:C411×C101=410 combinations. Thus, we obtained the final prediction sample set (Supplementary 3; drug combination prediction dataset of integrated Western and TCM drugs).The results revealed that out of 1000 random sample predictions, 34, 62, and 114 pairs of Western and TCM drugs were predicted to be combination drugs more than 500, 100–500, and 1–100 times, and 200 groups were predicted to be noncombination drugs. In this study, we considered predictions of more than 100 times to be potential anti-RA Western and TCM combination drugs, with higher numbers of positive sample predictions indicating a greater potential for drug combination. Therefore, we selected the 15 groups that were predicted more than 700 times for subsequent cell experiments. The combination drug information is displayed in Table5 (Supplementary 4; prediction results of all Western and TCM combinations).Table 5
Combination data for 15 most predicted drug combinations.
GroupTCM ingredientAnti-RA drugCount1EugenolAcetaminophen9932(-)-EpicatechinCelecoxib9573Rosmarinic acidCelecoxib9464Chenodeoxycholic acidAcetaminophen9295EugenolDiclofenac9126EugenolNaproxen9107MyricetinCelecoxib9088Chenodeoxycholic acidDiclofenac8989Chenodeoxycholic acidNaproxen89710RheinCelecoxib89411AnetholeNaproxen85412AnetholeDiclofenac85313NobiletinCelecoxib83214FisetinCelecoxib79115RheinHydroxychloroquine770


### 3.6. Experimental Verification Results

Because RAW264.7 cells were used for experimental verification, we first determined whether a single small molecule had an inhibitory effect on RAW264.7 cells. We then tested the Western and TCM drug combinations. Finally, the CI value was calculated using CompuSyn software to determine whether the small molecule combinations of the two drugs were synergetic.We then applied the method described in Section2.5. Based on observations of the proliferation effect of the 15 small molecules on the RAW264.7 cells, myricetin, nobiletin, fisetin, celecoxib, and hydroxychloroquine significantly reduced the survival rate of the RAW264.7 cells. The IC50 values of the myricetin, nobiletin, fisetin, celecoxib, and hydroxychloroquine were> 100, 23.93, 31.19, >100, and 79.43 μmol/L, respectively, as presented in Figure 6.Figure 6
Effects of myricetin, nobiletin, fisetin, celecoxib, hydroxychloroquine, and rhein on the viability of RAW264.7 cells.
(a)(b)(c)(d)(e)(f)We further tested combinations of Western and TCM containing these five small molecules (a total of seven pairs) based on the prediction results of the machine learning model and calculated the CI values, as displayed in Table6. We found that of the seven pairs of Western and TCM combinations, five of the pairs had CI values less than 0.7, indicating that these drugs have good synergy.Table 6
Effect of seven predicted pairs of Western and TCM combinations on RAW264.7 cell viability.
CombinationTCM (μmol/L)Drug (μmol/L)IC50RAW264.7 cell inhibition rate (%)CIRhein+celecoxib>100+100CI<110010039.31±2.92CI=0.372010035.68±1.49CI=0.49410038.60±0.80CI=0.391002020.08±1.22CI=0.40100414.22±1.59CI=0.17202020.17±1.49CI=0.40Nobiletin+celecoxib8.83+100CI<110010096.20±0.24CI=0.011002061.50±1.67CI=0.17Myricetin+celecoxib100+4CI<110010064.05±0.94CI=0.062010037.13±2.13CI=0.46410037.96±1.58CI=0.411002053.92±1.59CI=0.03100449.97±0.93CI=0.02Fisetin+celecoxib10.89+100CI<110010099.06±0.03CI=0.02410057.80±1.89CI=0.88Rhein+hydroxychloroquine100+57.15CI<11002086.88±0.04CI=0.28Epicatechin+celecoxib>100+100CI>1Rosmarinic acid+celecoxib>100+100CI>1


## 3.1. Collection and Collation of Combined Drug Data

We obtained 946 pairs (comprising 816 drugs) of combinatorial drug data reported by the FDA from the DrugCombDB database. We further downloaded all existing small molecule drug data from the DrugBank database. After name matching and sorting, some small molecule drug data without key information, such as targets, were removed. In total, we included 488 drugs and their corresponding 546 pairs of combinations. The data obtained for each small molecule drug mainly included DrugBank ID, canonical SMILES, and known targets (Supplementary1; pretreated FDA-certified combination drug dataset).


## 3.2. Construction of Combination Drug Training Set

First, using the method described in Section2.2, we obtained the molecular fingerprint similarity value (SFP), the ADME index similarity value (SADME), and the minimum, maximum, median, and mean of the similarity (SSeqmin,SSeqmax,SSeqmed,andSSeqmean, respectively) of the drug target sequence between the two small molecule drugs. The minimum, maximum, median, and mean values of the GO functional similarity of the drug targets wereSGOmin,SGOmax,SGOmed,andSGOmean, respectively. The KEGG pathway similarity value (SPathway); the minimum, maximum, median, and mean of the shortest path of different drug targets on the human PPI network (SPPImin,SPPImax,SPPImed,andSPPImean, respectively); and the proximity (Snp) of the different drug targets were also obtained. Thus, we obtained a total of 16 characteristic variables to predict combinations of Western and TCM drugs.In this study, we evaluated 16 characteristic variables between two drugs for 546 pairs of combined drugs and removed those for which characteristic variables could not be calculated or were missing. We obtained a total of 404 positive samples of combined drugs (comprising 488 drugs) reported by the FDA. Due to a lack of effective noncombination drug data, we combined the 488 drugs in pairs, removed the existing combination drug samples and those missing data, and obtained a total of 45603 noncombination drugs as negative samples (Supplementary2; positive and negative sample training dataset for all drug combinations).Finally, to construct the training set for the TCM-combined drug samples, because the number of noncombined drugs should be much larger than that of combined drugs and to avoid large deviation in positive and negative samples, we set the ratio of positive to negative samples as 1 : 5; the number of positive samples was 404, and the number of negative samples randomly extracted from the 45603 noncombined drug samples was404×5=2020. Through this, a training sample set was formed containing 2424 samples, and the training set was repeated 1000 times to randomly generate 1000 prediction training sets of combined drugs, as displayed in Figure 2.Figure 2
Training set of prediction samples comprised 1000 training sets of prediction samples randomly selected from 2020 samples of noncombination drugs and 404 samples of combination drugs.


## 3.3. Analysis and Screening of Characteristic Variables in Combination Drug Training Set

At-test and a Mann–Whitney nonparametric test were used to determine whether significant differences existed between the 16 characteristic variables in the combined drug group and those in the noncombined drug group (p<0.05). The results revealed that SFP,SADME,SSeqmin,SGOmax,SGOmed,SGOmean,andSPPImax did not significantly differ between the combined drug and noncombined drug groups (p>0.05). Significant differences were found in the SSeqmax,SSeqmed,SSeqmean,SPathway,SPPImin,SPPImed,SPPImean,andSnp between the two statistical tests (p<0.05). Statistical differences were found in SGOmin in the nonparametric tests (p<0.05). All statistical tests were performed in IBM SPSS Statistics (version 26.0), and their data are presented in Table 1. The violin diagrams of the distribution of characteristic variables in different groups were completed using ggplot2 (version 3.3.5) and are displayed in Figure 3.Table 1
Results oft-test and Mann–Whitney nonparametric test for 16 characteristic variables between different groups.
FeatureNoncombination groupCombination groupt-testMann–Whitney testSFP0.783±0.0870.783±0.0910.9140.664SADME3.347±2.9943.160±2.6970.1670.199SSeqmin0.014±0.0450.020±0.0700.1090.400SSeqmax0.066±0.1890.173±0.342∗∗∗∗SSeqmed0.021±0.0570.030±0.084∗∗∗SSeqmean0.026±0.0620.049±0.105∗∗∗∗SGOmin0.236±0.1750.221±0.1930.114∗∗SGOmax0.523±0.2390.550±0.2880.0610.092SGOmed0.374±0.1720.369±0.2050.6490.532SGOmean0.376±0.1660.374±0.1940.8250.644SPathway0.007±0.0290.014±0.039∗∗∗∗SPPImin0.867±0.0920.821±0.142∗∗∗∗SPPImax0.931±0.1160.934±0.0940.5000.889SPPImed0.917±0.0430.910±0.050∗∗∗∗SPPImean0.100±0.0640.121±0.091∗∗∗∗Snp−0.487±0.215−0.436±0.216∗∗∗∗∗p<0.05and∗∗p<0.01.Figure 3
Violin diagram of value distribution of 16 characteristic variables in training sample set.Anrfe function using the caret package was used to screen the 16 characteristic variables. A random forest model was used as the selection function, which was functions=rfFuncs, and the method was cross-validation, which was method=“cv.” The others parameters were set to default. The training set was repeated 1000 times to extract feature variables from 1000 randomly generated training set samples, and each extracted feature variable set was recorded. The frequency of each feature variable in the 1000 feature selections was counted; the results are displayed in Figure 4.Figure 4
Occurrence of 16 feature variables from 1000 feature extractions.As illustrated in Figure4, eight characteristic variables, namely, SPathway,SGOmin,SGOmax,SSeqmin,SSeqmax,SSeqmean,SPPImean,andSnp, had appearance frequencies of over 70% after 1000 characteristic screenings. SPathwayandSGOmin were identified as key characteristic variables for predicting combined drugs. The other eight indicators had appearance frequencies below 70%; for example, SFP had a rate of less than 10%. Thus, we removed the group of indicators with frequencies lower than 70% and included the top eight indicators as feature variables for subsequent machine learning.


## 3.4. Comparison and Selection Results of Machine Learning Models

The five included machine learning models were subjected to three replicates of 10-fold cross-validation under 1000 random samples to evaluate the predictive classification performance of different models. The results are displayed in Table2, and the ROC and PR curves of the five models are presented in Figure 5.Table 2
Prediction evaluation results of five machine learning classification models under 1000 random sample training sets.
ModelAccuracyPrecisionRecallF1 scoreAUROCAUPRk-NN0.827±0.0120.848±0.0120.965±0.0130.903±0.0080.651±0.0240.309±0.033NB0.813±0.0140.845±0.0120.949±0.0130.894±0.0090.617±0.0260.269±0.031SVM0.834±0.0110.838±0.0120.992±0.0060.909±0.0070.657±0.0270.330±0.038RF0.849±0.0120.861±0.0120.977±0.0080.915±0.0070.739±0.0240.460±0.039AdaBoost0.833±0.0110.834±0.0120.998±0.0040.909±0.0070.667±0.0270.317±0.038Figure 5
ROC and PR curves of five machine learning models for classification prediction.
(a)(b)Through evaluation and comparison of the predictive ability of the different models, we found that the classification index scores of the random forest model were generally better than those of the other prediction classifiers. Therefore, the random forest model was adopted as the combined drug prediction classifier for this study. The parameter was represented asntree=300.


## 3.5. Prediction Results of Integrative Medicine

We entered “Rheumatoid Arthritis” into the TCMSP database as a search term and identified 41 small molecule types in TCM ingredients that have been used to treat RA. We collected the target treatment information for each small molecule and verified it using the Uniprot database (https://www.uniprot.org/). We collected the molecular formula, PubChem CID, CAS number, and canonical SMILES data of each small molecule from PubChem (Table 3).Table 3
Small molecules of 41 anti-RA TCM ingredients collected from TCMSP database.
TCMSP IDNameMolecular formulaPubChem CIDCASMOL006505(-)-EpicatechinC15H14O672276490-46-0MOL006821(-)-Epigallocatechin gallateC22H18O1165064989-51-5MOL008680AcetaldehydeC2H4O17775-07-0MOL000475AnetholeC10H12O637563104-46-1MOL000008ApigeninC15H10O55280443520-36-5MOL002773Beta-caroteneC40H5652804897235-40-7MOL000358Beta-SitosterolC29H50O22228483-46-5MOL003973CaffeineC8H10N4O2251958-08-2MOL008842Chenodeoxycholic acidC24H40O410133474-25-9MOL000390DaidzeinC15H10O45281708486-66-8MOL000254EugenolC10H12O2331497-53-0MOL013179FisetinC15H10O65281614528-48-3MOL000392FormononetinC16H12O45280378485-72-3MOL000481GenisteinC15H10O55280961446-72-0MOL002467GingerolC17H26O444279323513-14-6MOL000666HexanalC6H12O618466-25-1MOL005916IrisolidoneC17H14O652817812345-17-7MOL000422KaempferolC15H10O65280863520-18-3MOL000305Lauric acidC12H24O23893143-07-7MOL000006LuteolinC15H10O65280445491-70-3MOL007990MilitarinC34H46O1717163858139-23-4MOL002008MyricetinC15H10O85281672529-44-2MOL003493NaphthaleneC10H893191-20-3MOL003403NicotineC10H14N28959454-11-5MOL005828NobiletinC21H22O872344478-01-3MOL006214ProgesteroneC21H30O2599457-83-0MOL012297PuerarinC21H20O952818073681-99-0MOL000098QuercetinC15H10O75280343117-39-5MOL012744ResveratrolC14H12O3445154501-36-0MOL002268RheinC15H8O610168478-43-3MOL011865Rosmarinic acidC18H16O8528179220283-92-5MOL006356SorbitolC6H14O6578050-70-4MOL007154Tanshinone IIAC19H18O3164676568-72-9MOL003186TripterineC29H38O412272434157-83-0MOL003187TriptolideC20H24O610798538748-32-2MOL004932Uralsaponin AC42H62O16128229103000-77-7MOL000511Ursolic acidC30H48O36494577-52-1MOL000635VanillinC8H8O31183121-33-5MOL000173WogoninC16H12O55281703632-85-9MOL009357Yakuchinone AC20H24O313314578954-23-1MOL009358Yakuchinone BC20H22O3644036581840-57-5The information for drugs commonly used in RA treatment was also obtained from DrugBank. In this study, acetaminophen, celecoxib, choline salicylate, diclofenac, hydroxychloroquine, indometacin, leflunomide, methotrexate, naproxen, sulfasalazine, and 10 other common drugs were selected, as presented in Table4.Table 4
Information from DrugBank database for 10 common anti-RA small molecule drugs.
IDNameMolecular formulaPubChem CIDCASDB00316AcetaminophenC8H9NO21983103-90-2DB00328IndometacinC19H16ClNO4371553-86-1DB00482CelecoxibC17H14F3N3O2S2662169590-42-5DB00563MethotrexateC20H22N8O512694159-05-2DB00586DiclofenacC14H11Cl2NO2303315307-86-5DB00788NaproxenC14H14O315639122204-53-1DB00795SulfasalazineC18H14N4O5S5339599-79-1DB01097LeflunomideC12H9F3N2O2389975706-12-6DB01611HydroxychloroquineC18H26ClN3O3652118-42-3DB14006Choline salicylateC12H19NO4546863502016-36-6The aforementioned 41 small molecule TCM (collected from TCMSP) and 10 common anti-RA Western medicine drugs were combined in pairs to obtain potential drug combinations:C411×C101=410 combinations. Thus, we obtained the final prediction sample set (Supplementary 3; drug combination prediction dataset of integrated Western and TCM drugs).The results revealed that out of 1000 random sample predictions, 34, 62, and 114 pairs of Western and TCM drugs were predicted to be combination drugs more than 500, 100–500, and 1–100 times, and 200 groups were predicted to be noncombination drugs. In this study, we considered predictions of more than 100 times to be potential anti-RA Western and TCM combination drugs, with higher numbers of positive sample predictions indicating a greater potential for drug combination. Therefore, we selected the 15 groups that were predicted more than 700 times for subsequent cell experiments. The combination drug information is displayed in Table5 (Supplementary 4; prediction results of all Western and TCM combinations).Table 5
Combination data for 15 most predicted drug combinations.
GroupTCM ingredientAnti-RA drugCount1EugenolAcetaminophen9932(-)-EpicatechinCelecoxib9573Rosmarinic acidCelecoxib9464Chenodeoxycholic acidAcetaminophen9295EugenolDiclofenac9126EugenolNaproxen9107MyricetinCelecoxib9088Chenodeoxycholic acidDiclofenac8989Chenodeoxycholic acidNaproxen89710RheinCelecoxib89411AnetholeNaproxen85412AnetholeDiclofenac85313NobiletinCelecoxib83214FisetinCelecoxib79115RheinHydroxychloroquine770


## 3.6. Experimental Verification Results

Because RAW264.7 cells were used for experimental verification, we first determined whether a single small molecule had an inhibitory effect on RAW264.7 cells. We then tested the Western and TCM drug combinations. Finally, the CI value was calculated using CompuSyn software to determine whether the small molecule combinations of the two drugs were synergetic.We then applied the method described in Section2.5. Based on observations of the proliferation effect of the 15 small molecules on the RAW264.7 cells, myricetin, nobiletin, fisetin, celecoxib, and hydroxychloroquine significantly reduced the survival rate of the RAW264.7 cells. The IC50 values of the myricetin, nobiletin, fisetin, celecoxib, and hydroxychloroquine were> 100, 23.93, 31.19, >100, and 79.43 μmol/L, respectively, as presented in Figure 6.Figure 6
Effects of myricetin, nobiletin, fisetin, celecoxib, hydroxychloroquine, and rhein on the viability of RAW264.7 cells.
(a)(b)(c)(d)(e)(f)We further tested combinations of Western and TCM containing these five small molecules (a total of seven pairs) based on the prediction results of the machine learning model and calculated the CI values, as displayed in Table6. We found that of the seven pairs of Western and TCM combinations, five of the pairs had CI values less than 0.7, indicating that these drugs have good synergy.Table 6
Effect of seven predicted pairs of Western and TCM combinations on RAW264.7 cell viability.
CombinationTCM (μmol/L)Drug (μmol/L)IC50RAW264.7 cell inhibition rate (%)CIRhein+celecoxib>100+100CI<110010039.31±2.92CI=0.372010035.68±1.49CI=0.49410038.60±0.80CI=0.391002020.08±1.22CI=0.40100414.22±1.59CI=0.17202020.17±1.49CI=0.40Nobiletin+celecoxib8.83+100CI<110010096.20±0.24CI=0.011002061.50±1.67CI=0.17Myricetin+celecoxib100+4CI<110010064.05±0.94CI=0.062010037.13±2.13CI=0.46410037.96±1.58CI=0.411002053.92±1.59CI=0.03100449.97±0.93CI=0.02Fisetin+celecoxib10.89+100CI<110010099.06±0.03CI=0.02410057.80±1.89CI=0.88Rhein+hydroxychloroquine100+57.15CI<11002086.88±0.04CI=0.28Epicatechin+celecoxib>100+100CI>1Rosmarinic acid+celecoxib>100+100CI>1


## 4. Discussion

Machine learning methods have been widely used in many fields of research, such as text classification [23], sentiment analysis [24], breast cancer diagnosis [25], web page classification [26], and image generation [27]. Many researchers have used computational and modeling methods to identify potential combinatorial drugs [17, 28–30], and the use of machine learning to predict drug combinations has notably increased in recent years [31–34]. However, nearly all of this research has focused of combinations of Western drugs. Combinations of Western and TCM drugs have rarely been reported, especially within the context of RA treatment. TCM is a complex system composed of many small molecules with synergistic effects. In recent years, the research of TCM small molecule is still the focus of modern pharmacology. We mainly focused on the prediction of the combination of small molecule with clear pharmacological effects in TCM and Western medicines that have been approved drugs. Therefore, in this study, we examined the small molecule characteristics of TCM ingredients that differed from those in Western combination drugs. We developed a method for complete modeling and prediction to screen anti-RA combination drugs.For feature variable construction, we selected drug structures and ADME features. However, the results of the statistical analyses and feature variable screening revealed that drug structure similarity and ADME characteristics had no significant effect on the classification of combination drugs. Notably, sequence similarity, GO functional similarity, and the distance of drug targets within the PPI network of the small molecule drug targets were found to be related to the combination of drugs, with the similarity value of the drug target pathway (SPathway) being identified as a key feature in combination drugs, both statistically and in feature screening. Therefore, the synergistic mechanism of two or more drugs within a pathway may serve as a key reference for future developments of combined drugs.In constructing the sample training set, we were unable to collect sufficient noncombined drug data (namely the number of negative samples) due to limited research data. Therefore, in the original data, the number of positive samples was only 1% that of negative samples, which created a serious imbalance in the number of positive and negative samples. If the model had been constructed using the original datasets, the disproportionate number of negative samples would have caused overfitting, and the predictions would have been biased toward the classifications of the larger number of samples. This would have greatly reduced the normalizing ability of the models. However, simply reducing the sample size may have excluded potential drug combinations in the 45603 samples.To effectively and reasonably overcome this challenge, we proposed a method in which multiple training sets were constructed using random sampling performed 1000 times. For each instance of modeling, the positive samples remained unchanged, and a limited number of negative samples were randomly selected from the total negative samples; 1000 training sets were randomly composed. The machine learning model was trained by and predicted the 1000 training sets. The final potential combination drugs were then determined on the basis of the percentage of each sample predicted to be a combination drug (positive) in the 1000 learning sessions. However, because the number of noncombination drugs is generally much larger than that of combination drugs but large deviations in numbers of positive and negative samples would cause prediction deviation in the machine learning model, in this study, the ratio of positive and negative samples was set to 1 : 5. Through these changes, we were able to prevent imperfect data sampling of combined drugs.In this study, from the prediction and evaluation results of five different machine learning classification models in 1000 random sample training sets, the accuracy, precision,F1 score, AUROC, and AUPR values of RF are the highest among the five machine learning models. AdaBoost performs best on recall, followed by SVM and RF. The prediction scores of k-NN and NB on these six machine learning evaluation indexes were lower, suggesting that these two models may not be suitable for combination drug prediction. The prediction performance of SVM and AdaBoost is second only to that of RF. Therefore, through the evaluation and comparison of the predictive ability of different models, we found that the RF classification index scores are generally better than other predictive classifiers. Therefore, the RF model is used as the final combined drug prediction classifier, where the parameter ntree=300. We used the RF classifier to predict 410 groups of Western and TCM drug combinations to treat RA and thus preliminarily predicted 96 pairs of potential Western and TCM drug combinations to treat RA. We used RAW264.7 cells to experimentally verify the effects of a single drug small molecule on the proliferation of RAW264.7 cells. We then verified the effects of combined drugs on the cells. The results revealed that five of the seven selected pairs of combined Western and TCM drugs had high synergistic effects, indicating that our model had high accuracy and practicability. The five optimal combinations were Rhein and Celecoxib, Nobiletin and Celecoxib, Myricetin and Celecoxib, Fisetin and Celecoxib, and Rhein and Hydroxychloroquine.Notably, celecoxib was found to have a higher synergistic therapeutic effect when combined with several small molecule TCM ingredients, namely, myricetin, rhein, nobiletin, and fisetin. In addition, rhein did not demonstrate an inhibitory effect on RAW264.7 cell proliferation when applied alone. However, when combined with celecoxib (CI<0.7), the combined inhibition rate of 100 μmol/L rhein and 20 μmol/L celecoxib was 35.68%, which was higher than the rate of each of the two alone. Rhien had high synergy with hydroxychloroquine.Notably, triptolide combined with Western drugs, such as methotrexate, has previously been reported to produce favorable results in past experiments [35]. However, the number of predicted combinations containing triptolide in our model was not high; it predicted triptolide combined with diclofenac (75 times), naproxen (75 times), and celecoxib (24 times) but did not predict triptolide combined with methotrexate (0 times). This may be due to triptolide’s reported clinical toxicity [36]. Although triptolide has been clinically reported to improve treatment outcomes, it may not be the best option for drug combinations and may even lead to more systemic side effects [37]. Our prediction results may support this hypothesis. The screening of potential combination of TCM and Western medicine based on the machine learning combined with experimental validation methods can improve efficiency and reduce research costs.In conclusion, we constructed combination drug characteristics based on existing combination drug data and characteristics from TCM drug research. Key characteristic indexes were screened, and a random forest classification model was used to predict potential combinations of Western and TCM drugs for anti-RA treatment. Through experimental verification, we preliminarily demonstrated that this method of combination therapy for anti-RA treatment can be applied to clinical practice. In the future, we will conduct vivo experiments and clinical trials to validate the findings of this study. In addition, our machine learning model for integrated Western and TCM drug combination can be applied for predicting not only anti-RA drug combinations but also TCM combinations for other diseases.



---

*Source: 6086388-2023-02-15.xml*