Phar-LSTM: a pharmacological representation-based LSTM network for drug–drug interaction extraction

Datasets

Text mining and natural language processing have recently benefitted the pharmacological industry. The DDIExtraction 2011 (Segura-Bedmar, Martínez Fernández & Sánchez Cisneros, 2011) and DDIExtraction 2013 (Segura-Bedmar, Martínez Fernández & Herrero Zazo, 2013) challenge tasks are held to promote the research of DDI extraction by providing benchmark datasets and enabling researchers to compare their methods fairly. The DDIExtraction 2011 challenge focuses on the binary classification of DDIs, that is, deciding whether the given document contains DDIs. For the DDIExtraction 2013 challenge, researchers must identify five DDI types: ADVICE, EFFECT, MECHANISM, INT, and NEGATIVE, which can be considered to be a multi-class classification problem.

The DDIExtraction 2011 dataset includes 579 documents about 14,949 drugs from DrugBank. These DrugBank documents contain rich chemical and pharmaceutical information. There are 5,806 sentences containing 3,160 DDIs (binary) in the DDIExtraction 2011 dataset. The DDIExtraction 2013 dataset has 784 documents from DrugBank and 233 abstract documents from MedLine, with a total of 5,021 DDIs (five specific DDI types). We selected task 9.2 as the testing dataset. More details of the datasets, including the training and testing information, can be found in Segura-Bedmar, Martínez Fernández & Sánchez Cisneros (2011) and Segura-Bedmar, Martínez Fernández & Herrero Zazo (2013).

Embedding based pharmacological representation

To obtain useful pharmacological information from the biomedical literature, we present a newly defined embedding strategy to convert the raw input (biomedical documents) into high-dimensional vectors, that is, the pharmacological representations. It should be noted that a major difference between DDI extraction and other natural language processing tasks is that the two target drug entities in the DDI instance should be fully considered since the target drug pair contains important pharmacological information. For a document containing n drugs, there are $C_n^2$ DDI candidates. A document may contain more than one DDI instance, and all DDI candidates in the same document are expected to differ from each other. A common way to represent a DDI is “drug blinding”, that is, replacing the two target drugs with “drug1” and “drug2”, and the other drugs in the document are represented as “durg0” (Liu et al., 2016). However, this drug blinding strategy may discard some valuable pharmacological information contained in the target drugs (e.g., the distinguishing information between different target drugs).

The pharmacological representation framework is shown in Fig. 2. First, we tokenized the documents into token sequences. It should be noted that different from the drug blinding strategy, the target drugs are not replaced with the words “drug1” and “drug2”. Therefore, more pharmacological information is extracted. The t th token unit was mapped to the $x_t^{\left(\mathrm{token}\right)}$ using the word embedding strategy. We transformed each token into a d dimensional vector through random encoding as inspired by Wieting & Kiela (2019), wherein the values for each dimension are in $\left[-1/\sqrt{d},1/\sqrt{d}\right]$. In practice, we set d = 400 in our experiments.

Second, to extract the target drug related information, a set of filters are introduced to obtain the corresponding information in four aspects from the token sentence. The target drug filter determines whether a token is a target drug or not (1: True; 0: False), and the neighbor filter determines whether a token is a neighbor of the target drug (1: True; 0: False). Whether a token exists between a pair of target drugs (1: True; 0: False) is determined by the path filter. The importance filter measures the degree of the token associated with the nearest target drug. The closer the distance, the higher the degree (i.e., more importance is attached).

We define the metric I = 1/(r + 1)², where r is the distance from the token to the nearest target drug (e.g., if the token is the target drug, r = 0; if the token is located near to the target drug, r = 1). Similar to the vector space model (Salton, Wong & Yang, 1975), we characterize the target drug related information of the t th token as a four-dimensional vector $x_t^{\left(\mathrm{Phar}\right)}$.

Finally, a given document $D=\left(w_0,w_1,\dots ,w_m\right)$ with m words is represented by $D^{\left(\mathrm{input}\right)}=\left(\left(x_0^{\left(\mathrm{token}\right)},x_0^{\left(\mathrm{Phar}\right)}\right),\left(x_1^{\left(\mathrm{token}\right)},x_1^{\left(\mathrm{Phar}\right)}\right),\dots ,\left(x_n^{\left(\mathrm{token}\right)},x_n^{\left(\mathrm{Phar}\right)}\right)\right)$, which we call pharmacological representation. Usually, n is larger than m because a document may be tokenized by splitting using punctuation, such as “-” and “.”.

LSTM-based multi-task learning

Most previous studies tackled two tasks (i.e., the coarse-grained task and fine-grained task) separately. Because multi-task learning may learn the common features of the related tasks that would benefit each task (Caruana, 1997), here, we present an LSTM-based multi-task learning scheme as shown in Fig. 3 (flowing from the bottom up). The input layer converts the raw input into the pharmacological representations. For the recurrent layer, we adopt a special Recurrent Neural Network (RNN) structure (LSTM) (Van Houdt, Mosquera & Nápoles, 2020), which can store the previous information for a long time in data processing. Note that each recurrent block of the recurrent layer can be assigned a different LSTM variant, which is illustrated in Fig. 4. The multi-task output layer feeds the common features extracted from the hidden layer into the corresponding SoftMax classifiers. We pretrain the parameters of the neural network(except the multi-task output layer) and then fine-tune in the classification stage.

The Phar-LSTM block contains four gates and a cell state, whose variation that we use here (Fig. 4A) is formulated as (1)${\displaystyle z_t=\left(H_{t-1},x_t^{\left(\mathrm{token}\right)}\right),}$ (2)${\displaystyle P_t=\left(g\left(W^P{x}^{\left(\mathrm{Phar}\right)}+b^P\right),z_t\right),}$ (3)${\displaystyle C_t=g\left(W^I{P}_t+b^I\right)\circ C_{t-1}+g\left(W^F{P}_t+b^F\right)\circ h\left(W^C{p}_t+b^C\right),}$ (4)${\displaystyle O_t=g\left(W^O{P}_t+b^O\right),}$ (5)${\displaystyle H_t=O_t\circ h\left(C_t\right).}$

In the above, $x_t^{\left(\mathrm{token}\right)}$ and $x_t^{\left(\mathrm{Phar}\right)}$ represent the corresponding vectors (i.e., the pharmacological representation) generated by the t th token using the embedding strategy (see Embedding based pharmacological representation). W^P, W^C, W^I, W^F, and W^O are the weight matrices for the pharmacological gate, the cell state, the input gate, the forget gate, and the output gate separately, and b^P, b^C, b^I, b^F, and b^O, respectively, are the corresponding bias units. The functions g and h are activation functions. The sigmoid function is usually used as g for the four gates, and the tangent function is typically used as h for the cell state. The ∘ denotes point-wise multiplication.

Other LSTM variants

Different structures of LSTM may influence the results; therefore, we studied the DDI extraction performance of the models with different LSTM variants. The LSTM variants can be derived by modifying the gates, activation function, and connections. The derived 10 variants are shown in Fig. 4, and the details are as follows:

(1) No pharmacological gate (NPG)

Graves & Schmidhuber (2005) originally proposed the LSTM, which is also known as vanilla LSTM. Phar-LSTM can be transformed to the vanilla LSTM by removing the pharmacological gate: (6)${\displaystyle C_t=g\left(W^I{z}_t+b^I\right)\circ C_{t-1}+g\left(W^F{z}_t+b^F\right)\circ h\left(W^C{z}_t+b^C\right).}$

(2) No input gate (NIG)

By removing the input gate, we obtain a lighter version of C_t.C_t conveys less information to the next node: (7)${\displaystyle C_t=C_{t-1}+g\left(W^F{P}_t+b^F\right)\circ h\left(W^C{P}_t+b^C\right).}$

(3) No Input Activation Function (NIAF)

By removing the activation function of I_t, we obtain a “wilder” version of I_t since I_t is no longer confined to $\left[-1,1\right]$ by the sigmoid function: (8)${\displaystyle C_t=\left(W^I{P}_t+b^I\right)\circ C_{t-1}+g\left(W^F{P}_t+b^F\right)\circ h\left(W^C{P}_t+b^C\right).}$

(4) No forget gate (NFG)

Gers, Schmidhuber & Cummins (2000) first proposed a variant of LSTM by adding a forget gate, which enabled the LSTM to better forget the history information. By removing the forget gate, we obtain a lighter version of C_t. The C_t of NFG can remember more information because the function of the forget gate is to restrain the useless information from persisting in the history. (9)${\displaystyle C_t=g\left(W^I{P}_t+b^I\right)\circ C_{t-1}+h\left(W^C{P}_t+b^C\right).}$

(5) No forget activation function (NFAF)

Similar to the NIAF, we obtain a “wilder” version of F_t by removing the activation function of F_t: (10)${\displaystyle C_t=g\left(W^I{P}_t+b^I\right)\circ C_{t-1}+\left(W^F{P}_t+b^F\right)\circ h\left(W^C{P}_t+b^C\right).}$

(6) No output gate (NOG)

Similar to the NIG, by removing the output gate, we obtain a lighter version of H_t, and C_t conveys less information to the next node: (11)${\displaystyle H_t=h\left(C_t\right).}$

(7) No output activation function (NOAF)

Similar to the NIAF, by removing the activation function of O_t, we obtain a “wilder” version of O_t because O_t is no longer confined to $\left[-1,1\right]$ by the sigmoid function: (12)${\displaystyle O_t=W^O{P}_t+b^O.}$

(8) Coupled input and forget gate (CIFG)

Instead of separately calculating what should be forgotten and what should be inputted as new information, the CIFG combines the two steps. The CIFG forgets only when inputting something in its place and inputs new values to the state only when forgetting something older: (13)${\displaystyle C_t=g\left(W^I{P}_t+b^I\right)\circ C_{t-1}+\left(1-g\left(W^I{P}_t+b^I\right)\right)\circ h\left(W^C{P}_t+b^C\right).}$

(9) Peephole (P)

Gers & Schmidhuber (2000) argued that the cell state should control the gates in order to learn precise timings. Therefore, we add connections from the cell to the gates in Phar-LSTM, which are named as Peephole, to make precise timings easier to learn: (14)${\displaystyle C_t=g\left(W^I\left(P_t,C_{t-1}\right)+b^I\right)\circ C_{t-1}+g\left(W^F\left(P_t,C_{t-1}\right)+b^F\right)\circ h\left(W^C{P}_t+b^C\right),}$ (15)${\displaystyle O_t=g\left(W^O\left(P_t,C_t\right)+b^O\right).}$

(10) Full gate recurrence (FGR)

The LSTM (Van Houdt, Mosquera & Nápoles, 2020) consists of cell state and input and output gates and does not include the forget gate and peephole connections. A hybrid of real-time recurrent learning (Robinson & Fallside, 1987) and backpropagation through time (Werbos, 1988) is used for training. In this case, only the gradient of the cell state was propagated back, and the gradient for the other recurrent connections was truncated. FGR means that all the gates received recurrent inputs from the previous time step: (16)${\displaystyle I_t=g\left(W^I\left[\begin{array}{c}\hfill P_t,\hfill \\ \hfill I_{t-1},\hfill \\ \hfill F_{t-1},\hfill \\ \hfill O_{t-1},\hfill \\ \hfill C_{t-1}\hfill \end{array}\right]+b^I\right),}$ (17)${\displaystyle F_t=g\left(W^F\left[\begin{array}{c}\hfill P_t,\hfill \\ \hfill I_{t-1},\hfill \\ \hfill F_{t-1},\hfill \\ \hfill O_{t-1},\hfill \\ \hfill C_{t-1}\hfill \end{array}\right]+b^F\right),}$ (18)${\displaystyle O_t=g\left(W^O\left[\begin{array}{c}\hfill P_t,\hfill \\ \hfill I_{t-1},\hfill \\ \hfill F_{t-1},\hfill \\ \hfill O_{t-1},\hfill \\ \hfill C_{t-1}\hfill \end{array}\right]+b^O\right).}$

Classification and training

We use the SoftMax classifier for classification. Let k denote the number of DDI types. The output $o\in R^{\left|k\right|}$ is the probabilities of each class to which S belongs. (19)${\displaystyle y=argmax\left(\frac{1}{exp\left(W{H}_n+b\right)}\left[\begin{array}{c}\hfill exp\left(W_1{H}_n+b\right)\hfill \\ \hfill exp\left(W_2{H}_n+b\right)\hfill \\ \hfill \cdots \hfill \\ \hfill exp\left(W_k{H}_n+b\right)\hfill \end{array}\right]\right).}$

We use the cross-entropy (De Boer et al., 2005) cost function and ridge regularization (Hoerl & Kennard, 1970) as the optimization objective. For the i th instance, y⁽ⁱ⁾ denotes the output. The cross-entropy cost is (20)${\displaystyle J=-\left({\sum }_{j=1}^{k}1\left\{y^{\left(i\right)}=j\right\}log\frac{exp\left(W_j{H}_n+b\right)}{{\sum }_{l=1}^{k}exp\left(W{H}_n+b\right)}\right)+\frac{\lambda }{2}{∥W∥}^2,}$ where $1\left\{\blacksquare \right\}$ is the indicator function, such that $1\left\{\text{a true statement}\right\}=1$ and $1\left\{\text{a false statement}\right\}=0$. We optimize the parameters of the objective function J with Rmsprop (Chowdhury & Lavelli, 2011), which is a variant of mini-batch stochastic gradient descent. During each training step, the gradient of J is calculated. Then, all the parameters are adjusted according to the gradient. After the end of training, we have a model that is able to predict two drugs’ interactions when a sentence about these drugs is given.

Evaluation metrics

In this section, we describe the evaluation metrics used in our experiments. For the DDIExtraction 2011 and DDIExtraction 2013 corpuses, Precision (P), Recall (R), F-score (F), and Accuracy (Acc) are widely used as the evaluation metrics (Asada, Miwa & Sasaki, 2023). Since the DDIExtraction 2013 corpus is a multi-class classification problem, we adopt the micro-average and macro-average strategy to score the overall performance on the five classes.

To obtain the F-score, the contingency table (or confusion matrix) is built first, in which each row of the matrix represents the instances in a predicted class and each column represents the instances in an actual class. The contingency table enables us to obtain the true positive (TP), false positive (FP), false negative (FN), and true negative (TN). Based on that, the precision, recall, F-score, and accuracy can be defined as follows: (21)${\displaystyle P=\frac{\text{TP}}{\text{TP}+\text{FP}},}$ (22)${\displaystyle R=\frac{\text{TP}}{\text{TP}+\text{FN}},}$ (23)${\displaystyle F=2\times \frac{P\times R}{P+R},}$ (24)${\displaystyle \text{Acc}=\frac{\text{TP}+\text{TN}}{\text{TP}+\text{TN}+\text{FP}+\text{FN}}.}$

For the DDIExtraction 2013 corpus, there are five P, R, and F-score values for each class since there are five different classes (the five DDI types: ADVICE, EFFECT, MECHANISM, INT, and NEGATIVE). Each DDI type is evaluated separately. Moreover, to measure the overall performance, two commonly used metrics, i.e., micro-averaged F-score (CLA) and macro-averaged F-scores (MAVG), are calculated. The CLA is calculated by constructing a global contingency table and then calculating the precision and recall, and the MAVG is calculated by first calculating the precision and recall for each type and then taking the average of those results.

To evaluate the scalability of our method, we propose a metric to evaluate the performance gap of the models between the two corpuses under the assumption that if a model has good performance in scalability, it would not only have a high average F-score but also have less variance. For instance, the F-scores of model A are 0.65 and 0.65 on the DDIExtraction 2011 and 2013 corpuses, respectively, and those of model B are 0.60 and 0.70. Although the average F-scores of model A and model B are both 0.65, model A would be considered to have better scalability than model B. Based on this, the metric can be defined using the 1-standard deviation of F-scores as (25)${\displaystyle 1-\sigma =1-\sqrt{\begin{array}{c}\hfill \frac{1}{2}{\left(F_{2011}-\frac{F_{2011}+F_{2013}}{2}\right)}^2+\frac{1}{2}{\left(F_{2013}-\frac{F_{2011}+F_{2013}}{2}\right)}^2\hfill \end{array}}.}$

To evaluate the consistency, the training process for 200 epochs of each learning model is shown as a boxplot, and Welch’s t-test at a significance level of alpha = 0.05 was used to determine whether the mean test set performance of a learning model was significantly different from that of Phar-LSTM.

To evaluate the reproducibility, we first ran the training process for 10 times using different random seeds and obtained the boxplot of the overall training process to show the gap between each run. Based on the boxplot, we further calculated the sum of the variance of the F-scores of each epoch and the sum of the standard deviation of the F-scores of each epoch to measure the differences among the 10 runs.

Hyperparameter settings

Based on previous research and experience, the Phar-LSTMs were trained by an RMSprop optimizer with a loss function of cross entropy and a learning rate of 0.001. Dropout layers were added to each of the embedding layers and hidden layers with a ratio of 0.2. For each run, the number of training epochs were 200 and the batch sizes were 32. All the experiments were run on GeForce GTX-1080 and took 9.3 h on average to complete.

Scalability

To evaluate the scalability of our method, experiments were conducted on both the DDIExtraction 2011 and DDIExtraction 2013 datasets. Before 2013, most studies were evaluated on the DDIExtraction 2011 dataset. After 2013, most research has focused on evaluating the methods on the DDIExtraction 2013 dataset. As far as we know, few methods (Björne, Kaewphan & Salakoski, 2013; Jiang, Gu & Jiang, 2017; Thomas et al., 2013) have been evaluated on both datasets.

We first compared the performance of our scheme with the traditional methods as well as the deep learning-based method on the DDIExtraction 2011 dataset. The results are shown in Table 1. The traditional methods are typically based on manually extracted features or kernels. For example, Björne et al. (2011b) leveraged many syntactic-based features, including tokens, dependency types, POS tags, text, and stems. Similarly, Thomas et al. (2013) combined an all-path-graph kernel, a shallow linguistic kernel, and a k-band shortest path spectrum kernel, which were all derived from syntactic analysis. Other methods such as FBK-HLT (Chowdhury et al., 2011) and LIMSI-FBK (Björne et al., 2011a) used either features or kernels or both.

These features and kernels are highly dependent on third-party tools such as syntactic parsing, which makes the method sensitive to the quality of the parsing results and the expertise of researchers in designing features or kernels. Therefore, although the heuristic idea of using features and kernels can be helpful to other researchers, the models themselves may not have good scalability.

It can be observed from Table 1 that the Phar-LSTM scheme achieved the best performance, with the F-score of 0.6880. Another deep learning-based method, Skeleton-LSTM (Jiang, Gu & Jiang, 2017) (F-score: 0.6769), also performed significantly better than other traditional methods, which illustrates the superiority of deep learning-based methods for the coarse-grained task of DDI extraction.

Table 2 shows the results of our scheme in comparison with the baselines on the DDIExtraction 2013 dataset. From Table 2, we can observe that Phar-LSTM achieved the best performance in terms of MAVG and CLA (0.708 and 0.716, respectively). Skeleton-LSTM (Jiang, Gu & Jiang, 2017) had similar performance as Phar-LSTM and performed significantly better than other baselines. One posssible reason may be that Skeleton-LSTM and Phar-LSTM both use the end-end-learning framework (i.e., feed the raw input into the neural network and produce the output directly), which would capture some latent features because the features are automatically extracted by the neural network rather than by third-party tools.

To further evaluate the scalability of our scheme, three baselines were chosen for comparing with our defined metric, and the results are shown in Table 3. Note that many methods developed for addressing the coarse-grained DDI extraction task may not be applicable for the fine-grained task of DDI extraction. It can be observed in Table 3 that the Phar-LSTM scheme achieved the highest (1 − σ) value (0.986) among all methods, which demonstrates the scalability of our method.

Consistency

To evaluate the consistency of the Phar-LSTM scheme, we compared Phar-LSTM with 10 different variants of LSTM, in which the number of epochs for each variant were set to 200. Welch’s t-test at a significance level of alpha = 0.05 was used to determine whether the performance of each variant was significantly different from another. A summary of the results of the different methods with 200 epochs is shown in Fig. 5. The boxplots of the variants that differ significantly from Phar-LSTM are highlighted in blue.

It can be observed from Fig. 5 that Phar-LSTM generally achieved the best performance on both datasets. Moreover, the F-scores of Phar-LSTM for most epochs were relatively stable, which indicates the consistency of our method. Another observation based on Fig. 5 is that removing the output gate (NOG) or the activation functions (NOAF, NIAF, and NFAF) significantly hurt the performance on the two datasets. The ability to output information and the activation of the perceptron appear to be critical for the LSTM architecture. This is probably because the output value of the hidden layers cannot be constrained without the activation function and therefore fails to train the parameters. If the output gate is removed from the LSTM unit, although $x_t^{\left(\text{token}\right)}$ and $x_t^{\left(\text{Phar}\right)}$ can be integrated to the hidden layer H_t by C_t, the original information of and $x_t^{\left(\text{token}\right)}$ and $x_t^{\left(\text{Phar}\right)}$ is diluted during the calculation of H_t. And then, the final SoftMax regression built on the hidden layer of the last unit captures little information of the input, which leads to the failure of training.

On the contrary, although removing the input gate (NIG) or forget gate (NFG) or coupling them into one gate (CIFG) can decrease the F-score, comparing with the NOG, the input information is still integrated through the other gate (e.g., the forget gate for the NIG and the input gate for the NFG). Therefore, the parameters can be trained successfully.

Similarly, removing the pharmacological gate (NPG) generally decreases the F-score more than the NIG, NFG, and CIFG do. This illustrates that the pharmacological gate contains important information to represent the DDI than other gates do. This proves that the Phar-LSTM scheme indeed improved the DDI extraction performance by incorporating the pharmacological gate.

Both adding the Peephole (P) and the full gate recurrence (FGR) decrease the performance while increasing the computational complexity. We generally advise against using them for DDI extraction.

Reproducibility

Due to the random seed mechanism and the implementation of the GPU training architecture, the training process is usually unreproducible. To evaluate the reproducibility of Phar-LSTM, we ran our scheme for many times to check the differences of the outputs. The boxplots of Phar-LSTM’s performance for 10 runs with different epochs on the two datasets are shown in Fig. 6, from which we can see that the performances of Phar-LSTM for most epochs are close.

Some specific epochs can be observed for DDIExtraction 2011, such as the 19th, 37th, 55th, 56th, and 72th epochs. For these specific epochs, the performance of different runs differed from each other. However, there were no such specific epochs for DDIExtraction 2013. One possible reason may be the data distribution. DDIExtraction 2011 is smaller, and the annotation strategy is different from that of DDIExtraction 2013. Researchers should be aware of these specific epochs. The best way to check the model is using a validation set. By conducting experiments on the validation set and drawing the boxplot of the learning curve of different runs, researchers can easily find these specific epochs and improve their extraction system.

Another finding is that both the climb stages (0 < epoch < 10) of DDIExtraction 2011 and DDIExtraction 2013 blurred. The reason is that the initial random states of the parameters are different, which may cause the performances to differ during the climb stage. However, after the climb stage, the performances with the following epochs are much closer. This means that the Phar-LSTM scheme can adapt to different initial random seeds.

To further compare the reproducibility of Phar-LSTM with other variants, we summed the variance and the standard deviation of each epoch for the 10 runs. The metrics (variance and standard deviation) can indicate the overall reproducibility. From Table 4, we can see that Phar-LSTM, NOG, NIAF, and NFAF had better scores than the other models, and the conclusion is consistent with Fig. 5. However, from Fig. 5, we can see that the NOG, NIAF and NFAF had poor F-scores. Although the results of the three models can be reproduced easily, the value of the three models is low. Phar-LSTM reaches a good balance between high reproducibility and high F-score.