SafeDrug: Dual Molecular Graph Encoders for Recommending Effective and Safe Drug Combinations

1 Introduction

Today abundant health data such as longitudinal electronic health records (EHR) and massive biomedical data available on the web enable researchers and doctors to build better predictive models for clinical decision making [Choi et al., 2016a; Xiao et al., 2018]. Among others, recommending effective and safe medication combinations is an important task, especially for helping patients with complex health conditions [Shang et al., 2019; Shang et al., 2018]. The primary goal of recommending medication combination is to customize a safe combination of drugs for a particular patient based on the patient's health conditions. Earlier medication recommendation models are instance-based, which are only based on the current hospital visit [Zhang et al., 2017; Wang et al., 2017]. As a result, a patient with newly diagnosed hypertension will likely be treated the same as another patient who has suffered chronic uncontrolled hypertension. Such a limitation affects the safety and utility of the recommendations. To overcome this issue, Longitudinal methods such as [Le et al., 2018; Shang et al., 2019] are proposed. They leverage the temporal dependencies within longitudinal patient history to provide a more personalized recommendation. However, they still suffer from the following limitations.

• Inadequate Medication Encoding. Existing works [Zhang et al., 2017; Shang et al., 2019] often represent medications using one-hot encoding. Each drug is considered a (binary) unit, ignoring that drugs in their meaningful molecular graph representation encode important drug properties, such as efficacy and safety profiles. Also, molecule substructures are correlated with functionalities. Such knowledge can be useful for improving accuracy and safety in medication recommendations.
• Implicit and Non-controllable DDI Modeling. Some existing works model drug-drug interactions (DDI) via soft or indirect constraints, like knowledge graphs (KGs) [Wang et al., 2017; Mao et al., 2019] and reinforcement postprocessing [Zhang et al., 2017]. These implicit handling of DDIs results in non-controllable rates in the final recommendation or sub-optimal recommendation accuracy.

To address these, we propose a DDI-controllable drug recommendation model, named SafeDrug, to leverage drugs' molecule structures and explicitly model DDIs. We argue it is beneficial to incorporate molecule structures in medication recommendation. Our SafeDrug has the following contributions:

• Dual Molecular Encoders to Capture Global and Local Molecule Patterns. SafeDrug model firstly learns patient representation, which is fed into dual molecule encoders to capture the global pharmacological properties [Brown and Fraser, 1868] and the local sub-structural patterns of a drug [Huang et al., 2020b]. Globally, a message-passing neural network (MPNN) encoder is constructed to pass molecular information messages for the whole drug layer by layer. Drug connectivity information is well captured in multiple hops. Locally, a bipartite encoder segments drug molecules into substructures, each possibly associated with

small functional groups. In this work, the substructure representation is fed into an effective masked neural network. The final output of the model is obtained by element-wise integration of the global and local encoded embeddings.

• DDI Controllable Loss Function. Inspired by proportional-integral-derivative (PID) controller [Astrom and Hagglund, 1995], we design a new technique to adaptively combine supervised loss and unsupervised DDI constraints with considerable flexibility. During the training, the negative DDI signal would be emphasized and backpropagated if the DDI rate of individual samples is above a certain threshold/target. In the experiment, the adaptive gradient descending can balance the model accuracy and final DDI level. With a preset target, SafeDrug model can provide reliable drug combinations to meet different levels of DDI requirements.
• Comprehensive Evaluations. We follow previous works [Zhang et al., 2017; Shang et al., 2019; Le et al., 2018] and evaluate SafeDrug on a benchmark medication recommendation dataset extracted from the MIMIC-III database [Johnson et al., 2016]. SafeDrug outperforms the best baselines relatively by 19.43% improvement in DDI reduction, 2.88% in Jaccard similarity, and 2.14% in F1 measure. Besides, SafeDrug requires much fewer parameters than previous deep learning based drug recommendation approaches with 14% reduction in training time and around 2× speed-up during inference.

Data processing files and codes are released in Github 1 . A full version of the paper can be found in arXiv 2 .

2 RelatedWorks

2.1 Medication Recommendation

Existing medication recommendation algorithms can be categorized into instance-based and longitudinal recommendation methods. Instance-based approaches focus on the patient's current health status. For example, LEAP [Zhang et al., 2017] extracts feature information from the current encounter and adopts a multi-instance multi-label drug recommendation setting. Longitudinal approaches are proposed to leverage the temporal dependencies within the clinical history, see [Choi et al., 2016a; Xiao et al., 2018]. Among them, RETAIN [Choi et al., 2016b] utilizes a two-level RNN with reverse time attention to model the longitudinal information. GAMENet [Shang et al., 2019] adopts memory neural networks and stores historical drug information as references for further prediction. These methods either model the final drug combination as multi-label binary classification [Choi et al., 2016b; Shang et al., 2019; Shang et al., 2018] or by sequential decision making [Zhang et al., 2017]. Despite their initial success, these existing works still have the following limitations. Other important data such as drug molecule structures can help augment drug recommendation but have not been utilized by existing models. Also, existing works model DDI inadequately, either via soft or indirect constraints.

1 https://github.com/ycq091044/SafeDrug

2 https://arxiv.org/abs/2105.02711

This paper proposes a DDI-controllable drug recommendation model that leverages drugs' molecule structures and models DDIs effectively.

2.2 Molecule Representation

Molecule graph representation remains an important topic due to the association between molecule structure and properties, e.g., efficacy and safety. Early days, molecular descriptors [Mauri et al., 2006] and drug fingerprint [Rogers and Hahn, 2010; Duvenaud et al., 2015] were commonly used to represent drug molecules. Deep learning models were developed in recent years to learn molecule representation and model molecule substructures (a set of connected atoms). For example, [Huang et al., 2020b] developed a sequential pattern mining method on SMILES strings and modeled a pair of drugs as a set of their substring representation. Also, [Huang et al., 2020a] proposed to directly model molecule graphs using graph-based neural network models.

In this paper, inspired by [Huang et al., 2020a], SafeDrug is proposed to capture both the molecular global and local information during the recommendation.

3 ProblemFormulation

Electrical Health Records (EHR). Patient EHR data captures comprehensive medical histories of patients in the format of longitudinal vectors of medical codes (e.g., diagnoses, procedures, and drugs). Formally, EHR for patient j can be represented as a sequence Xj = [x(1) j , x(2) j , . . . , x(Vj) j ], where Vj is the total number of visits for patient j. The i-th entry of Xj documents the i-th visit of patient j, specified by x(i) j = [d(i) j , p(i) j ,m(i) j ], where d(i) j ∈ {0, 1}|D| , p(i) j ∈ {0, 1}|P| , m(i) j ∈ {0, 1}|M| , are multi-hot diagnosis, procedure and medication vectors, respectively. D,P,M are the corresponding element sets, and | · | denotes the cardinality.

Safe Drug Combination Recommendation. We focus on drug recommendation while reducing/controlling the drugdrug interaction (DDI) rates in the prediction. We use a symmetric binary adjacency matrix D ∈ {0, 1}|M|×|M| to denote DDI relation. Dij = 1 means that the interaction between drug i and j has been reported while Dij = 0 indicates a safe co-prescription. The paper aims at learning a drug recommendation function, f(·), such that for each visit, e.g., the t-th visit, given the patient's up-to-now diagnosis and procedure sequences [d(1) , d(2) , . . . , d(t) ] and [p(1) , p(2) , . . . , p(t) ], along with the DDI matrix D, f(·) can generate a drug recommendation,

The objective function consists of two parts: (i) extracting real drug combinations, m(t) , as supervision to penalize ˆm(t) ; (ii) using matrix D to derive unsupervised DDI constraints on ˆm(t) .

4 The SafeDrugModel

As illustrated in Figure 1, our SafeDrug model comprises of four components: (1) a longitudinal patient representation module that learns patient representation from their EHR data; (2) a global message passing neural network (MPNN) encoder that takes patient representation as input and outputs a global drug vector with each entry quantifying the similarity between the patient representation and each of the drug representations; (3) in parallel, a bipartite encoder also takes the same patient representation and outputs a local drug vector, which encodes drugs' molecular substructural functionalities; (4) finally, the global and local drug representation vectors are element-wise combined in the medication representation module, where the final drug output is obtained from element-wise thresholding.

4.1 Longitudinal Patient Representation

From longitudinal EHR data, patient health can be encoded by their diagnosis and procedure information.

Diagnosis Embedding. To fully utilize rich diagnosis information, we design an embedding table, Ed ∈ R|D|×dim, where each row stores the embedding vector for a particular diagnosis and dim is the dimension of the embedding space. Given a multi-hot diagnosis vector d(t) ∈ {0, 1}|D| , we project the corresponding diagnosis codes into the embedding space by vector-matrix dot product, which essentially takes the summation of each diagnosis embedding,

During the training, Ed is learnable and shared between each visit and every patient.

Procedure Embedding. Similarly, a shared procedure table Ep ∈ R|P|×dim is also designed to encode the associated procedure vector p(t) (also a multi-hot vector),

Two embedding vectors, d(t) e , p(t) e ∈ Rdim, encode the patient's current health condition. In fact, a health snapshot might be insufficient for treatment decisions. For example, it is inappropriate to prescribe the same medicine for one patient with a newly diagnosed hypertension and another patient who has chronic hypertension for more than two years.

Here we use RNN models to model dynamic patient history. Specifically, we use two separate RNNs to obtain the hidden diagnosis and procedure vectors, d(t) h and p(t) h , in case that one of the sequences might be inaccessible in practice,

where d(t) h , p(t) h ∈ Rdim. And d(0) h , p(0) h are all-zero vectors. Patient Representation. We then concatenate the diagnosis embedding d(t) h and procedure embedding p(t) h into a more compact patient representation h(t) . We follow a common and effective approach to first concatenate two vectors as a double-long vector, and then apply a feed-forward neural network, NN1(·) : R2dim ↦→ Rdim, (W1 is the parameter, and # is the concatenation operation),

To this end, we obtain an overall patient representation h(t) ∈ Rdim, from which, we next generate a safe drug recommendation by comprehensively modeling the drug molecule databases.

4.2 Dual Molecular Graph Encoders

To investigate drug properties and their dependencies, we exploit the global molecule connectivity features and the local molecule functionalities. In the following sections, we introduce two new molecular graph encoders, in parallel, to comprehensively encode the drug molecules.

(I) Global MPNN Encoder

We encode drug molecule data using a message passing neural network (MPNN) operator with learnable fingerprints, aiming at convolving and pooling atom information across single molecule graphs into vector representations.

To begin with, we collect a set of all appeared atoms, C = {ai}, and design a learnable atom embedding table, Ea ∈ R|C|×dim, where each row looks up the initial embedding/fingerprint for a particular atom. Our design of the MPNN model is very different from previous works [Coley et al., 2017; Yang et al., 2019], which mostly use atom descriptors, e.g., atommass and chirality, as the initial features. However, for molecule structure modeling, atom-atom connectivity is more important than individual atoms.

Given a drug molecule graph (atoms a0, . . . , an as the vertices and atom-atom bonds as edges) with adjacency matrix A and the initial atom fingerprints from Ea, i.e., {y(0) 0 , y(0) 1 , . . . , y(0) n }, we formulate the layer-wise message passing over the graph as,

where l is the layer index, z(l+1) i is the encoded messages from the neighbors of atom i at the l-th iteration, y(l+1) i is the hidden state of atom i. During message passing, z(l+1) i and y(l+1) i with each vertex atom i are updated using message function, MESSAGEl(·), and vertex update function, UPDATEl(·), and W(l) is the layer-wise parameter matrix. After applying message passing for L layers, the global representation of a drug molecule is pooled by a readout function, which calculates the average of all atom fingerprints,

We apply the same MPNN encoder with shared parameters for each single drug moledule (there are |M| different drugs in total). We collect the MPNN embeddings of all drug molecules into a drug memory Eg ∈ R|M|×dim, where each row is the embedding of a drug.

Patient-to-drug Matching. When treating patient representation h(t) as a query, we aim at searching the most relevant drugs from the memory Eg. The health-to-drug matching score is calculated by the dot product of the MPNN representation and the patient representation h(t) . For all drugs, we operate each row of Eg with h(t) , followed by a sigmoid function σ1(·),

Every element of m(t) r stores a matching score for one drug. In light of the recent successful post-LN architecture [Wang et al., 2019] in Transformer [Vaswani et al., 2017], the matching score is further parameterized by a feed-forward neural network, NN2(·) : Rdim ↦→ Rdim, followed by layernormalization (LN) step with residual connection,

(II) Local Bipartite Encoder

Molecules with similar structures shall be mapped into nearby space using the MPNN encoder. However, a pair of drugs that have significant overlap on structure domain might behave differently in terms of DDI interaction or other functional activities [Ryu et al., 2018]. In fact, a drug's functionality is better reflected by molecule substructures, which are defined as a set of connected atoms that may have associated functional groups. This section focuses on the substructure composition of molecules and builds a substructure-to-drug bipartite architecture for capturing drug local functionality and their dependencies.

Molecule Segmentation. To capture the relations between drug molecules and their substructures, we first adopt the well-known molecule segmentation method, breaking retrosynthetically interesting chemical substructures (BRICS) [Degen et al., 2008], which preserves the most critical drug functional groups and more crucial bonds. SafeDrug can easily support other segmentation methods such as RECAP, ECFP. After the segmentation, each molecule is decomposed into a set of chemical substructures. In this setting, a pair of drugs could intersect on a mutual set of substructures. We denote the overall substructure set as S. Based on the many-tomany relations, it is ready to construct a bipartite architecture, represented by a mask matrix H ∈ {0, 1}|S|×|M| : Hij = 1 indicates that drug j contains substructure i.

Bipartite Learning. We hope that H can help to derive a drug representation on its functionality information, given the input, h(t) . To do so, we first apply a feed-forward neural network, NN3(·) : Rdim ↦→ R|S| , for dimensionality transformation, followed by a sigmoid function, σ2(·). In fact, NN3(·) transforms the patient representation h(t) into a local functional vector,

Each entry of m(t) f quantifies the significance of the corresponding functionality. Therefore, m(t) f can be viewed as an anti-disease functionality combination for the current patient. Next, we are motivated to generate the recommendation of drugs that can cover all these anti-disease functionalities while considering to prevent DDIs.

We therefore perform network pruning and design a masked 1-layer neural network, NN4(·) : R|S| ↦→ R|M| , where the parameter matrix is masked by the bipartite architecture H, essentially through matrix element-wise product ⊙. During the training process, NN4(·) will learn to map the functionality vector m(t) f into the corresponding local drug representation,

We show that the mask H would enable the model to (i) have far fewer parameters (due to the sparsity of H); (ii) prevent DDIs by avoiding a co-prescription of interacted drugs (we demonstrate it in Appendix).

4.3 Medication Representation

The final medication representation is given by treating the global drug matching vector as attention signals to further adjust the local drug functionality representation. We finally use a sigmoid function σ3(·) to scale the output,

where ⊙ represents the element-wise product. By a threshold value δ, we can achieve the recommended drug combination, i.e., a multi-hot vector ˆm(t) , by picking those entries where the value is greater than δ.

4.4 Model Training and Inference

SafeDrug is trained end-to-end. We simultaneously learn Ed, Ep, parameters in RNNd and RNNp, network parameters W1, Ea, layer-wise parameters in MPNN {W(i) }, parameters in layer normalization (LN), andW2,W3 andW4.

Multiple Loss. In this paper, the recommendation task is formulated as multi-label binary classification. Suppose |M| is the total number of drugs. We use m(t) to denote the target multi-hot drug recommendation, ˆo(t) for the output realvalued drug representation (opposed to the predicted multihot vector ˆm(t) after applying threshold δ). We consider the prediction of each drug as a sub-problem and use binary cross-entropy (BCE) as a loss,

where m(t) i and ˆo(t) i are the i-th entries. To make the result more robust, we also adopt multi-label hinge loss to ensure that the truth labels have at least 1 margin larger than others,

Based on the DDI adjacency matrix D, we also design the adverse DDI loss over ˆo(t) as,

where · is the product between scalars. Note that the above loss functions are defined for a single visit. During the training, loss back propagation will be conducted at patient level by the averaged losses across all visits.

Controllable Loss Function. A standard approach to training with multiple loss functions is by the weighted sum of the loss measuring terms [Dosovitskiy and Djolonga, 2019],

where α and β are usually pre-defined hyperparameters.

In this scenario, two prediction loss, Lbce and Lmulti, are compatible, and thus we select α from the validation set. We observe that adverse DDIs also exist in the dataset (doctors might mistakenly prescribe the interacted drugs). Therefore, by training with ground truth labels, the correct learning process might conversely increase DDI. Thus, we consider adjusting β during the training process by a proportional-integral-derivative (PID) controller [An et al., 2018]

to balance the prediction loss and DDI loss. For simplicity, we only utilize the proportional error signal as negative feedback when the DDI rate of the recommended drugs is above certain thresholds.

In safe drug recommendation, the key is to maintain a low DDI, so we set a DDI acceptance rate γ ∈ (0, 1) for the loss function. Our motivation is that if the patient-level DDI is lower than the threshold, γ, then we consider only maximizing the prediction accuracy; otherwise, β will adjust adaptively to reduce DDI as well.

The controllable factor β follows the piece-wise strategy,

where Kp is a correcting factor for the proportional signal.

We show in the experiment that γ can be viewed as an upper bound of the output DDI rate. By presetting γ, our SafeDrug can meet different level of DDI requirements.

The model inference phase basically follows the same pipeline as training. We apply a threshold δ = 0.5 on the output drug representation in Eqn. (14) and select the drugs, corresponding to entries where the value is greater than δ, as the final recommendations. We structure the overall algorithm in Appendix.

5 Experiments

Dataset and Metrics. The experiments are carried out on MIMIC-III [Johnson et al., 2016]. The statistics of the postprocessed data is reported in Table 1. We use five efficacy metrics: DDI rate, Jaccard similarity, F1 score, PRAUC, and # of medications, to evaluate the recommendation efficacy and three metrics: # of parameters, training time, and inference time, to evaluate the complexity of the models. Due to space limitation, we push other details to Appendix.

Baselines. We compare SafeDrug with the following baselines from different perspectives: standard Logistic Regression (LR), multi-label classification approach: Ensemble Classifier Chain (ECC), [Read et al., 2009], RNN-based approach: RETAIN [Choi et al., 2016b], instance-based approach: LEAP [Zhang et al., 2017], longitudinal memorybased approaches: DMNC [Le et al., 2018] and GAMENet [Shang et al., 2019]. [Shang et al., 2018] is not considered a baseline since it requires extra ontology data. We also compared SafeDrug with its two model variants: SafeDrugL is for the model with only bipartite encoder (local, L) and SafeDrugG for the model with only the MPNN encoder

(global, G). In implementation, we use m(t) l and m(t) g as the final drug representation for them, respectively.

Performance Comparison. From Table 2, it is easy to see the drug sequence generation models LEAP, DMNC, and ECC yield poor results among the baselines. We , therefore, conclude that formulating the drug recommendation task as multi-label prediction might be more straightforward and effective. LR and RETAIN do not consider the DDI information, so they output undesirable DDI rate. The recommendations given by GAMENet are based on historical combinations, however, real medication records usually contain high DDI rate, e.g., around 8% in MIMIC-III, so it also provides a high DDI rate. Moreover, the performance of two variants, SafeDrugL and SafeDrugG, show that both of the molecular encoders do help with the prediction. Overall, our SafeDrug also outperforms the other baselines with significantly lower DDI rate and better accuracy.

Significance Test. The above results are given by 10 rounds of bootstrapping sampling. Next, we conduct T-test on each metric for two-tailed hypothesis. The p-values are shown in the parenthesis. Commonly, a p-value smaller than 0.05 would be considered as significant. We also use (--) for cases where the result of SafeDrug is worse than the baseline. To conclude, our model can beat all compared baselines in almost every metric with p=0.001 significance level.

Comparison of Model Efficiency. Model complexity of SafeDrug and other DL baselines are evaluated in Table 3. Since RETAIN is not specialized for drug recommendation, improvements in Table 3 are over LEAP, DMNC, and GAMENet. Efficiency wise, our model has relatively lower space and time complexity over LEAP, DMNC, and GAMENet. Also, it is much more efficient during inference. Leap and DMNC adopt sequential modeling and recommend drugs one by one, which is time-consuming. GAMENet stores a large memory bank, and thus it requires larger space. By comparison, we conclude that SafeDrug is efficient and flexible for real deployment.

Analysis of DDI Controllability. We evaluate DDI controllability and show that the average DDI rate of our output can be controlled by threshold γ. The ground truth DDI rate in MIMIC-III is 0.0808. We consider that learning-based drug recommendation approaches should optimize the final recommendations by processing EHR and drug-drug relation databases, so the output DDI rate is supposed to be lower than human experts'. In this section, we test for a series of targeted DDI γ ranging from 0 to 0.08. For each γ, we train a separate model. The experiment is conducted for 5 times. After convergence, we show the average metrics in Table 4.

From the result, the DDI rates are controlled and upperbounded by γ in most cases, which is consistent with the design. When γ becomes larger, more drugs are allowed in one combination, and SafeDrug is more accurate. When γ is too small (< 0.02), model accuracy will drop dramatically. The parameter γ provides a way for doctors to control the tradeoff between DDI rates and accuracy in recommendations.

6 Conclusion

In this paper, we propose SafeDrug accurate and safe drug recommendation. Our model extracts and encode rich molecule structure information to augment the drug recommendation task. We design an adaptive loss controller to provide additional flexibility. We evaluated SafeDrug using benchmark data and showed better accuracy and efficiency.

Acknowledgements

This work was in part supported by the National Science Foundation award SCH-2014438, PPoSS 2028839, IIS1838042, the National Institute of Health award NIH R01 1R01NS107291-01 and OSF Healthcare.