STRING (Szklarczyk et al., 2023) is a database of known and predicted PPIs. We used STRING to get the PPI network of S1R, as shown in Figure 2A; we marked the correlation scores of proteins related to S1R in the network, among which the scores of dopamine D2 receptor (DRD2) and binding-immunoglobulin protein (BIP) are highest, 0.983 and 0.990, respectively, so we picked them as primary targets for network pharmacology analysis. We obtained the sequences of S1R (Q99720), DRD2 (P14416), and BIP (P11021) from the UniProt repository (Consortium, 2019). In addition, we obtained S1R (PDB ID: 5HK1) (Schmidt et al., 2016), DRD2 (6 PDB ID: LUQ) (Fan et al., 2020), and BIP (PDB ID: 3LDN) (Macias et al., 2011) from the RCSB Protein Data Bank (PDB) (Berman et al., 2000), which are 2.51 Å, 3.10 Å, and 2.20 Å, respectively, and their structures are shown in Figure 2A.

The half-inhibitory concentration (IC₅₀) refers to the concentration of the drug or inhibitor required to inhibit half of the specified biological process, and the inhibition constant Ki reflects the inhibitory strength of the inhibitor on the target. The smaller the value, the stronger the inhibitory ability. pIC₅₀ is the negative logarithm of the IC₅₀ value, which is usually used to characterize the activity of molecules in drug screening. The formula for converting IC₅₀ values to pIC₅₀ values is p I C 50 = − log 10 I C 50 . (1)

We obtained data on inhibitors of S1R, DRD2, and BIP and their binding abilities to their targets from the ChEMBL database (Gaulton et al., 2012). Although both IC₅₀ and Ki can reflect the activity of the inhibitor, for data consistency, we screened the inhibitor data with IC₅₀ as the subsequent drug–target affinity (DTA) training data on the HNN. Similarly, under the premise of ensuring the number of datasets, we screened the data whose source description was scientific literature and excluded other data. Figure 2B shows the simplified molecular input line entry system (SMILES) length distribution and binding force distribution of the three protein inhibitors. The inhibitor distribution of S1R and DRD2 showed a Gaussian distribution trend, while the inhibitor distribution of BIP was relatively sparse.

The drug screening library used in this study comes from FDA-approved drugs in the DrugBank database (Wishart et al., 2008). DrugBank is a comprehensive pharmaceutical knowledge bank that provides pharmacists, pharmacologists, health professionals, and drug researchers with free academic resources to help advance drug development and clinical practice. We chose DrugBank as the screening bank because it contains extensive drug information and a list of FDA-approved drugs, which can be used to screen potential drugs for the treatment of AD. These drugs have been proven to be safe and effective treatments in human clinical trials, so they are expected to be used in the treatment of AD. We selected FDA-approved drugs in the DrugBank database as screening libraries, and a total of 2509 drug molecules were available for drug repurposing studies.

The overview of the HNNDTA framework proposed in this study is shown in Figure 3. First, we used a network pharmacology approach to find other targets in the same pathway as the AD target S1R, namely, DRD2 and BIP. We searched the ChEMBL website for inhibitor data for these three targets. The target protein is encoded as a one-dimensional target embedding, and the drug molecule is encoded as a one-dimensional drug embedding. The two encoding vectors are spliced in zero dimension, and after the calculation of the deep neural network (DNN), the final DTA is obtained, which can be expressed as follows: D T A = D N N c a t v p , v d , (2) where the function c a t a , b represents the splicing operation of the 1D a and b vectors, and v _p and v _d represent the encoding vectors of the target protein and the drug molecule, respectively. In this paper, there are 13 kinds of target encoders and 9 kinds of drug encoders, all of which are built by DeepPurpose (Huang et al., 2020). A suitable combined model will produce better prediction accuracy. During the training phase, the dataset was randomly divided into independent training, validation, and test sets in a ratio of 7:1:2. The training set was used to train the model, while the validation and test sets were used to evaluate its performance. Due to the nature of our HNNDTA framework, which was trained on datasets specific to individual targets, it exhibits higher predictive accuracy for single targets. We have observed that models trained on single targets exhibit higher accuracy than those trained on mixed-target datasets.

The drug encoder receives SMILES sequences as input. The Morgan encoder first uses the ECFP (Rogers and Hahn, 2010) algorithm to generate the feature representation sequence of the circular substructure of the drug, with a length of 1,024 bits. A multi-layer perceptron (MLP) then processes the sequence of feature representations to obtain a vector representation that can be fed into a neural network. The Morgan encoder is expressed as follows: f morg S M I L E S = M L P E C F P S M I L E S . (3)

Similar to the Morgan encoder, the daylight encoder also uses the ECFP algorithm to generate a feature sequence based on the channel substructure of the drug, which is used as the input of the multi-layer perceptron to generate a feature sequence with a length of 2048 bits. The daylight encoder is represented as follows: f day S M I L E S = M L P E C F P S M I L E S . (4)

The PubChem encoder (Kim et al., 2019) generates feature sequences using handcrafted important substructures and then generates a feature sequence with a length of 881 bits through a multi-layer perceptron. The PubChem encoder is represented as follows: f pub S M I L E S = M L P S u b s t r u c t u r e S M I L E S . (5)

The rdkit_2d_normalize encoder (Reczko and Bohr, 1994) generates a feature sequence with a length of 200 bits according to the global pharmacophore of the drug and then normalizes the feature sequence by fitting the cumulative density function of a given molecule sample. The rdkit_2d_normalize encoder is represented as follows: f rdkit S M I L E S = M L P N o r m a l i z e F e a t u r e F e a t u r e = G l o b a l P h a r m a c o p h o r e S M I L E S . (6)

The extended reeb graph (ErG) method (Stiefl et al., 2006) mixes the simplified graph and the binding attribute pair to generate a feature sequence and uses the node description of the drug carrier type to encode the relevant molecular properties; the encoded features are obtained after the MLP calculation vector. The ErG coder is expressed as follows: f erg S M I L E S = M L P G r a p h F e a t u r e F e a t u r e = S c a f f o l d − B a s e d N o d e D e s c r i p t o r S M I L E S . (7)

MLP obtains the output value through feedforward propagation and updates the model parameters through reverse transmission so that the model output value gradually approaches the real value. The output of the MLP forward propagation is expressed as follows: y = A C ∑ i = 1 M l ω i , l • A C ∑ i = 1 M l − 1 ω i , l − 1 • … , (8) where AC is the activation function and the typical activation function is the modified linear unit ReLU; M _l is the number of neurons in the lth layer network, ω _i,l is the weight of the ith neuron in the lth layer network; and the termination condition of … in the aforementioned formula is the first layer of the neural network, that is, the input layer. The reverse transfer uses the Adam optimizer to update the model weights. The underlying algorithm is the gradient descent method. The update on the weight of ω _i,l can be expressed as follows: ω i , l ← ω i , l − η ∂ E ∂ ω i , l , (9) where E is the difference between the predicted value and the real value, ∂ E ∂ ω i , l is the partial derivative of E to ω _i,l , and η is the learning rate.

The input to a target encoder is the amino acid sequence of the target. The signature sequence generated by the amino acid composition (AAC) coder is 8420 positions in length, where each position is consistent with the maximum length of overlapping subsequences (k-mers) of one amino acid. The amino acid composition coder is expressed as follows: A A C i = f i L , i = 1,2 , … , 20 , (10) where f _i represents the number of occurrences of amino acid i in the protein and L represents the length of the amino acid sequence. The AAC encoder concatenates 20 AAC values for each position in the amino acid sequence to obtain a signature sequence of 8420 elements in length.

The pseudo amino acid composition (PseAAC) encoder adds the hydrophobic and hydrophilic pattern information on the protein based on AAC to generate a 30-bit feature vector representation. The pseudo-amino acid composition encoder is expressed as follows: P s e A A C i , j = ∑ k = 1 L f k , i w k , j ∑ k = 1 L f k , i , i = 1,2 , … , 20 ; j = 1,2 , … , 30 , (11) where f _k,i represents the frequency of amino acid i in the kth position in the protein sequence and w _k,j represents the weight of the pattern of the kth amino acid and the relative position j. The PseAAC encoder concatenates 30 PseAAC values for each position in the amino acid sequence, resulting in a feature vector of 30 elements in length.

The conjoint triad (ConTriad) encoder (Shen et al., 2007) forms a 7-letter alphabet based on amino acid triplet features, generating a feature vector with a length of 343 elements. The ConTriad encoder is expressed as follows: C o n T r i a d i = ∑ j = 1 7 f j , i w j L − 2 , i = 1,2 , … , 3430 , (12) where f _j,i indicates that the three adjacent amino acids in the protein sequence are converted into a number according to the 7-letter alphabet, the ith element indicates the frequency of the jth triplet appearing in the protein sequence, and w _j is the weight of the jth triplet. The ConTriad encoder concatenates 343 ConTriad values for each position in the amino acid sequence, resulting in a feature vector of 343 elements in length.

The quasi-sequential encoder consists of a 100-element feature vector of quasi-sequential descriptors (Chou, 2000). The feature vectors generated by the aforementioned manual feature encoder will be further processed as input to MLP to obtain the feature vector of the target. The quasi-sequential encoder is expressed as follows: Q u a s i S e q i = ∑ j = 1 N ρ j d i j , i = 1 , 2 , … , 100 , (13) where ρ _j represents the weight of the jth quasi-sequential descriptor and d _ij is the distance between the ith amino acid and the jth sequence descriptor. The quasi-sequential encoder concatenates 100 QuasiSeq values for each position in the amino acid sequence, resulting in a feature vector of 100 elements in length.

The aforementioned drug and target feature extraction methods are based on prior chemical knowledge and manual transformation, so these encoders cannot be mixed. The encoders introduced in this section are general-purpose encoders based on DNNs, including convolutional neural networks (CNNs) (Krizhevsky et al., 2017), gated recurrent units (GRUs) (Chung et al., 2014), long short-term memory (LSTM) (Hochreiter and Schmidhuber, 1997), and transformers (Vaswani et al., 2017). These neural networks treat amino acid sequences as one-dimensional data.

The CNN encoder is a multilayer 1D CNN (Krizhevsky et al., 2017). After encoding the amino acid sequence character by character, the obtained deep feature vector will pass through multiple 1D convolutional layers and finally pass through the one-dimensional maximum pooling layer to obtain the output of the target feature vector. The output of the 1D convolutional layer is the result of convolving the input with the convolution kernel, which can be expressed as follows: o u t = i n p u t ⊗ k e r n e l , (14) where ⊗ represents a convolution operation. Assuming that the convolution kernel size is 2k + 1, k ∈ N ⁺, the ith convolution output can be expressed as follows: o u t i = ∑ j = i − k i + k ∑ a = 1 2 k + 1 i n p u t j • k e r n e l a . (15)

GRU and LSTM encoders are types of recurrent neural networks. In both networks, each node will get an output based on the state at the last moment and the current input and update the state of the node. This can solve the problem of traditional convolutional networks without long-term memory to a certain extent. Specifically, the SMILES sequence or amino acid sequence will first pass through the CNN for feature extraction and then use the output of the CNN as the input of the recurrent network.

The transformer encoder applies a self-attention mechanism (Vaswani et al., 2017). Due to the computational time and memory cost of the transformer, amino acid sequences are decomposed into moderately sized protein substructures, such as motifs, and each segmentation is then treated as a token and fed into a self-attention-based encoder. If a SMILES sequence or amino acid sequence is treated as a sentence, cut into several meaningful phrases, and encoded into several vectors with the same number of phrases, denoted as x, then the output of the transformer can be expressed as x 1 = norm x + a t t n x , m a s k , o u t = norm x 1 + f e e d f o r w a r d x 1 , (16) where attn is a self-attention function, mask is a Boolean value about whether the input x is eliminated, feedforward is a feedforward neural network, and norm is a layer normalization operation.

In mathematical statistics, mean-squared error (MSE) is a method used to measure the difference between the predicted and real values. It calculates the mean of the squared difference between predicted and true values, which is the expected value of the squared difference between predicted and true values. The smaller the value of the MSE, the higher the prediction accuracy of the prediction model. Assuming there are n samples, MSE can be expressed by the following formula: M S E = 1 n ∑ i = 1 n y i − y ̂ i 2 , (17) where y _i and y ̂ i are the true and estimated values of the ith sample, respectively. In this paper, the MSE is used to evaluate the accuracy of the model to predict the binding affinity of the drug to the target.

Harrell’s C-index (also known as the concordance index, CI) is a widely used metric for evaluating the performance of risk models. It is commonly employed in survival analysis, especially when dealing with censored data (Harrell et al., 1982). The C-index measures the degree of concordance between predicted and observed rankings of survival times. It serves as an indicator of the model’s accuracy with values closer to 1, indicating a higher level of consistency between the predicted outcomes and the actual observed outcomes.

Suppose the data are represented by vectors ( T i ̃ , Δ i , X i 1 , … , X i p ) , i = 1 , … , n , where T i ̃ is a possibly right-censored continuous survival time and ( X i 1 , … , X i p ) T is a vector of predictor variables. It is assumed that T i ̃ is the minimum of the true survival time T _i and an independent continuous censoring time C _i . The variable Δ _i ≔I(T _i < = C _i ) indicates whether T _i has been fully observed (Δ _i = 1) or not (Δ _i = 0). A one-dimensional score η i ∈ R is estimated for each observation i = 1, …, n, by averaging the cumulative hazard estimates over all trees and all time points. The concordance index is given by C = ∑ i , j I T i ̃ > T j ̃ ⋅ I η j > η i ⋅ Δ j ∑ i , j I T i ̃ > T j ̃ ⋅ Δ j , (18) where the indices i and j refer to pairs of observations in the sample (Schmid et al., 2016).