Computational drug repositioning based on side-effects mined from social media

Mining Twitter for drug side-effects

We used the SoMeDoSEs pharmacovigilance system (Plachouras, Leidner & Garrow, under review) to extract reports of drug side-effects from Twitter over a 6 month period between January and June 2014. SoMeDoSEs works by first applying topic filters to identify tweets that contain keywords related to drugs, before applying volume filters which remove tweets that are not written in English, are re-tweets or contain a hyperlink to a web page, since these posts are typically commercial offerings. Side-effects were then mapped to an entry in the FDA Adverse Event Reporting System. Tweets that pass these filters are then classified by a linear SVM to distinguish those that mention a drug side-effect from those that do not. The SVM classifier uses a number of natural language features including unigrams and bigrams, part-of-speech tags, sentiment scores, text surface features, and matches to gazetteers related to human body parts, side-effect synonyms, side-effect symptoms, causality indicators, clinical trials, medical professional roles, side effect-triggers and drugs.

For each gazetteer, three features were created: a binary feature, which is set to 1 if a tweet contains at least one sequence of tokens matching an entry from the gazetteer, the number of tokens matching entries from the gazetteer, and the fraction of characters in tokens matching entries from the gazetteer. For side-effect synonyms we used the Consumer Health Vocabulary (CHV) (Zeng et al., 2005), which maps phrases to Unified Medical Language System concept universal identifiers (CUI) and partially addresses the issue of misspellings and informal language used to discuss medical issues in tweets. The matched CUIs were also used as additional features.

To develop the system, 10,000 tweets which passed the topic and volume filters were manually annotated as mentioning a side-effect or not, resulting in a Cohen’s Kappa for inter-annotator agreement on a sample of 404 tweets annotated by two non-healthcare professional of 0.535. Using a time-based split of 8,000 tweets for training, 1,000 for development, and 1,000 for testing, the SVM classifier that used all the features achieved a precision of 55.0%, recall of 66.9%, and F1 score of 60.4% when evaluated using the 1,000 test tweets. This is statistically significantly higher than the results achieved by a linear SVM classifier using only unigrams and bigrams as features (precision of 56.0%, recall of 54.0% and F1 score of 54.9%). One of the sources of false negatives was the use of colloquial and indirect expressions by Twitter users to express that they have experienced a side-effect. We also observed that a number of false positives discuss the efficacy of drugs rather than side-effects.

Twitter data

Over the 6 month period, SoMeDoSEs typically identified ∼700 tweets per day that mentioned a drug side-effect, resulting in a data set of 620 unique drugs and 2,196 unique side-effects from 108,009 tweets, once drugs with only a single side-effect were excluded and drug synonyms had been resolved to a common name using exact string matches to entries in World Drug Index (Thomson Reuters, 2015b), which worked for approximately half of the data set with the remainder matched manually. We were also careful to remove indications that were falsely identified as side-effects using drug indications from Cortellis Clinical Trials Intelligence (Thomson Reuters, 2015a). We used this data to construct a 2,196 row by 620 column matrix of binary variables X, where x ∈ {0, 1}, indicating whether each drug was reported to cause each side-effect in the Twitter data set.

Calculating the sample covariance matrix

Using this data, we are able to form the sample covariance matrix S for binary variables as follows (Allen, 1997), such that element S_i,j gives the covariance of drug i with drug j: (1)${\displaystyle \begin{array}{ll}{S}_{i,j}\hfill & =\frac{1}{n-1}\sum _{k=1}^n\left(x_{ki}-\overline{x_i}\right)\left(x_{kj}-\overline{x_j}\right)\hfill \\ \hfill & =\frac{1}{n-1}\sum _{k=1}^n{x}_{ki}{x}_{kj}-\overline{x_i}\overline{x_j}\hfill \end{array}}$ where $\overline{x_i}=\frac{1}{n}\sum _{k=1}^n{x}_{ki}$ and x_ki is the k-th observation (side-effect) of variable (drug) X_i. It can be shown than the average product of two binary variables is equal to their observed joint probabilities such that: (2)${\displaystyle \frac{1}{n-1}\sum _{k=1}^n{x}_{ki}{x}_{kj}=P\left(X_j=1|X_i=1\right)}$ where P(X_j = 1|X_i = 1) refers to the conditional probability that variable X_j equals one given that variable X_i equals one. Similarly, the product of the means of two binary variables is equal to the expected probability that both variables are equal to one, under the assumption of statistical independence: (3)${\displaystyle \overline{x_i}\overline{x_j}=P\left(X_i=1\right)P\left(X_j=1\right).}$ Consequently, the covariance of two binary variables is equal to the difference between the observed joint probability and the expected joint probability: (4)${\displaystyle S_{i,j}=P\left(X_j=1|X_i=1\right)-P\left(X_i=1\right)P\left(X_j=1\right).}$ Our objective is to find the precision or concentration matrix θ by inverting the sample covariance matrix S. Using θ, we can obtain the matrix of partial correlation coefficients ρ for all pairs of variables as follows: (5)${\displaystyle {\rho }_{i,j}=-\frac{{\theta }_{i,j}}{\sqrt{{\theta }_{i,i}{\theta }_{j,j}}}.}$ The partial correlation between two variables X and Y given a third, Z, can be defined as the correlation between the residuals R_x and R_y after performing least-squares regression of X with Z and Y with Z, respectively. This value, denotated ρ_x,y|z, provides a measure of the correlation between two variables when conditioned on the third, with a value of zero implying conditional independence if the input data distribution is multivariate Gaussian. The partial correlation matrix ρ, however, gives the correlations between all pairs of variables conditioning on all other variables. Off-diagonal elements in ρ that are significantly different from zero will therefore be indicative of pairs of drugs that show unique covariance between their side-effect profiles when taking into account (i.e., removing) the variance of side-effects profiles amongst all the other drugs.

Shrinkage estimation

For the sample covariance matrix to be easily invertible, two desirable characteristics are that it should be positive definite, i.e., all eigenvalues should be distinct from zero, and well-conditioned, i.e., the ratio of its maximum and minimum singular value should not be too large. This can be particularly problematic when the sample size is small and the number of variables is large (n < p) and estimates of the covariance matrix become singular. To ensure these characteristics, and speed up convergence of the inversion, we condition the sample covariance matrix by shrinking towards an improved covariance estimator T, a process which tends to pull the most extreme coefficients towards more central values thereby systematically reducing estimation error (Ledoit & Wolf, 2003), using a linear shrinkage approach to combine the estimator and sample matrix in a weighted average: (6)${\displaystyle S^{\prime }=\alpha T+\left(1-\alpha \right)S}$ where α ∈ {0, 1} denotes the analytically determined shrinkage intensity. We apply the approach of Schäfer and Strimmer, which uses a distribution-free, diagonal, unequal variance model which shrinks off-diagonal elements to zero but leaves diagonal entries intact, i.e., it does not shrink the variances (Schäfer & Strimmer, 2005). Shrinkage is actually applied to the correlations rather than the covariances, which has two distinct advantages: the off-diagonal elements determining the shrinkage intensity are all on the same scale, while the partial correlations derived from the resulting covariance estimator are independent of scale.

Graphical lasso for sparse inverse covariance estimation

A useful output from the covariance matrix inversion is a sparse ρ matrix containing many zero elements, since, intuitively, we know that relatively few drug pairs will share a common mechanism of action, so removing any spurious correlations is desirable and results in a more parsimonious relationship model, while the non-zero elements will typically reflect the correct positive correlations in the true inverse covariance matrix more accurately (Jones et al., 2012). However, elements of ρ are unlikely to be zero unless many elements of the sample covariance matrix are zero. The graphical lasso (Friedman, Hastie & Tibshirani, 2008; Banerjee, ElGhaoui & d’Aspremont, 2008; Hastie, Tibshirani & Wainwright, 2015) provides a way to induce zero partial correlations in ρ by penalizing the maximum likelihood estimate of the inverse covariance matrix using an ℓ₁-norm penalty function. The estimate can be found by maximizing the following log-likelihood using the block coordinate descent approach described by Friedman, Hastie & Tibshirani (2008): (7)${\displaystyle logdet\theta -\mathrm{tr}\left(S^{\prime }\theta \right)-\lambda \Vert \theta {\Vert }_1.}$ Here, the first term is the Gaussian log-likelihood of the data, tr denotes the trace operator and ‖θ‖₁ is the ℓ₁-norm—the sum of the absolute values of the elements of θ, weighted by the non-negative tuning paramater λ. The specific use of the ℓ₁-norm penalty has the desirable effect of setting elements in θ to zero, resulting in a sparse matrix, while the parameter λ effectively controls the sparsity of the solution. This contrasts with the use of an ℓ₂-norm penalty which will shrink elements but will never reduce them to zero. While this graphical lasso formulation is based on the assumption that the input data distribution is multivariate Gaussian, Banerjee, ElGhaoui & d’Aspremont (2008) showed that the dual optimization solution also applies to binary data, as is the case in our application.

It has been noted that the graphical lasso produces an approximation of θ that is not symmetric, so we update it as follows (Sustik & Calderhead, 2012): (8)${\displaystyle \theta \leftarrow \frac{\left(\theta +{\theta }^T\right)}{2}.}$ The matrix ρ is then calculated according to Eq. (5), before repositioning predictions for drug i are determined by ranking all other drugs according to their absolute values in ρ_i and assigning their indications to drug i.

Recovering known indications

To evaluate our method we have attempted to predict repositioning targets for indications that are already known. If, by exploiting hindsight, we can recover these, then our method should provide a viable strategy with which to augment existing approaches that adopt an integrated approach to drug repositioning (Emig et al., 2013). Figure 1A shows the performance of the method at identifying co-indicated drugs at a range of λ values, resulting in different sparsity levels in the resulting ρ matrix. We measured the percentage at which a co-indicated drug was ranked amongst the top 5, 10, 15, 20 and 25 predictions for the target drug, respectively. Of the 620 drugs in our data set, 595 had a primary indication listed in Cortellis Clinical Trials Intelligence, with the majority of the remainder being made up of dietary supplements (e.g., methylsulfonylmethane) or plant extracts (e.g., Agaricus brasiliensis extract) which have no approved therapeutic effect. Rather than removing these from the data set, they were left in as they may contribute to the partial correlation between pairs of drugs that do have approved indications.

Results indiciate that the method achieves its best performance with a λ value of 10⁻⁹ where 42.41% (243/595) of targets have a co-indicated drug returned amongst the top 5 ranked predictions (Fig. 1A). This value compares favourably with both a strategy in which drug ranking is randomized (13.54%, standard error ±0.65), and another in which drugs are ranked according to the Jaccard index (28.75%). In Ye, Liu & Wei (2014), a related approach is used to construct a repositioning network based on side-effects extracted from the SIDER database, Meyler’s Side Effects of Drugs, Side Effects of Drugs Annual, and the Drugs@FDA database (Kuhn et al., 2010; Aronson, 2015; Ray, 2014; US Food and Drug Administraction, 2016), also using the Jaccard index as the measure of drug–drug similarity. Here, they report an equivilent value of 32.77% of drugs having their indication correctly predicted amongst the top 5 results. While data sets and underlying statistical models clearly differ, these results taken together suggest that the use of side-effect data mined from social media can certainly offer comparable performance to methods using side-effect data extracted from more conventional resources, while the use of a global statistical model such as the graphical lasso does result in improved performance compared to a pairwise similarity coefficient such as the Jaccard index.

To further investigate the influence of the provenance of the data, we mapped our data set of drugs to ChEMBL identifiers (Gaulton et al., 2012; Bento et al., 2014) which we then used to query SIDER for side-effects extracted from drug labels. This resulted in a reduced data set of 229 drugs, in part due to the absence of many combination drugs from SIDER (e.g. the antidepressant Symbyax which contains olanzapine and fluoxetine). Using the same protocol described above, best performance of 53.67% (117/229) was achieved with a slightly higher λ value of 10⁻⁶. Best performance on the same data set using side-effects derived from Twitter was 38.43% (88/229), again using a λ value of 10⁻⁹, while the randomized strategy achieved 12.05% (standard error ± 1.14), indicating that the use of higher quality side-effect data from SIDER allows the model to achieve better performance than is possible using Twitter data. Perhaps more interestingly, combining the correct predictions between the two datasources reveals that 30 are unique to the Twitter model, 59 are unique to the SIDER model, with 58 shared, supporting the use side-effect data mined from social media to augment conventional resources.

We also investigated whether our results were biased by the over-representation of particular drug classes within our data set. Using Using Cortellis Clinical Trials Intelligence, we were able to identify the broad class for 479 of the drugs (77.26%) in our data set. The five largest classes were benzodiazepine receptor agonists (3/14 drugs returned amongst the top 5 ranked predictions), analgesics (6/12), H₁-antihistamines (8/11), cyclooxygenase inhibitors (9/11), and anti-cancer (2/11). This indicates that the over-representation of H₁-antihistamines and cyclooxygenase inhibitors did result in a bias, and to a lesser extent analgesics, but that the overall effect of these five classes was more subtle (28/59 returned amongst the top 5 ranked predictions, 47.46%).

The best performance of our approach at the top 5 level is achieved when the resulting ρ matrix has a sparsity of 35.59% (Figs. 1B and 2) which both justifies the use of the ℓ₁-norm penalized graphical lasso, and generates a graphical model with approximately a third of the parameters of a fully dense matrix, while the comparable performance at λ values between 10⁻¹² and 10⁻⁷ also indicates a degree of robustness to the choice of this parameter. Beyond the top 5 ranked predictions, results are encouraging as the majority of targets (56.02%) will have a co-indicated drug identified by considering only the top 10 predictions, suggesting the method is a feasible strategy for prioritisation of repositioning candidates.

Predicting proposed indications of compounds currently in clinical trials

While the previous section demonstrated our approach can effectively recover known indications, predictions after the fact are—while useful—best supported by more forward-looking evidence. In this section, we use clinical trial data to support our predictions where the ultimate success of our target drug is still unknown. Using Cortellis Clinical Trials Intelligence, we extracted drugs present in our Twitter data set that were currently undergoing clinical trials (ending after 2014) for a novel indication (i.e., for which they were not already indicated), resulting in a subset of 277 drugs currently in trials for 397 indications. Figure 3 shows the percentage at which a co-indicated drug was ranked amongst the top 5, 10, 15, 20 and 25 predictions for the target. Similar to the recovery of known indications, best performance when considering the top 5 ranked predictions was achieved with a λ value of 10⁻⁹, resulting in 16.25% (45/277) of targets having a co-indicated drug, which again compares well to a randomized strategy (5.42%, standard error ± 0.32) or a strategy using the Jaccard index (10.07%).

Recovery of proposed clinical trial indications is clearly more challenging than known indications, possibly reflecting the fact that a large proportion of drugs will fail during trials and therfore many of the 397 proposed indications analysed here will in time prove false, although the general trend in performance as the sparsity parameter λ is adjusted tends to mirror the recovery of known indications. Despite this, a number of interesting predictions with a diverse range of novel indications are made that are supported by experimental and clinical evidence; a selection of 10 of the 45 drugs where the trial indication was correctly predicted are presented in Table 1. We further investigated three repositioning candidates with interesting pharmacology to understand their predicted results.

Oxytocin

Oxytocin is a nonapeptide hormone that acts primarily as a neuromodulator in the brain via the specific, high-affinity oxytocin receptor—a class I (Rhodopsin-like) G-protein-coupled receptor (GPCR) (Gimpl et al., 2002). Currently, oxytocin is used for labor induction and the treatment of Prader-Willi syndrome, but there is compelling pre-clinical evidence to suggest that it may play a crucial role in the regulation of brain-mediated processes that are highly relevant to many neuropsychiatric disorders (Feifel, 2012). A number of animal studies have revealed that oxytocin has a positive effect as an antipsychotic (Feifel & Reza, 1999; Lee et al., 2005), while human trials have revealed that intranasal oxytocin administered to highly symptomatic schizophrenia patients as an adjunct to their antipsychotic drugs improves positive and negative symptoms significantly more than placebo (Feifel et al., 2010; Pedersen et al., 2011). These therapeutic findings are supported by growing evidence of oxytocin’s role in the manifestation of schizophrenia symptoms such as a recent study linking higher plasma oxytocin levels with increased pro-social behavior in schizophrenia patients and with less severe psychopathology in female patients (Rubin et al., 2010). The mechanisms underlying oxytocin’s therapeutic effects on schizophrenia symptoms are poorly understood, but its ability to regulate mesolimbic dopamine pathways are thought to be responsible (Feifel, 2012). Here, our method predicts schizophrenia as a novel indication for oxytocin based on its proximity to chlorpromazine, which is currently used to treat schizophrenia (Fig. 4). Chlorpromazine also modulates the dopamine pathway by acting as an antagonist of the dopamine receptor, another class I GPCR. Interestingly, the subgraph indicates that dopamine also has a high partial correlation coefficient with oxytocin, adding further support to the hypothesis that oxytocin, chlorpromazine and dopamine all act on the same pathway and therefore have similar side-effect profiles. Side-effects shared by oxytocin and chlorpromazine include hallucinations, excessive salivation and anxiety, while shivering, weight gain, abdominal pain, nausea, and constipation are common side-effects also shared by other drugs within the subgraph. Currently, larger scale clinical trials of intranasal oxytocin in schizophrenia are underway. If the early positive results hold up, it may signal the beginning of an new era in the treatment of schizophrenia, a field which has seen little progress in the development of novel efficacious treatments over recent years.

Ramelteon

Ramelteon, currently indicated for the treatment of insomnia, is predicted to be useful for the treatment of bipolar depression (Fig. 5). Ramelteon is the first in a new class of sleep agents that selectively binds the MT₁ and MT₂ melatonin receptors in the suprachiasmatic nucleus, with high affinity over the MT₃ receptor (Owen, 2006). It is believed that the activity of ramelteon at MT₁ and MT₂ receptors contributes to its sleep-promoting properties, since these receptors are thought to play a crucial role in the maintenance of the circadian rhythm underlying the normal sleep-wake cycle upon binding of endogenous melatonin. Abnormalities in circadian rhythms are prominent features of bipolar I disorder, with evidence suggesting that disrupted sleep-wake circadian rhythms are associated with an increased risk of relapse in bipolar disorder (Jung et al., 2014). As bipolar patients tend to exhibit shorter and more variable circadian activity, it has been proposed that normalisation of the circadian rhythm pattern may improve sleep and consequently lead to a reduction in mood exacerbations. Melatonin receptor agonists such as ramelteon may have a potential therapeutic effect in depression due to their ability to resynchronize the suprachiasmatic nucleus (Wu et al., 2013). In Fig. 5, evidence supporting the repositioning of ramelteon comes from ziprasidone, an atypical antipsychotic used to treat bipolar I disorder and schizophrenia (Nicolson & Nemeroff, 2007). Ziprasidone is the second-ranked drug by partial correlation coefficient; a number of other drugs used to treat mood disorders can also be located in the immediate vicinity including phenelzine, a non-selective and irreversible monoamine oxidase inhibitor (MAOI) used as an antidepressant and anxiolytic, milnacipran, a serotonin–norepinephrine reuptake inhibitor used to treat major depressive disorder, and tranylcypromine, another MAOI used as an antidepressant and anxiolytic agent. The co-location of these drugs in the same region of the graph suggests a degree of overlap in their respective mechanistic pathways, resulting in a high degree of similarity between their side-effect profiles. Nodes in this subgraph also have a relatively large degree indicating a tighter association than for other predictions, with common shared side-effects including dry mouth, sexual dysfunction, migraine, and orthostatic hypotension, while weight gain is shared between ramelteon and ziprasidone.

Meloxicam

Meloxicam, a nonsteroidal anti-inflammatory drug (NSAID) used to treat arthritis, is predicted to be a repositioning candidate for the treatment of non-Hodgkin lymphoma, via the mobilisation of autologous peripheral blood stem cells from bone marrow. By inhibiting cyclooxygenase 2, meloxicam is understood to inhibit generation of prostaglandin E₂, which is known to stimulate osteoblasts to release osteopontin, a protein which encourages bone resorption by osteoclasts (Rainsford, Ying & Smith, 1997; Ogino et al., 2000). By inhibiting prostaglandin E₂ and disrupting the production of osteopontin, meloxicam may encourage the departure of stem cells, which otherwise would be anchored to the bone marrow by osteopontin (Reinholt et al., 1990). In Fig. 6, rituximab, a B-cell depleting monoclonal antibody that is currently indicated for treatment of non-Hodgkin lymphoma, is the top ranked drug by partial correlation, which provides evidence for repositioning to this indication. Interestingly, depletion of B-cells by rituximab has recently been demonstrated to result in decreased bone resorption in patients with rheumatoid arthritis, possibly via a direct effect on both osteoblasts and osteoclasts (Wheater et al., 2011; Boumans et al., 2012), suggesting a common mechanism of action between meloxicam and rituximab. Further evidence is provided by the fifth-ranked drug clopidogrel, an antiplatelet agent used to inhibit blood clots in coronary artery disease, peripheral vascular disease, cerebrovascular disease, and to prevent myocardial infarction. Clopidogrel works by irreversibly inhibiting the adenosine diphosphate receptor P2Y12, which is known to increase osteoclast activity (BonekEY Watch, 2012). Similarly to the ramelteon subgraph, many drugs in the vicinity of meloxicam are used to treat inflammation including diclofenac, naproxen (both NSAIDs) and betamethasone, resulting in close association between these drugs, with shared side-effects in the subgraph including pain, cramping, flushing and fever, while swelling, indigestion, inflammation and skin rash are shared by meloxicam and rituximab.

While the side-effects shared within the subgraphs of our three examples are commonly associated with a large number of drugs, some of the side-effects shared by the three drug pairs such as hallucinations, excessive salivation and anxiety are somewhat less common. To investigate this relationship for the data set as a whole, we calculated log frequencies for all side-effects and compared these values against the normalized average rank of pairs where the side-effect was shared by both the query and target drug. If we assume that a higher ranking in our model indicates a higher likelihood of drugs sharing a protein target, this relationship demonstrates similar properties to the observations of Campillos et al. (2008) in that there is a negative correlation between the rank and frequency of a side-effect. The correlation coefficient has a value of −0.045 which is significantly different from zero at the 0.001 level, although the linear relationship appears to break down where the frequency of the side-effect is lower than about 0.025.