11institutetext: Jamia Millia Islamia University, New Delhi, India
11email<EMAIL_ADDRESS>
11email<EMAIL_ADDRESS>
22institutetext: The Chinese University of Hong Kong, Hong Kong SAR
22email<EMAIL_ADDRESS>
# A Methodology for the Prediction of Drug Target Interaction using CDK
Descriptors

Tanya Liyaqat 11 Tanvir Ahmad 11 Chandni Saxena 22

###### Abstract

Detecting probable Drug Target Interaction (DTI) is a critical task in drug
discovery. Conventional DTI studies are expensive, labor-intensive, and take a
lot of time, hence there are significant reasons to construct useful
computational techniques that may successfully anticipate possible DTIs.
Although certain methods have been developed for this cause, numerous
interactions are yet to be discovered, and prediction accuracy is still low.
To meet these challenges, we propose a DTI prediction model built on molecular
structure of drugs and sequence of target proteins. In the proposed model, we
use Simplified Molecular-Input Line-Entry System (SMILES) to create CDK
descriptors, Molecular ACCess System (MACCS) fingerprints, Electrotopological
state (Estate) fingerprints and amino-acid sequences of targets to get Pseudo
Amino Acid Composition (PseAAC). We target to evaluate performance of DTI
prediction models using CDK descriptors. For comparison, we use benchmark data
and evaluate models’ performance on two widely used fingerprints, MACCS
fingerprints and Estate fingerprints. The evaluation of performances shows
that CDK descriptors are superior at predicting DTIs. The proposed method also
outperforms other previously published techniques significantly.

###### Keywords:

Drug Target InteractionsCatBoostCDK descriptorsMolecular fingerprints

##  1 Introduction

Drug target interaction (DTI) is a prominent task in drug discovery and
research. It entails detecting possible links among chemical compounds and
protein targets which acts as a guide in the preliminary phases of drug
discovery and developmental research. Experiments carried out in wet labs are
labor intensive and require a significant amount of money [22]. According to
statistics, each novel molecular entity takes around 1.8 billion USD and the
authorization of a novel drug application usually requires at least 9 years
[8]. As a result, high-efficiency computational prediction techniques to
investigate drug target interactions based on Machine Learning (ML) and Deep
Learning (DL) have sparked a lot of attention in recent years [2]. The bonding
of a medicine to a target’s location resulting in the alteration of its
functioning is considered as drug target interaction. Any chemical molecule
that causes an alteration in the body’s physiology when swallowed, ingested,
or inhaled is referred to as a medication or medicine. On the other hand,
targets consist of elements as nucleic acids or lipids, that are intended to
modify. Ion channels, enzymes, nuclear receptors, and G-protein coupled
receptors are among the most popular biological targets. To treat illness and
ailments, the medicine inhibits the target’s function in order to prevent
certain catalytic processes from occurring in the human body. This is
accomplished by preventing it from interacting with particular enzymes known
as substrates. The drug discovery procedure that detects novel therapeutic
molecules for targets relies heavily on DTI prediction [24]. Feature-based
computational techniques for DTI prediction have gained significant attention
over the years. The availability of the structural information of chemical
compounds in the form of fingerprints or descriptors has played an important
role. However, most studies consider fingerprints over descriptors. Hence, it
becomes important to compare performance and identify better alternative. We
provide more details about feature-based techniques in Section 2.

Considering the widely accepted ability of structure information of molecules,
we aim to evaluate the performance of CDK descriptors against two widely used
fingerprints, namely Molecular ACCess System (MACCS) and Electrotopological
state (Estate) fingerprints. The proposed model utilizes Pseudo amino acid
composition derived using amino-acid sequences of targets via _iFeature
webserver_ [4]. We use drug Simplified Molecular-Input Line-Entry System
(SMILES) to obtain CDK descriptors, MACCS fingerprints and Estate
fingerprints. The purpose here is to evaluate the impact of employing CDK
descriptors for DTI prediction. We compare models’ performance against two
frequently used fingerprints, MACCS fingerprints and Estate fingerprints on
benchmark data. This work mainly focus on extracting and feature processing,
followed by a systematic prediction methodology based on machine learning. For
example, in this case, we utilize the Categorical Boosting (CatBoost)
classifier to make predictions. For validation, we compare our proposed model
to several recently proposed models. The results reveal that the proposed DTI
prediction model identifies drug-target interactions more accurately using CDK
descriptors than MACCS and Estate fingerprints.
We organize the paper as follows. Section 2 offers an overview of computation
approaches to DTI predictions and highlights recent methods closely related to
our work. Section 3 provides the details about the datasets and feature
encodings. Section 4 describes our proposed methodology and a brief overview
of the CatBoost algorithm. Evaluation metrics and performance results are
presented in Section 5 and Section 6 respectively. Finally, we conclude the
paper in Section 7.

##  2 Computation Approaches to DTI Predictions

In this section, we provide an overview of computation approaches to DTI
predictions and highlight some closely related work to our proposed
methodology. The computational strategies for the prediction of DTIs can be
broadly divided into ligand-based, docking-based, and chemogenomic approaches
[11].

Ligand-based. The rationale behind ligand-based techniques is that identical
compounds bind to identical biological targets and have identical features. It
starts with a single molecule or a group of chemicals known to be effective
against the target and it is further guided by the structure-activity
relationships. However, there are certain drawbacks to this strategy. Because
the protein sequence information is not employed in the prediction process,
discovering new interactions reduces the connection across the identified
ligand and protein families [9].

Docking-based. The docking-based approach, on the other hand, uses the 3D
shape of proteins and chemical compounds to determine their possibilities of
interaction [14, 6]. However, specific proteins like the membrane proteins
have unknown 3D structures that make it less applicable [20].

Figure 1: A brief taxonomy of computational approaches to DTI prediction

Chemogenomic. The chemogenomic approaches use drug and protein information
together to anticipate interactions. To infer probable interactions, a shared
subspace is created by unifying the biochemical space of drugs and the genome
space of targets. The main benefit of this method is that it utilizes a
significant amount of biological data that is freely accessible from public
repositories [33]. Chemogenomic approaches are roughly divided into network-
based methods and feature-based methods. Network-based approaches integrate
data like drug-drug interactions, protein-protein interactions, drug-disease
interactions, and drug-target interactions from multiple sources into a single
unified framework to boost DTI prediction [3, 12, 16, 17, 35]. For instance,
Wan et al. [28] devised an end-to-end technique entitled _NeoDTI_ to combine
data from omics networks and learn topology that preserves the information of
drugs and targets. Recent years have seen a fast growth of ML models based on
knowledge graphs (KG). Mohammad et al. [19] suggest _triModel_ , a model based
on Knowledge Graph (KG) embeddings to derive novel DTI from the model’s scores
built by learning embeddings about drugs and targets from multi-modal
heterogeneous data. On the other hand, feature-based approaches represent each
drug target pair as an array of descriptors. Drugs and proteins are
transformed into corresponding descriptors based on their chemical properties.
Integration of individual features of drugs and targets forms the input to
these approaches as a 1D array [15, 21, 34].

Most researchers prefer feature-dependent computation models to predict DTIs
focussed on structural information of drugs based on molecular fingerprints
that are bit strings indicating the existence of a specific substructure. For
example, Han et al. [27] present an automatic learning system called _LRF-
DTIS_ retrieving drug and target characteristics in the form of _PsePSSM_ and
molecular fingerprints employing lasso for feature selection, smote for
handling imbalance, and random forest (RF) for predicting interactions. Wenyu
et al. [18] suggest _PDTI-ESSB_ turning all drug molecules into molecular
substructure fingerprints and representing protein sequences as multiple
features to express their evolutionary, sequential, and structural
information. To prevent the drawbacks of sparseness and dimensionality curse,
Wang et.al. [30] propose _MSPEDTI_ which uses a Convolutional Neural Network
to derive relevant low-dimensional features from the sequence and structural
information defined in the form of PSSM and molecular fingerprints. Likewise,
Wang et al. [31] propose a similar method using _PsePSSM_ and _PubChem_
fingerprints with feature weighted Rotation Forest. Sajadi et al. [25]
introduce an interesting approach to handling the sparsity in the interaction
matrix of drugs and targets through drug fingerprints. Wang et al. [29] use
_MACCS_ and _PAAC_ for predicting DTIs with a novel method to create negative
samples. Another way of capturing the structural information is molecular
descriptors which are theoretically derived properties representing the
physical, chemical, and topological characteristics of drugs. The majority of
researchers use benchmark data to assess the effectiveness of their methods.
We utilize the suggested benchmark data to evaluate our proposed method. In
Section 3, we give a brief description about the data and feature
representation for drugs and proteins.

##  3 Materials and Methods

###  3.1 Dataset Description

The gold-standard datasets used in this investigation are enzyme, Ion Channel
(IC), Nuclear Receptor (NR) and G-Protein Coupled Receptor (GPCR), have been
compiled by Yamanishi et al. [33] from the SuperTarget [10], BRENDA [26],
DrugBank [32] and KEGG [13] repositories. The number of DTI pairs in these
datasets after deleting the redundant information are $2926$, $635$, $1476$,
and $90$, respectively. We consider each of these combinations as a positive
interaction. To ensure a balance between positive and negative interactions,
random under-sampling is applied to all pairs of DTIs to generate negative
interactions. Table 1 displays the statistical data regarding these benchmark
datasets.

Table 1: Qualitative details of the datasets. Dataset |  Drugs |  Targets |  Known Interaction |  Ratio of Imbalance
---|---|---|---|---
Ion channel |  201 |  204 |  1476 |  0.036
Enzyme |  445 |  664 |  2926 |  0.010
Nuclear Receptor |  54 |  26 |  90 |  0.068
GPCR |  223 |  95 |  635 |  0.031

###  3.2 Drug Feature Representation

In this study, molecular descriptors and molecular fingerprints convert the
molecular structures of drugs into numerical form. For molecular descriptors,
we use CDK descriptors derived from
ChemDes-1-1-1http://www.scbdd.com/chemdes/, for molecular fingerprints, namely
MACCS and Estate, we use the RDKit
library000https://www.rdkit.org/docs/GettingStartedInPython.html. The concept
behind fingerprints is to define molecular structure using a library of
molecular substructures, that translates a chemical compound into a bit vector
of 0’s and 1s by detecting whether a compound has a particular substructure or
not. If the compound has that substructure, then only the bit at that position
is put to 1 which results in a characterization of the molecular structure as
a binary string. The number of molecular substructures in MACCS and Estate is
166 and 79 respectively, which is the size of the final binary bit vector. CDK
descriptors, on the other hand, is an open-source platform for detecting and
categorizing compounds using descriptor classes. The final size of the CDK
descriptors is 275 and consists of autocorrelation descriptors, connectivity
descriptors, constitutional descriptors, kappa descriptors, molecular
properties, topological descriptors, WHIM descriptors, CPSA descriptors,
geometrical descriptors, and quantum chemical descriptors.

###  3.3 Protein Feature Representation

PseAAC is a parallel-correlation-based method for delineating protein
information that generates 20 + D features. There are various models based on
amino acid compositions (AAC) that lack target sequence-order knowledge. Pse-
AAC, introduced in [7] can be used to express both AAC and AA sequence order
data. This technique is common in bioinformatics and related areas. Pse-AAC
combines the core features of AAC with certain additional parts indicating a
set of protein correlation factors and helps to improve model performance for
multiple tasks. The expression for the features of PseAAC is as shown in Eq.
1.

$\centering
P=[p_{1},p_{2},p_{3},...p_{19},p_{20},...,p_{20+\lambda}],\quad(\lambda<S)\
,\@add@centering$ (1)

where S corresponds to the size of the target sequence under consideration.
The components are represented as shown below:

$P_{j}=\frac{F_{j}}{\sum_{j=1}^{20}F_{j}+W\sum_{k=1}^{\lambda}\psi_{j}},\quad
1<<20\ ,\\\ $ (2)
$P_{j}=\frac{W\psi_{j}}{\sum_{j=1}^{20}F_{j}+W\sum_{k=1}^{\lambda}\psi_{j}},\quad
20+1<<20+\lambda\ .$ (3)

Here, P stands for a vector of features, and W stands for a weight factor of
0.05. $F_{j}$ shows amino acid occurrence frequency, $\psi$ is the sequence
correlation factor and $\lambda$ represents information about the protein
sequence order. In Eq. 2, for a single protein, AAC is expressed by 20
components, whereas sequence order is expressed by $20+1$ to $20+\lambda$
components, referred to as Pse-AAC shown in Eq. 3. Finally, for each protein
sequence, we obtained a 50-D feature vector by setting the value lambda to 30.
Rcpi [1], iFeature [4], and iLearn [5] are among the useful software tools for
encoding protein sequences. Section 4 defines the work flow of DTI prediction
and shows the use of PAAC in our proposed methodology.

Figure 2: Proposed Methodology

##  4 Proposed Methodology

The proposed workflow of our DTI prediction model is as shown in Fig. 2. It is
comprised of three main steps. The description of each step is as follows:

  * •

Step 1. At first, we transform data into the corresponding values for which we
represent each drug and target using their SMILES and amino acid sequences,
respectively. Further, we use them to encode drug features and target
features. For drugs, we use CDK descriptors, MACCS fingerprints and Estate
fingerprints as features. For target protein, we use encoded features in the
form of Pseudo amino acid composition (PseAAC).

  * •

Step 2. In the next step, we process features from the input data. For each
DTI pair, the drug features are concatenated with the target features to
produce the final 1D array of features. The size of the array varies according
to the drug feature under consideration. For instance, in the case of CDK
descriptors, the length turns out to be $325$. For MACCS and Estate
fingerprints, it is $216$ and $129$, respectively. To normalize the range of
independent attributes, we use the StandardScaler() function that converts the
training and testing data to scaled representation. It accelerates
calculations in machine learning techniques.

  * •

Step 3. In the final step, we provide the processed features and existing DTIs
to a classification model as input. We use the CatBoost algorithm to train a
classifier. After this stage, we measured the performance of the trained model
on a test dataset. We provide further details about the CatBoost algorithm in
the following subsection.

###  4.1 CatBoost

Gradient boosted trees and Random forest is robust machine models for tabular
heterogeneous data. CatBoost classifier is another open-source gradient
boosting framework released in 2017 [23]. Although mainly designed to handle
category features, CatBoost works on numerical and text data. According to the
literature, it outperforms boosting algorithms such as XGBoost and LightGBM on
a variety of datasets and has a substantially shorter prediction time. The
technique is well-known for its use of ordered boosting to counteract
overfitting and the use of symmetric trees for faster execution. During model
training, a sequence of decision trees is built one after the other, with each
succeeding tree having a smaller loss. In other words, each decision tree
learns from its predecessor and influences the subsequent trees to increase
its performance, culminating in a powerful learner. Gradient boosting trees
are good at dealing with numerical data but struggle with categorical
features. Both strategies require a large amount of memory and are
computationally expensive. CatBoost, as a solution, uses target-based
statistics to address categorical features, saving time and resources. To
overcome overfitting, the CatBoost method employs an ordered boosting
mechanism. After numerous boosting stages, traditional gradient boosting
approaches use all the training samples to develop a prediction model. This
strategy causes a prediction shift in the created model, resulting in a
distinct form of target leakage problem. CatBoost eliminates the
aforementioned challenge by utilizing an ordered boosting architecture.
Furthermore, unlike standard learning classifiers, the CatBoost approach
gracefully handles overfitting by employing many permutations of the training
dataset, which emerges as the main reason for deploying its intelligence in
the current work. Five-fold cross-validation is used to avoid data overfitting
and for the appropriate application of CatBoost in the current DTI problem. We
evaluate the performance of the created model on multiple performance metrics
defined in the next section.

##  5 Evaluation Parameters

To evaluate and compare our proposed methodology with other methods, we use
different metrics such as Precision, Sensitivity (or Recall), Accuracy,
Mathews Correlation Coefficient (MCC), Area Under Curve (AUC), and Area Under
Precision-Recall( AUPR). The AUC and AUPR graphs are good choices for
unbalanced data. Hence, most research uses it as a comparison criterion. The
ROC curve demonstrates how well the trained model performs at different
cutoffs. False positive rates are compared with actual true positive rates to
form the curve. The AUC values range from 0 to 1, with higher values
suggesting an effective model. As the AUC summarizes the curve with a range of
cutoff values as a single score, the AUPR also does the same. The difference
is that it shows the precision (y-axis) and recall (x-axis) for various
probability cutoffs. The following are the definitions of other evaluation
parameters used in this study:

$\displaystyle Precision$ $\displaystyle=\frac{TP}{TP+FP}\ ,$ (4)
$\displaystyle Sensitivity$ $\displaystyle=\frac{TP}{TP+FP}\ ,$ (5)
$\displaystyle Accuracy$ $\displaystyle=\frac{TP+TN}{TP+FP+TN+FN}\ ,$ (6)
$\displaystyle MCC$ $\displaystyle=\frac{TP*TN-
FP*FN}{(TP+FP)(TN+FN)(TN+FP)(TN+FN)}\ ,$ (7)

where TP means correctly labeled positive observations, FP implies wrongly
labeled negative observations, TN means correctly labeled negative
observations, and FN implies wrongly labeled positive observations.

##  6 Performance Evaluation

To evaluate the prediction capacity of our proposed approach, we use 5-fold
cross-validation on all DTI networks, including IC, Enzymes, NR, and GPCR. We
roughly divide the data into 5-folds and use four folds for the model training
and the left-out fold for testing the model’s performance. We repeated the
process five times to get five different result sets and reported the mean
values. Table 2 displays the experimental findings of five-fold cross-
validation on various performance evaluators across all datasets.

Table 2: Obtained results using CDK descriptors, MACCS fingerprints, and
Estate fingerprints across all datasets for different performance metrics.
Dataset |  Feature |  Accuracy |  Precision |  Recall |  MCC |  AUC |  AUPR
---|---|---|---|---|---|---|---
Enzyme |  CDK |  0.8988 |  0.9015 |  0.8954 |  0.7977 |  0.9533 |  0.9597
MACCS |  0.8975 |  0.9013 |  0.89262 |  0.7915 |  0.9548 |  0.9595
Estate |  0.8822 |  0.8854 |  0.8783 |  0.7646 |  0.9468 |  0.9521
GPCR |  CDK |  0.8409 |  0.8304 |  0.8567 |  0.6826 |  0.8995 |  0.8967
MACCS |  0.8314 |  0.8178 |  0.8551 |  0.6650 |  0.9063 |  0.9029
Estate |  0.8196 |  0.8104 |  0.8346 |  0.6400 |  0.8818 |  0.8713
IC |  CDK |  0.8929 |  0.8864 |  0.8997 |  0.7848 |  0.9539 |  0.9591
MACCS |  0.8872 |  0.8793 |  0.8977 |  0.7750 |  0.9464 |  0.94767
EState |  0.8523 |  0.8409 |  0.8692 |  0.7058 |  0.9235 |  0.9174
NR |  CDK |  0.8056 |  0.7788 |  0.8556 |  0.6155 |  0.8445 |  0.8376
MACCS |  0.75 |  0.7560 |  0.7778 |  0.5101 |  0.8269 |  0.8342
EState |  0.7611 |  0.7837 |  0.7333 |  0.5332 |  0.8429 |  0.8327

Table 2 shows that we obtain the best accuracy of 0.8988 and the second best
accuracy of 0.8929 using CDK descriptors on the enzyme and ion channel
dataset, respectively. Results also show that the performance of CDK
descriptors is better in most cases compared to MACCS and Estate fingerprints,
which confirms that CDK descriptors are better at predicting the DTIs compared
to these two widely used molecular fingerprints. Fig. 3 shows the performance
curves of the area under the ROC curve across all benchmark datasets.

(a)

(b)

(c)

(d)

Figure 3: Comparative analysis of the area under the ROC curves obtained using
CDK descriptors, MACCS Fingerprints and Estate fingerprints across all
datasets.

###  6.1 Comparison with previous methods

To evaluate the proposed model’s capacity for predicting DTIs in a logical
manner, we contrasted it with earlier techniques using the gold-standard
dataset and chose the AUC as the assessment measure since it best captures the
model’s performance. Table 3 aggregates the AUC values from prior approaches
such as Yamanishi [34], DTCWT [21], Yang [15], Elastic net [27], RoFDT [31],
AutoDTI++ [25] and MSPEDTI [30]. Table 3 clearly shows that the proposed model
outperformed the prior technique across benchmark datasets. This shows that
using CDK descriptors with the CatBoost classifier can significantly improve
the capacity to anticipate DTIs.

Table 3: Comparison of the proposed methodology with previous studies with respect to AUC values. Methods |  Enzyme |  GPCR |  IC |  NR
---|---|---|---|---
Yamanishi [34] |  0.821 |  0.811 |  0.692 |  0.814
Elastic net [27] |  0.8605 |  0.7785 |  0.804 |  0.8418
Yang [15] |  0.9529 |  0.8878 |  0.925 |  0.8487
DTCWT [21] |  0.9498 |  0.8775 |  0.9270 |  0.7755
RoFDT [31] |  0.9172 |  0.8557 |  0.8827 |  0.7531
AutoDTI++ [25] |  0.90 |  0.86 |  0.91 |  0.87
MSPEDTI [30] |  0.9437 |  0.8802 |  0.9088 |  0.8663
Proposed Model |  0.9533 |  0.8995 |  0.9539 |  0.8445

##  7 Conclusion

We presented a machine learning-based prediction model for DTIs. We use random
under-sampling to deal with the imbalance of the datasets. To encode features
such as CDK descriptors, MACCS fingerprints, and Estate fingerprints for drugs
and PseAAC for targets, we use drug chemical structures and amino acid
sequences. The objective of this research is to evaluate the impact of
employing CDK descriptors for DTI prediction. We compare its performance
against two frequently used fingerprints, namely MACCS fingerprints, and
Estate fingerprints. The experimental findings reveal that CDK descriptors
outperform the other two commonly used fingerprints. We use five-fold cross-
validation criteria to get the results. The provided methodology is both
practical and effective in forecasting DTIs, according to the comparative
outcomes. We intend to expand our research in the future by taking into
account novel target feature descriptors coupled with deep learning
techniques.

## References

  * [1] Cao, D.S., Xiao, N., Xu, Q.S., Chen, A.F.: Rcpi: R/bioconductor package to generate various descriptors of proteins, compounds and their interactions. Bioinformatics 31(2), 279–281 (2015)
  * [2] Chen, R., Liu, X., Jin, S., Lin, J., Liu, J.: Machine learning for drug-target interaction prediction. Molecules 23(9), 2208 (2018)
  * [3] Chen, Z.H., You, Z.H., Guo, Z.H., Yi, H.C., Luo, G.X., Wang, Y.B.: Prediction of drug–target interactions from multi-molecular network based on deep walk embedding model. Frontiers in Bioengineering and Biotechnology 8, 338 (2020)
  * [4] Chen, Z., Zhao, P., Li, F., Leier, A., Marquez-Lago, T.T., Wang, Y., Webb, G.I., Smith, A.I., Daly, R.J., Chou, K.C., et al.: ifeature: a python package and web server for features extraction and selection from protein and peptide sequences. Bioinformatics 34(14), 2499–2502 (2018)
  * [5] Chen, Z., Zhao, P., Li, F., Marquez-Lago, T.T., Leier, A., Revote, J., Zhu, Y., Powell, D.R., Akutsu, T., Webb, G.I., et al.: ilearn: an integrated platform and meta-learner for feature engineering, machine-learning analysis and modeling of dna, rna and protein sequence data. Briefings in bioinformatics 21(3), 1047–1057 (2020)
  * [6] Cheng, A.C., Coleman, R.G., Smyth, K.T., Cao, Q., Soulard, P., Caffrey, D.R., Salzberg, A.C., Huang, E.S.: Structure-based maximal affinity model predicts small-molecule druggability. Nature biotechnology 25(1), 71–75 (2007)
  * [7] Chou, K.C.: Prediction of protein cellular attributes using pseudo-amino acid composition. Proteins: Structure, Function, and Bioinformatics 43(3), 246–255 (2001)
  * [8] Dickson, M., Gagnon, J.P.: Key factors in the rising cost of new drug discovery and development. Nature reviews Drug discovery 3(5), 417–429 (2004)
  * [9] Ezzat, A., Wu, M., Li, X.L., Kwoh, C.K.: Computational prediction of drug–target interactions using chemogenomic approaches: an empirical survey. Briefings in bioinformatics 20(4), 1337–1357 (2019)
  * [10] Günther, S., Kuhn, M., Dunkel, M., Campillos, M., Senger, C., Petsalaki, E., Ahmed, J., Urdiales, E.G., Gewiess, A., Jensen, L.J., et al.: Supertarget and matador: resources for exploring drug-target relationships. Nucleic acids research 36(suppl_1), D919–D922 (2007)
  * [11] Jacob, L., Vert, J.P.: Protein-ligand interaction prediction: an improved chemogenomics approach. bioinformatics 24(19), 2149–2156 (2008)
  * [12] Ji, B.Y., You, Z.H., Jiang, H.J., Guo, Z.H., Zheng, K.: Prediction of drug-target interactions from multi-molecular network based on line network representation method. Journal of Translational Medicine 18(1), 1–11 (2020)
  * [13] Kanehisa, M., Goto, S., Hattori, M., Aoki-Kinoshita, K.F., Itoh, M., Kawashima, S., Katayama, T., Araki, M., Hirakawa, M.: From genomics to chemical genomics: new developments in kegg. Nucleic acids research 34(suppl_1), D354–D357 (2006)
  * [14] Li, H., Gao, Z., Kang, L., Zhang, H., Yang, K., Yu, K., Luo, X., Zhu, W., Chen, K., Shen, J., et al.: Tarfisdock: a web server for identifying drug targets with docking approach. Nucleic acids research 34(suppl_2), W219–W224 (2006)
  * [15] Li, Y., Liu, X.z., You, Z.H., Li, L.P., Guo, J.X., Wang, Z.: A computational approach for predicting drug–target interactions from protein sequence and drug substructure fingerprint information. International Journal of Intelligent Systems 36(1), 593–609 (2021)
  * [16] Liu, Y., Wu, M., Miao, C., Zhao, P., Li, X.L.: Neighborhood regularized logistic matrix factorization for drug-target interaction prediction. PLoS computational biology 12(2), e1004760 (2016)
  * [17] Luo, Y., Zhao, X., Zhou, J., Yang, J., Zhang, Y., Kuang, W., Peng, J., Chen, L., Zeng, J.: A network integration approach for drug-target interaction prediction and computational drug repositioning from heterogeneous information. Nature communications 8(1), 1–13 (2017)
  * [18] Mahmud, S.H., Chen, W., Meng, H., Jahan, H., Liu, Y., Hasan, S.M.: Prediction of drug-target interaction based on protein features using undersampling and feature selection techniques with boosting. Analytical biochemistry 589, 113507 (2020)
  * [19] Mohamed, S.K., Nováček, V., Nounu, A.: Discovering protein drug targets using knowledge graph embeddings. Bioinformatics 36(2), 603–610 (2020)
  * [20] Opella, S.J.: Structure determination of membrane proteins by nuclear magnetic resonance spectroscopy. Annual review of analytical chemistry (Palo Alto, Calif.) 6, 305 (2013)
  * [21] Pan, J., Li, L.P., You, Z.H., Yu, C.Q., Ren, Z.H., Chen, Y.: Prediction of drug–target interactions by combining dual-tree complex wavelet transform with ensemble learning method. Molecules 26(17), 5359 (2021)
  * [22] Paul, S.M., Mytelka, D.S., Dunwiddie, C.T., Persinger, C.C., Munos, B.H., Lindborg, S.R., Schacht, A.L.: How to improve r&d productivity: the pharmaceutical industry’s grand challenge. Nature reviews Drug discovery 9(3), 203–214 (2010)
  * [23] Prokhorenkova, L., Gusev, G., Vorobev, A., Dorogush, A.V., Gulin, A.: Catboost: unbiased boosting with categorical features. Advances in neural information processing systems 31 (2018)
  * [24] Sachdev, K., Gupta, M.K.: A comprehensive review of feature based methods for drug target interaction prediction. Journal of biomedical informatics 93, 103159 (2019)
  * [25] Sajadi, S.Z., Zare Chahooki, M.A., Gharaghani, S., Abbasi, K.: Autodti++: deep unsupervised learning for dti prediction by autoencoders. BMC bioinformatics 22(1), 1–19 (2021)
  * [26] Schomburg, I., Chang, A., Ebeling, C., Gremse, M., Heldt, C., Huhn, G.: & schomburg, d.(2004). brenda, the enzyme database: updates and major new developments. Nucleic Acids Research D 32
  * [27] Shi, H., Liu, S., Chen, J., Li, X., Ma, Q., Yu, B.: Predicting drug-target interactions using lasso with random forest based on evolutionary information and chemical structure. Genomics 111(6), 1839–1852 (2019)
  * [28] Wan, F., Hong, L., Xiao, A., Jiang, T., Zeng, J.: Neodti: neural integration of neighbor information from a heterogeneous network for discovering new drug–target interactions. Bioinformatics 35(1), 104–111 (2019)
  * [29] Wang, C., Wang, W., Lu, K., Zhang, J., Chen, P., Wang, B.: Predicting drug-target interactions with electrotopological state fingerprints and amphiphilic pseudo amino acid composition. International Journal of Molecular Sciences 21(16), 5694 (2020)
  * [30] Wang, L., Wong, L., Chen, Z.H., Hu, J., Sun, X.F., Li, Y., You, Z.H.: Mspedti: Prediction of drug–target interactions via molecular structure with protein evolutionary information. Biology 11(5), 740 (2022)
  * [31] Wang, Y., Wang, L., Wong, L., Zhao, B., Su, X., Li, Y., You, Z.: Rofdt: Identification of drug–target interactions from protein sequence and drug molecular structure using rotation forest. Biology 11(5), 741 (2022)
  * [32] Wishart, D.S., Knox, C., Guo, A.C., Cheng, D., Shrivastava, S., Tzur, D., Gautam, B., Hassanali, M.: Drugbank: a knowledgebase for drugs, drug actions and drug targets. Nucleic acids research 36(suppl_1), D901–D906 (2008)
  * [33] Yamanishi, Y., Araki, M., Gutteridge, A., Honda, W., Kanehisa, M.: Prediction of drug–target interaction networks from the integration of chemical and genomic spaces. Bioinformatics 24(13), i232–i240 (2008)
  * [34] Yamanishi, Y., Kotera, M., Kanehisa, M., Goto, S.: Drug-target interaction prediction from chemical, genomic and pharmacological data in an integrated framework. Bioinformatics 26(12), i246–i254 (2010)
  * [35] Zheng, X., Ding, H., Mamitsuka, H., Zhu, S.: Collaborative matrix factorization with multiple similarities for predicting drug-target interactions. In: Proceedings of the 19th ACM SIGKDD international conference on Knowledge discovery and data mining. pp. 1025–1033 (2013)


remarkRemark hypothesisHypothesis claimClaim HDG methods for solving the two-
fluid plasma modelA. Ho and U. Shumlak [][nocite]ex_supplement

# Hybridizable discontinuous Galerkin methods for solving the two-fluid plasma
model††thanks: Submitted to the editors December 1, 2023. The information,
data, or work presented herein was funded in part by the Air Force Office of
Scientific Research under award numbers FA9550-15-1-0271. This material is
also based upon work supported by the National Science Foundation under Grant
No. PHY-2108419.

Andrew Ho Computational Plasma Dynamics Lab, Aerospace and Energetics Research
Program, University of Washington, Seattle, WA. Uri Shumlak22footnotemark: 2

###### Abstract

The two-fluid plasma model has a wide range of timescales which must all be
numerically resolved regardless of the timescale on which plasma dynamics
occurs. The answer to solving numerically stiff systems is generally to
utilize unconditionally stable implicit time advance methods. Hybridizable
discontinuous Galerkin (HDG) methods have emerged as a powerful tool for
solving stiff partial differential equations. The HDG framework combines the
advantages of the discontinuous Galerkin (DG) method, such as high-order
accuracy and flexibility in handling mixed hyperbolic/parabolic PDEs with the
advantage of classical continuous finite element methods for constructing
small numerically stable global systems which can be solved implicitly. In
this research we quantify the numerical stability conditions for the two-fluid
equations and demonstrate how HDG can be used to avoid the strict stability
requirements while maintaining high order accurate results.

###### keywords:

finite elements, hybridizable discontinuous Galerkin, computational plasma
physics, two-fluid plasma models

35L60, 65M12, 65M60

##  1 Introduction

The two-fluid plasma model is a very numerically stiff PDE system. The model
couples together the time evolution of an ion fluid, an electron fluid, and
Maxwell’s equations. As a result the characteristic dynamics leads to
timescales which can span two to four orders of magnitude, include both the
speed of light and various electron wave modes along with ion fluid dynamics.
The large span of timescales makes the PDE system a prime candidate for being
solved using implicit time stepping. However, there are several challenges in
actualizing an implicit numerical solution. The system is a mixed advection-
diffusion-reaction system consisting of six scalar fields and four vector
fields with a nominal 18 components to be evolved in 3D. The hyperbolic
characteristics makes using standard continuous Galerkin (CG) methods
challenging, and traditional discontinuous Galerkin (DG) methods are not
typically amenable with implicit time stepping schemes[20]. One solution is
the use of hybridizable discontinuous Galerkin (HDG) methods[4]. These methods
are designed to allow the implicit DG discretization to be partitioned into
three distinct phases, reducing the largest global system solve size as well
as producing a scheme which maintains local compactness. Nguyen and
Peraire[19] have demonstrated the applicability of HDG to a wide variety of
linear and non-linear PDEs, including Stoke’s equations, incompressible
Navier-Stokes equations, and the compressible Euler equations. They
demonstrated that the method can be formulated to achieve optimal $k+1$ orders
of convergence, where $k$ is the polynomial basis order. Lee et al.[15]
demonstrated that the linearized incompressible resistive MHD equations can be
solved using HDG, opening the door to modeling basic plasma systems using HDG.
In this paper we present techniques and tools for expanding the applicability
of HDG to large coupled PDE systems such as the two-fluid plasma model using
HDG. Section 2 introduces the two-fluid plasma system to be solved, and
discusses challenges associated with solving this particular system using
explicit Runge-Kutta discontinuous Galerkin (RKDG). LABEL:sec:hdg presents the
HDG discretization method and practical implementation details required for
solving complex PDE systems. Section 4 presents the verification of the
implementation by checking the convergence rates on linear systems. Finally
Section 5 presents results of solving the two-fluid plasma using HDG.

##  2 The two-fluid plasma model

The two-fluid plasma model is a subset of the general multi-fluid 5N-moment
plasma model[22]. This model treats each constituent particle species type as
their own fluid species. Coupling between the species occurs via local inter-
species collisions and long-range electromagnetic field coupling described
using Maxwell’s Equations. For the two-fluid model a single positively charge
ion species is coupled to an electron species. The exact form being considered
for this paper is given by

(1)

where $\rho_{\alpha}$, $\vec{p}_{\alpha}$, $U_{\alpha}$ are the mass,
momentum, and internal energy densities of species $\alpha$ and $\mu_{c}$,
$\kappa_{c}$, $R_{c}$, $Q_{c}$ are constants for the viscosity, thermal
diffusivity, interspecies momentum transfer, and interspecies heat transfer
coefficients respectively. $m_{\alpha}$, $Z_{\alpha}$, and $\gamma$ describe
the mass, charge state, and heat capacity ratio of species $\alpha$. For the
two-fluid model $\alpha$ is either $i$ for ions or $e$ for electrons, while
$\beta$ is the other species. The following additional substitutions are used
to simplify the notation:

(2) $\displaystyle\vec{u}_{\alpha}=\frac{\vec{p}_{\alpha}}{\rho_{\alpha}}$ (3)
$\displaystyle P_{\alpha}=U_{\alpha}(\gamma-1)$ (4) $\displaystyle
T_{\alpha}=\frac{m_{\alpha}P_{\alpha}}{\rho_{\alpha}}$ (5)
$\displaystyle\rho_{c}=\frac{Z_{i}}{m_{i}}\rho_{i}+\frac{Z_{e}}{m_{e}}\rho_{e}$
(6)
$\displaystyle\vec{j}=\frac{Z_{i}}{m_{i}}\vec{p}_{i}+\frac{Z_{e}}{m_{e}}\vec{p}_{e}$

The plasma fluid equations are then coupled to EM fields via Maxwell’s
equations. These are evolved using For a large PDE system such as the two-
fluid plasma model, deriving correct analytical Jacobians by hand is
practically impossible due to the large number of non-zero entries. For the
two-fluid plasma system this requires at least 1047 unique non-zero entries,
before taking into account boundary conditions which may double or triple the
number of analytical Jacobian non-zeroes. To solve this problem the Jacobian
along with the rest of the HDG code is generated using SymPy[17], a symbolic
computer algebra system. The tool is capable of automatically detecting and
eliminating unused degrees of freedom and handles multiple subdomains/boundary
conditions. The generated HDG code is coupled with a diagonally implicit
Runge-Kutta[11] temporal solver, and parallelized on CPUs using MPI.

###  3.2 Face-oriented numerical flux

A key requirement for the HDG numerical flux is that from the perspective of
any element the numerical flux must be a function of only the DOFs inside the
element and the DOF on the neighboring face. Each face is assigned a local
normal, tangent, and binormal. As a result the face normal will point outwards
for one element, while for the adjacent neighboring element it points inwards.
The non-uniqueness of the tangent and binormal directions is not an issue so
long as they are chosen such that $\hat{t}=\hat{n}\times\hat{b}$. Vector and
tensor quantities are appropriately rotated to the face coordinate system for
computing numerical fluxes.

Figure 3: Global and face-oriented coordinate systems

Computing the numerical flux becomes

(58)
$\displaystyle(\stackrel{{\scriptstyle\leftrightarrow}}{{\boldsymbol{T}}}\hat{n})^{*}=\pm\vec{T}_{n}^{*}$

Depending on the given system this could potentially lead to a reduction of
face DOF[1]. For example, take the linear diffusion equation

$\displaystyle\partial_{t}\phi-\operatorname{\nabla{\cdot}}\vec{\gamma}=$ (59)
$\displaystyle\vec{\gamma}=\operatorname{\nabla}\phi$

The numerical fluxes chosen are the class of central fluxes of the form

(60)
$\displaystyle(\stackrel{{\scriptstyle\leftrightarrow}}{{\boldsymbol{T}}}\hat{n})^{*}=\stackrel{{\scriptstyle\leftrightarrow}}{{\boldsymbol{T}}}_{\lambda}\hat{n}+\tau(\vec{q}^{-}-\vec{\lambda})$

For the linear diffusion equation this expands to

(61)
$\displaystyle(\stackrel{{\scriptstyle\leftrightarrow}}{{\boldsymbol{T}}}\hat{n})^{*}=\begin{bmatrix}\lambda_{\vec{\gamma}}\cdot\hat{n}^{-}+\tau(q^{-}_{\phi}-\lambda_{\phi})\\\
-\lambda_{\phi}\stackrel{{\scriptstyle\leftrightarrow}}{{\boldsymbol{I}}}\hat{n}^{-}+\tau(q^{-}_{\vec{\gamma}}-\lambda_{\vec{\gamma}})\end{bmatrix}$

where $\tau>0$ is a user chosen constant. Note that in practice so long as
$\tau\approx 1$ the method will converge at the same rate. In a global
coordinate system this expands to

(62)
$\displaystyle(\stackrel{{\scriptstyle\leftrightarrow}}{{\boldsymbol{T}}}\hat{n})^{*}=\begin{bmatrix}\lambda_{\gamma_{x}}\hat{n}^{-}_{x}+\lambda_{\gamma_{y}}\hat{n}^{-}_{y}+\lambda_{\gamma_{z}}\hat{n}^{-}_{z}+\tau(q^{-}_{\phi}-\lambda_{\phi})\\\
\begin{bmatrix}-\lambda_{\phi}\hat{n}^{-}_{x}+\tau(q^{-}_{\gamma_{x}}-\lambda_{\gamma_{x}})\\\
-\lambda_{\phi}\hat{n}^{-}_{y}+\tau(q^{-}_{\gamma_{y}}-\lambda_{\gamma_{y}})\\\
-\lambda_{\phi}\hat{n}^{-}_{z}+\tau(q^{-}_{\gamma_{z}}-\lambda_{\gamma_{z}})\end{bmatrix}\end{bmatrix}$

which requires four element scalar fields and four face scalar fields. In a
rotated frame this can be simplified to

(63)
$\displaystyle(\stackrel{{\scriptstyle\leftrightarrow}}{{\boldsymbol{T}}}\hat{n})^{*}=\begin{bmatrix}\pm\lambda_{\gamma_{n}}+\tau(q^{-}_{\phi}-\lambda_{\phi})\\\
\begin{bmatrix}\mp\lambda_{\phi}+\tau(q^{-}_{\gamma_{n}}-\lambda_{\gamma_{n}})\\\
0\\\ 0\end{bmatrix}\end{bmatrix}$

eliminating two face scalar fields. For the two-fluid plasma model with either
mixed parabolic-hyperbolic cleaning or purely parabolic cleaning this rotation
scheme requires a total of 42 element and 34 face scalar fields, eliminating
eight face scalar fields from the non-rotated scheme.

##  4 Validation of convergence rates for various linear systems

Convergence is evaluated for several model PDE systems including the linear
advection, diffusion equation, and wave equation. In all cases periodic
boundary conditions are applied on a 3D cube domain $\vec{r}\in[0,1]^{3}$ and
$t\in[0,0.1]$ with a fixed $\Delta t=5\times 10^{-3}$. The error is compared
against the analytical solution $a_{a}$ using an $8^{3}$ tensor product Gauss-
Lobatto quadrature to compute the L-2 error integral

(64) $\displaystyle\epsilon=\sqrt{\int(q_{h}-q_{a})^{2}dV}$

The linear advection system is defined by

$\displaystyle\partial_{t}q+\operatorname{\nabla{\cdot}}(\vec{a}q)=$ (65)
$\displaystyle\vec{a}^{T}=\begin{bmatrix}1&\frac{1}{2}&2\end{bmatrix}$ (66)
$\displaystyle
q_{a}(x,y,z,t)=\cos(2\pi(x-a_{x}t))\cos(2\pi(y-a_{y}t))\cos(2\pi(z-a_{z}t))$

Figure 4 demonstrates that the linear advection system converges at the
expected optimal $O(N+1)$ rates.

Figure 4: L-2 error norms for the linear advection equation converge at the
expected optimal $O(N+1)$ rates for degree $N$ basis. The observed $10^{-5}$
minimum error is from the temporal discretization error.

Note that at small $\Delta x$ and high basis order the simulation reaches a
saturation floor of $10^{-5}$ from the temporal discretization error. This can
be reduced by choosing a smaller $\Delta t$ or more accurate temporal
discretization method.

The linear diffusion system is defined by

$\displaystyle\partial_{t}q-\operatorname{\nabla{\cdot}}(k\operatorname{\nabla}q)=$

This can be decomposed into a coupled system of first order equations to fit
the form of LABEL:eq:hdg_prototype as

$\displaystyle\partial_{t}q-\operatorname{\nabla{\cdot}}(k\vec{u})=$ (67)
$\displaystyle\vec{u}-\operatorname{\nabla}q=0$ (68) $\displaystyle k=10^{-2}$
(69) $\displaystyle q_{a}(x,y,z,t)=\exp\left(-3(2\pi)^{2}kt\right)\sin(2\pi
x)\sin(2\pi x)\sin(2\pi z)$

where initial conditions are given by $q(x,y,z)=q_{a}(x,y,z,t=0)$. Figure 5
demonstrates optimal $O(N+1)$ convergence rates for $q$, however $\vec{u}$
converges with a mix of $O(N+1)$ and $O(N)$ rates depending if $N$ is odd or
even. This oscillating convergence rate is a known behavior of central
fluxes[12, 6]. Implementing a bidirectional upwinding flux for handling second
order spatial derivatives has been observed to resolve the sub-optimal
convergence rates for even $N$[4].

(a) $q$

(b) $u_{x}$

Figure 5: L-2 error norms for the linear diffusion equation. Note that $q$
converges at $O(N+1)$, while $\vec{u}$ has an odd-even oscillation between
$O(N+1)$ and $O(N)$ convergence rates. $u_{y}$ and $u_{z}$ have near identical
plots to $u_{x}$ and are omitted for brevity.

The linear wave equation is defined by

(70) $\displaystyle\partial_{t}q-v=0$ (71)
$\displaystyle\partial_{t}\vec{u}-\operatorname{\nabla}v=0$ (72)
$\displaystyle c=1$ (73) $\displaystyle q_{a}(x,y,z,t)=\cos\left(2\pi
t\sqrt{3c^{2}}\right)\sin(2\pi x)\sin(2\pi y)\sin(2\pi z)$

We can recover the traditional wave equation by solving Eq. 70 for $v$ and
substituting into Eq. 71 and LABEL:eq:wave3:

$\displaystyle\partial_{t}\vec{u}-\operatorname{\nabla}\partial_{t}q=0$
$\displaystyle\begin{aligned}
\partial_{t}(\vec{u}-\operatorname{\nabla}q)&=0&\rightarrow&&\vec{u}&=\operatorname{\nabla}q\end{aligned}$
$\displaystyle\partial_{tt}q-\operatorname{\nabla{\cdot}}\left(c^{2}\operatorname{\nabla}q\right)=$

We choose to not solve the wave equation in this form because of the presence
of second order time and spatial derivatives, which our implementation cannot
directly handle. Figure 6 demonstrates that all variables present in the
simulation converge at the optimal $O(N+1)$ rate. Unlike the diffusion
equation the use of central fluxes does not negatively impact convergence
rates despite using the same splitting method for the Laplacian operator[5].

(a) $q$

(b) $v$

(c) $u_{x}$

Figure 6: L-2 error norms for the linear wave equation. All variables converge
at the optimal $O(N+1)$ rate. Plots for $u_{y}$ and $u_{z}$ are nearly
identical to $u_{x}$ and omitted.

##  5 Application of HDG to the two-fluid plasma system

The two-fluid plasma system was used to test the behavior of the HDG code on
large coupled systems. A quasi-1D domain of $256\times 1\times 1$ parabolic
elements with $x\in[-5,5]$ is initialized with a two-fluid extension of the
magnetized Brio-Wu shock tube[23, 16] such that

(74) $\displaystyle\begin{aligned} \begin{bmatrix}n_{i}\\\ U_{i}\\\ n_{e}\\\
U_{e}\\\ B_{x}\\\ B_{y}\end{bmatrix}_{\text{left}}&=\begin{bmatrix}1\\\
7.5\times 10^{-5}\\\ 1\\\ 7.5\times 10^{-5}\\\ 7.5\times 10^{-3}\\\
10^{-2}\end{bmatrix}&\begin{bmatrix}n_{i}\\\ U_{i}\\\ n_{e}\\\ U_{e}\\\
B_{x}\\\ B_{y}\end{bmatrix}_{\text{right}}&=\begin{bmatrix}0.125\\\ 7.5\times
10^{-5}\\\ 0.125\\\ 7.5\times 10^{-6}\\\ 7.5\times 10^{-3}\\\
-10^{-2}\end{bmatrix}\end{aligned}$ (75) $\displaystyle\begin{aligned}
\frac{\delta_{p}}{L}&=1&\gamma&=\frac{5}{3}&m_{i}&=1&m_{e}&=\frac{1}{1836}&\frac{c_{0}}{V_{A}}&=c_{h}=c_{p}=1&Z_{i}&=-Z_{e}=10\end{aligned}$
(76) $\displaystyle\mu_{c}=\kappa_{c}=Q_{c}=R_{c}=10^{-5}$

The is simulated from $t\in[0,100]$. Figure 7 shows various plasma properties
at the end of the simulation. The plotted fields show similar features to the
work of Shumlak and Loverich[23], while being distinct from results produced
using standard MHD. Some examples of these differences is the presence of a
Whistler wave propagating to the left around $x=1.75$, and the solution
producing $T_{i}\neq T_{e}$.

(a) $\rho$

(b) $B_{y}$

(c) $T$

(d) $U$

Figure 7: Two-fluid Brio-Wu shocktube properties at $t=100$. Similar features
to those found in Shumlak and Loverich[23] are found, which are notably
distinct from the classical MHD results.

To validate the correct propagation speed of various important two-fluid
waves, a second simulation with small jumps is initialized to ensure linear
behavior. The domain is changed to $1024\times 1\times 1$ linear elements with
$x\in[-\pi,\pi]$, and

(77) $\displaystyle\begin{aligned} \begin{bmatrix}n_{i}\\\ U_{i}\\\ n_{e}\\\
U_{e}\\\ B_{x}\\\ B_{y}\end{bmatrix}_{\text{left}}&=\begin{bmatrix}10.0005\\\
7.500375\times 10^{-5}\\\ 10.0005\\\ 7.500375\times 10^{-5}\\\ 7.5\times
10^{-2}\\\ 5\times 10^{-5}\end{bmatrix}&\begin{bmatrix}n_{i}\\\ U_{i}\\\
n_{e}\\\ U_{e}\\\ B_{x}\\\
B_{y}\end{bmatrix}_{\text{right}}&=\begin{bmatrix}9.9995\\\ 7.499625\times
10^{-5}\\\ 9.9995\\\ 7.499625\times 10^{-5}\\\ 7.5\times 10^{-2}\\\ -5\times
10^{-5}\end{bmatrix}\end{aligned}$ (78) $\displaystyle\begin{aligned}
\frac{\delta_{p}}{L}&=1&\gamma&=\frac{5}{3}&m_{i}&=1&1836m_{e}&=1&\frac{c_{0}}{V_{A}}&=c_{h}=c_{p}=1&Z_{i}&=-Z_{e}=1\end{aligned}$
(79) $\displaystyle\mu_{c}=\kappa_{c}=Q_{c}=R_{c}=10^{-5}$

A mixed 2D spatial/temporal Fourier analysis is then performed on $E_{y}$[23].
This is plotted against the theoretical L and R mode wave propagation speeds,
described by the dispersion relation

(80)
$\displaystyle\frac{c_{0}^{2}k^{2}}{V_{A}^{2}\omega^{2}}=1-\frac{\omega_{pe}^{2}}{\omega(\omega\pm\omega_{ce})}$
(81) $\displaystyle\begin{aligned}
\omega_{pe}&=(\omega_{p}\tau)\sqrt{\frac{n_{e}Z_{e}^{2}}{m_{e}}}&\omega_{ce}&=\frac{L}{\delta_{p}}\frac{|Z_{e}|B}{m_{e}}\end{aligned}$

Figure 8: Fourier analysis of $E_{y}$ matches the circularly polarized two-
fluid plasma dispersion relation.

Figure 8 demonstrates that the Fourier analysis has three dominant modes
corresponding to the three branches of the L and R modes. This test
demonstrates the viability of using HDG to solve large-scale coupled plasma
differential equations and the validity of the implemented tools to ensure
correctness and simplify implementing such solvers.

Figure 9 shows the numerical stability conditions for an equivalent explicit
explicit RK4. For analysis purposes the implicit scheme used $\Delta
t=6.28\times 10^{-3}$. While an explicit RK4 DG solver can run a $N=0$ basis
with this timestep, higher degree basis require a smaller timestep. Define the
timestep ratio gain as the implicit timestep $\Delta t$ used divided by the
maximum stable explicit timestep. In this case there is a timestep ratio gain
of $2.2$ for $N=1$ and $4.4$ for $N=2$ over an explicit RKDG method. The
primary reason these ratios are low is that the original test conditions were
setup such that it is feasible to run and benchmark an explicit RKDG code.
Nevertheless, it is sufficient for demonstrating that the HDG method can
continue to operate accurately in a regime where explicit time stepping
methods become unstable. In other test problems we have observed stable and
accurate runs using timesteps greater than $60$ times larger than allowable
for the explicit RK4 DG method.

Figure 9: Stability conditions for a two-fluid explicit RK4 DG solver for the
above specified conditions. The red highlighted region is the region of
absolute stability for RK4 and the curves define the stability eigencurves of
the ideal two-fluid system generated using the method from Section 2. Implicit
simulations were run with $\Delta t=6.28\times 10^{-3}$ to achieve sufficient
temporal resolution for the Fourier dispersion analysis in Fig. 8. A degree 0
basis could be run explicitly with this timestep; however, an explicit
solution with a high degree basis is unstable using this timestep.

##  6 Conclusion

This work increases the viability of using the two-fluid plasma model, which
while more physically accurate than single fluid models is numerically
stiff[24]. We perform a detailed numerical stability analysis of solving the
ideal two-fluid system using classical RKDG. The numerically stiffness of the
two-fluid plasma model results from needing to temporally resolve the speed of
light, $c_{0}/V_{A}$, and the eigenvalues of the source terms,
$\stackrel{{\scriptstyle\leftrightarrow}}{{\boldsymbol{J}}}_{S}$, when bulk
plasma motion is dictated by timescales on the order of the ion sound speed.
We also demonstrate that hybridizable discontinuous Galerkin enables solving
large-scale PDE systems with implicit temporal solvers using a systematic
formulation to allow automatic generation of the HDG Schur complement. This
systematic process is shown to handle over 40 coupled differential equations
and to accurately compute the thousands of analytical non-zero terms of
Jacobian, as well as to identify and eliminate variables which analytically
vanish. Potential future extensions include implementing bidirectional
upwinding numerical fluxes[4] to improve convergence of parabolic/elliptical
terms, as well as investigate effective techniques for reducing the time it
takes to construct and solve the HDG Schur complement system.

## Supplementary Material

A copy of the code can be obtained from
https://bitbucket.org/helloworld922/dg.

## References

  * [1] T. Bui-Thanh, Construction and analysis of hdg methods for linearized shallow water equations, SIAM Journal on Scientific Computing, 38 (2016), pp. A3696–A3719.
  * [2] N. Chalmers and L. Krivodonova, A robust cfl condition for the discontinuous galerkin method on triangular meshes, Journal of Computational Physics, 403 (2020), p. 109095, https://doi.org/https://doi.org/10.1016/j.jcp.2019.109095, https://www.sciencedirect.com/science/article/pii/S0021999119308009.
  * [3] J. G. Charney, R. FjÖrtoft, and J. V. Neumann, Numerical integration of the barotropic vorticity equation, Tellus, 2 (1950), pp. 237–254, https://doi.org/10.3402/tellusa.v2i4.8607, https://doi.org/10.3402/tellusa.v2i4.8607, https://arxiv.org/abs/https://doi.org/10.3402/tellusa.v2i4.8607.
  * [4] B. Cockburn, B. Dong, J. Guzman, M. Restelli, and R. Sacco, A hybridizable discontinuous galerkin method for steady-state convection-diffusion-reaction problems, SIAM Journal on Scientific Computing, 31 (2009), pp. 3827–3846.
  * [5] B. Cockburn, N. Nguyen, and J. Peraire, Chapter 8 - hdg methods for hyperbolic problems, in Handbook of Numerical Methods for Hyperbolic Problems, R. Abgrall and C.-W. Shu, eds., vol. 17 of Handbook of Numerical Analysis, Elsevier, 2016, pp. 173–197, https://doi.org/https://doi.org/10.1016/bs.hna.2016.07.001, https://www.sciencedirect.com/science/article/pii/S1570865916300035.
  * [6] B. Cockburn and C.-W. Shu, The local discontinuous galerkin method for time-dependent convection-diffusion systems, SIAM Journal on Numerical Analysis, 35 (1997), pp. 2440–2463.
  * [7] A. Dedner, F. Kemm, D. Kroner, C.-D. Munz, T. Schnitzer, and M. Wesenberg, Hyperbolic divergence cleaning for the mhd equations, Journal of Computational Physics, 175 (2002), pp. 645–673.
  * [8] J. E. Dennis and R. B. Schnabel, Numerical Methods for Unconstrained Optimization and Nonlinear Equations, Society for Industrial and Applied Mathematics, 1996, https://doi.org/10.1137/1.9781611971200, https://epubs.siam.org/doi/abs/10.1137/1.9781611971200, https://arxiv.org/abs/https://epubs.siam.org/doi/pdf/10.1137/1.9781611971200.
  * [9] A. H. Glasser and X. Z. Tang, The sel macroscopic modeling code, Computer Physics Communications, (2004), pp. 237–243.
  * [10] R. J. Guyan, Reduction of stiffness and mass matrices, AIAA Journal, 3 (1965), pp. 380–380, https://doi.org/10.2514/3.2874, https://doi.org/10.2514/3.2874, https://arxiv.org/abs/https://doi.org/10.2514/3.2874.
  * [11] E. Hairer and G. Wanner, Solving Ordinary Differential Equations II, Springer-Verlag, 1980.
  * [12] J. S. Hesthaven and T. Warburton, Nodal Discontinuous Galerkin Methods, Springer, 2008.
  * [13] R. M. Kirby, S. J. Sherwin, and B. Cockburn, To cg or to hdg: A comparative study, Journal of Scientific Computing, 51 (2012), pp. 183–212.
  * [14] F. F.-K. Kuo, Network Analysis and Synthesis, John Wiley & Sons Inc, 2nd ed., 1966.
  * [15] J. J. Lee, S. Shannon, an. Bui-Thanh, and J. N. Shadid, Analysis of an hdg method for linearized incompressible resistive mhd equations, SIAM Journal of Numerical Analysis, 57 (2019).
  * [16] J. Loverich, A. Hakim, and U. Shumlak, A discontinuous galerkin method for ideal two-fluid plasma equations, Communications in Computational Physics, 9 (2011), pp. 240–268.
  * [17] A. Meurer, C. P. Smith, M. Paprocki, O. Čertík, S. B. Kirpichev, M. Rocklin, A. Kumar, S. Ivanov, J. K. Moore, S. Singh, T. Rathnayake, S. Vig, B. E. Granger, R. P. Muller, F. Bonazzi, H. Gupta, S. Vats, F. Johansson, F. Pedregosa, M. J. Curry, A. R. Terrel, v. Roučka, A. Saboo, I. Fernando, S. Kulal, R. Cimrman, and A. Scopatz, Sympy: symbolic computing in python, PeerJ Computer Science, 3 (2017), p. e103, https://doi.org/10.7717/peerj-cs.103, https://doi.org/10.7717/peerj-cs.103.
  * [18] C.-D. Munz, P. Ommes, and R. Schneider, A three-dimensional finite-volume solver for the maxwell equations with divergence cleaning on unstructured meshes, Computer Physics Communications, 130 (2000), pp. 83–117.
  * [19] N. C. Nguyen and J. Peraire, Hybridizable discontinuous galerkin methods for partial differential equations in continuum mechanics, Journal of Computational Physics, 231 (2012), pp. 5955–5988.
  * [20] J. Peraire and P.-O. Persson, The compact discontinuous galerkin (cdg) method for elliptic problems, SIAM Journal on Scientific Computing, 30 (2008), pp. 1806–1824.
  * [21] W. Schiesser, The Numerical Method of Lines, Academic Press, 1st ed., 1991.
  * [22] U. Shumlak, R. Lilly, N. Reddell, E. Sousa, and B. Srinivasan, Advanced physics calculations using a multi-fluid plasma model, Computer Physics Communications, 182 (2011), pp. 1767–1770, https://doi.org/https://doi.org/10.1016/j.cpc.2010.12.048, https://www.sciencedirect.com/science/article/pii/S001046551100004X.  Computer Physics Communications Special Edition for Conference on Computational Physics Trondheim, Norway, June 23-26, 2010.
  * [23] U. Shumlak and J. Loverich, Approximate riemann solver for the two-fluid plasma model, Journal of Computational Physics, 187 (2003), pp. 620–638.
  * [24] B. Srinivasan and U. Shumlak, Analytical and computational study of the ideal full two-fluid plasma model and asymptotic approximations for Hall-magnetohydrodynamics, Physics of Plasmas, 18 (2011), p. 092113, https://doi.org/10.1063/1.3640811, https://doi.org/10.1063/1.3640811, https://arxiv.org/abs/https://pubs.aip.org/aip/pop/article-pdf/doi/10.1063/1.3640811/13684706/092113_1_online.pdf.
  * [25] T. Warburton and T. Hagstrom, Taming the cfl number for discontinuous galerkin methods on structured meshes, SIAM Journal on Numerical Analysis, 46 (2008), pp. 3151–3180, https://doi.org/10.1137/060672601, https://doi.org/10.1137/060672601, https://arxiv.org/abs/https://doi.org/10.1137/060672601.