Assessing and Predicting Drug-Induced Anticholinergic Risks: An Integrated Computational Approach
DOI: 10.1177/2042098617725267
Xu D, Anderson HD, Tao A, Hannah KL, Linnebur SA, Valuck RJ, Culbertson VL,
Therapeutic Advances in Drug Safety, 2017.

Background: Anticholinergic (AC) adverse drug events (ADE) are caused by inhibition of muscarinic receptors as a result of designated or off-target drug-receptor interactions. In practice, AC toxicity is assessed primarily based on clinician experience. Objective: The goal is to evaluate a novel concept of integrating big pharmacological and healthcare data to assess clinical AC toxicity risks. Methods: Anticholinergic toxicity scores (ATS) were computed using drug-receptor inhibitions identified through pharmacological data screening. A longitudinal retrospective cohort study using medical claims data was performed to quantify AC clinical risks. ATS was compared to two previously reported toxicity measures. A quantitative structure-activity relationship (QSAR) model was established for rapid assessment and prediction of AC clinical risks. Results: 25 common medications, 575,228 exposed and unexposed patients were analyzed. Our data indicated that A
TS is more consistent with the trend of AC outcomes than other toxicity methods. Incorporating drug pharmacokinetic parameters to ATS yielded a QSAR model with excellent correlation to AC incident rate (R2=0.83) and predictive performance (cross validation Q2=0.64). Good correlation and predictive performance (R2=0.68 / Q2=0.29) were also obtained for an M2 receptor-specific QSAR model and Tachycardia, an M2 receptor-specific ADE. Conclusions: Albeit using a small medication sample size, our pilot data demonstrated the potential and feasibility of a new computational AC toxicity scoring approach driven by underlying pharmacology and big data analytics. Follow-up work is under way to further develop the ATS scoring approach and clinical toxicity predictive model using a large number of medications and clinical parameters.

Experimental Identification and Computational Characterization of a Novel Extracellular Metalloproteinase Produced by Clostridium sordellii
DOI: 10.1039/C6RA27654G
Aldape MJ, Tao A, Heeney DD, McIndoo ER, French JM, Xu D. RSC Advances, 2017, 7, 13928-13938.

Clostridium sordellii is a lethal pathogen for both animals and humans. Severe capillary leakage, toxic shock syndrome, and an extreme leuke moid reaction (LR), are hallmark features of C. sordellii infections and contribute to its high mortality rate. Here we report the discovery of a previously unknown and uncharacterized metalloproteinase of C. sordellii (referred as Mcs1) that cleaves human vascular cell adhesion molecule (VCAM)-1 in vitro, an adhesion molecule critical to hematopoietic precursor retention and l eukocyte diapedesis. We successfully identified the open reading frame encoding Mcs1 within the ATCC 9714 genome and developed an Δmcs1 mutant strain using the ClosTron mutagenesis technology. No VCAM-1 proteolysis was observed from exotoxins collected from mutant strain cultures. Using advanced protein structural modeling and molecular dynamics simulation techniques, the 3D molecular structure and conformational features of Mcs1 were also characterized. Our data demonstrates that Mcs1 proteolytic activity is controlled by the electrostatic interactions between Glu113 and Arg227 residues and the gating motions within its cleft region. This pilot interdisciplinary investigation provided crucial experimental evidence of the existence of Mcs1 in C. sordellii and molecular insights into its 3D structure and proteolytic activity. These findings have the potential to help advance new therapeutics and diagnostics against deadly C. sordellii infections. Follow-up in vitro and in vivo work is under way to further characterize Mcs1 enzymatic kinetics and its role in C. sordellii pathogenesis.

Routine Microsecond Molecular Dynamics Simulations with AMBER on GPUs
  DOI: 10.1021/ct200909j
Götz AW, Williamson MJ, Xu D, Poole D, Le Grand S, Walker RC
, J. Chem. Theory Comput. 2012, 8(5): 1542–1555.

We present an implementation of generalized Born implicit solvent all-atom classical molecular dynamics (MD) within the AMBER program package that runs entirely on CUDA enabled NVIDIA graphics processing units (GPUs). We discuss the algorithms that are used to exploit  the processing power of the GPUs and show the performance that can be achieved in comparison  to simulations on conventional CPU clusters. The implementation supports three different  precision models in which the contributions to the forces are calculated in single precision  floating point arithmetics but accumulated in double precision (SPDP), or everything is computed  in single precision (SPSP), or double precision (DPDP). In addition to performance we  have focused on understanding the implications of the different precision models on the outcome  of implicit solvent MD simulations. We show results for a range of tests including the  accuracy of single point force evaluations and energy conservation as well as structural properties  pertainining to protein dynamics. The numerical noise due to rounding errors within the  SPSP precision model is sufficiently large to lead to an accumulation of errors which can result  in unphysical trajectories for long timescale simulations. We recommend the use of the mixedprecision  SPDP model since the numerical results obtained are comparable with those of the full double precision DPDP model and the reference double precision CPU implementation  but at significantly reduced computational cost. Our implementation provides performance for GB simulations in a single desktop that is on par with, and in some cases exceeds that of  traditional supercomputers.

Inhibition of Snowshoe Hare Succinate Dehydrogenase Activity as a Mechanism of Deterrence for Papyriferic Acid in Birch DOI: 10.1007/s10886-011-0039-9
Forbey JS, Pu X, Xu D, Kielland K, Bryant J. J. Chem. Ecol., 2011, 137(12):1285-1293.

The plant secondary metabolite papyriferic acid (PA) deters browsing by snowshoe hares (Lepus americanus) on the juvenile developmental stage of the Alaska paper birch (Betula neoalaskana). However, the physiological mechanism that reduces browsing remains unknown. We used pharmacological assays and molecular modeling to test the hypothesis that inhibition of succinate dehydrogenase (SDH) is a mode of action (MOA) of toxicity of PA in snowshoe hares. We tested this hypothesis by measuring the effect of PA on the activity of SDH in liver mitochondria isolated from wild hares. In addition, we used molecular modeling to determine the specific binding site of PA on SDH. We found that PA inhibits SDH from hares by an uncompetitive mechanism in a dose-dependent manner. Molecular modeling suggests that inhibition of SDH is a result of binding of PA at the ubiquinone binding sites in complex II. Our results provide a MOA for toxicity that may be responsible for the concentration-dependent anti-feedant effects of PA. We propose that snowshoe hares reduce the dose-dependent toxic consequences of PA by relying on efflux transporters and metabolizing enzymes that lower systemic exposure to dietary PA.

Advancements in Molecular Dynamics Simulations of Biomolecules on Graphical Processing Units
Xu D, Williamson MJ, Walker RC, Annu. Rep. Comp. Chem. 2010, 6:113-324.  DOI: 10.1016/S1574-1400(10)06001-9

Over the past few years competition within the computer game market coupled with the emergence of application programming interfaces to support general purpose computation on graphics processing units (GPUs) has led to an explosion in the use of GPUs for acceleration of scientific applications. Here we explore the use of GPUs within the context of condensed phase molecular dynamics (MD) simulations. We discuss the algorithmic differences that the GPU architecture imposes on MD codes, an overview of the challenges involved in using GPUs for MD, followed by a critical survey of contemporary MD simulation packages that are attempting to utilize GPUs. Finally we discuss the possible outlook for this field.

Symmetry Breaking in Linear ZnCl2+: A Theoretical Study  DOI:10.1016/j.molstruc.2010.03.082
Zou W, Xu D, Zajac P, Cooksy AL, Bersuker IB, Liu Y, Boggs JE, J. Mol. Struct. 2010, 978:263-268.

The four lowest excited states are calculated using the EOMIP-CCSD method and the spectroscopic constants are reported. It is found that all the states have linear structures, while some states are noncentrosymmetric with two equivalent global minima due to the pseudo Jahn–Teller effect. The computed barriers between the two minima are 88 and 169 cm−1, respectively, being much smaller and more reliable than the previous published multi-configurational results because of inaccurate treatment of the static correlation in the latter. For the ground state, anharmonic vibrational levels are also simulated using the finite element method. It shows that the configuration of the vibronic ground state is also noncentrosymmetric.

Distinct Glycan Topology for Avian and Human Sialopentasaccharide Receptor Analogues upon Binding Different Hemagglutinins: A Molecular Dynamics Perspective DOI:10.1016/j.jmb.2009.01.040
Xu D, Newhouse EI, Amaro RE, Pao HC, Cheng LS, Markwick PRL, McCammon JA, Li WW, Arzberger PW. J. Mol. Biol., 2009, 387:465-491.

Hemagglutinin (HA) binds to sialylated glycans exposed on the host cell surface in the initial stage of avi an influenza virus infection. It has been previously hypothesized that glycan topology plays a critical role in the human a daptation of avian flu virus es, such as the potentially pandemic H5N1. Comparative molecular dynamics studies are  complementary to experimental techniques, including glycan microarray, to understand the mechanism of species-specificity switch better. The examined systems comprise explicitly solv ated trimeric form s of avian H3, H5, and swine H9 in complex with avian and human glycan receptor analogues-LSTa (α-2,3-linked lactoseries tetrasaccharide a) and LSTc (α-2,6-linked lactoseries tetrasaccharide c), respectively. The glycans adopted distinct topological profiles with inducible torsional angles when bound to  different HAs. The corresponding receptor binding domain amino acid contact profiles were also distinct. Avian H5 was able to accommodate LSTc in a tightly “folded umbrella”-like topology through interactions with all five sugar residues. After considering conformational entropy, the relative binding free-energy changes, calculated using the molecular mechanics-generalized Born surface area technique, were in agreement with previous experimental findings and provided insights on electrostatic, van der Waals, desolvation, and entropic contributions to HA-glycan interactions. The topology profile and the relative abundance of free glycan receptors may influence receptor binding kinetics. Glycan composition and topological changes upon binding different HAs may be important determinants in species-specificity switch.
Mechanism of Glycan Receptor Recognition and Specificity Switch for Avian, Swine, and Human Adapted Influenza Virus Hemagglutinins: A Molecular Dynamics Perspective  DOI:10.1021/ja904052q
Newhouse EI, Xu D, Amaro RE, Pao HC, Markwick PRL, Wu KJ, Alam M, McCammon JA, Li WW, Arzberger PW. J. Am. Chem. Soc., 2009, 131 (47):17430–17442.

Hemagglutinins (HA’s) from duck, swine, and human influenza viruses have previously been shown to prefer avian and human glycan receptor analogues with distinct topological profiles, pentasaccharides LSTa (α-2,3 linkage) and LSTc (α-2,6 linkage), in comparative molecular dynamics studies. On the basis of detailed analyses of the dynamic motions of the receptor binding domains (RBDs) and interaction energy profiles with individual glycan residues, we have identified ~30 residue positions in the RBD that present distinct profiles with the receptor analogues. Glycan binding constrained the conformational space sampling by the HA. Electrostatic steering appeared to play a key role in glycan binding specificity. The complex dynamic behaviors of the major SSE and trimeric interfaces with or without bound glycans suggested that networks of interactions might account for species specificity in these low affinity and high avidity (multivalent) interactions between different HA and glycans. Contact frequency, energetic decomposition, and H-bond analyses revealed species-specific differences in HA−glycan interaction profiles, not readily discernible from crystal structures alone. Interaction energy profiles indicated that mutation events at the set of residues such as 145, 156, 158, and 222 would favor human or avian receptor analogues, often through interactions with distal asialo-residues. These results correlate well with existing experimental evidence, and suggest new opportunities for simulation-based vaccine and drug development.

Characterizing Loop Dynamics and Ligand Recognition in Human- and Avian-Type Influenza Neuraminidases via Generalized Born Molecular Dynamics and End-Point Free Energy Calculations DOI:10.1021/ja8085643
Amaro RE, Cheng X, Ivanov I, Xu D, Mccammon JA. J. Am. Chem. Soc., 2009, 131 (13):4702–4709.

The comparative dynamics and inhibitor binding free energies of group-1 and group-2 pathogenic influenza A subtype neuraminidase (NA) enzymes are of fundamental biological interest and relevant to structure-based drug design studies for antiviral compounds. In this work, we present seven generalized Born molecular dynamics simulations of avian (N1)- and human (N9)-type NAs in order to probe the comparative flexibility of the two subtypes, both with and without the inhibitor oseltamivir bound. The enhanced sampling obtained through the implicit solvent treatment suggests several provocative insights into the dynamics of the two subtypes, including that the group-2 enzymes may exhibit similar motion in the 430-binding site regions but different 150-loop motion. End-point free energy calculations elucidate the contributions to inhibitor binding free energies and suggest that entropic considerations cannot be neglected when comparing across the subtypes. We anticipate the findings presented here will have broad implications for the development of novel antiviral compounds against both seasonal and pandemic influenza strains.

Solving the Vibrational Schrodinger Equation on an Arbitrary Multidimensional Potential Energy Surface by the Finite Element Method DOI:10.1016/j.cpc.2009.06.010
Xu D
, Stare J, Cooksy AL. Comp. Phys. Comm., 2009, 180 (11):2079-209.

A computational protocol has been developed to solve the bounded vibrational Schrödinger equation for up to three coupled coordinates on any given effective potential energy surface (PES). The dynamic Wilson G-matrix is evaluated from the discrete PES calculations, allowing the PES to be parametrized in terms of any complete, minimal set of coordinates, whether orthogonal or non-orthogonal. The partial differential equation is solved using the finite element method (FEM), to take advantage of its localized basis set structure and intrinsic scalability to multiple dimensions. A mixed programming paradigm takes advantage of existing libraries for constructing the FEM basis and carrying out the linear algebra. Results are presented from a series of calculations confirming the flexibility, accuracy, and efficiency of the protocol, including tests on FHF−, picolinic acid N-oxide, trans-stilbene, a generalized proton transfer system, and selected model systems.

Ensemble-Based Virtual Screening Reveals Potential Novel Antiviral Compounds for Avian Influenza Neuraminidase
Cheng LS, Amaro RE, Xu D, Li WW, Arzberger PW, McCammon JA. J. Med. Chem., 2008, 51 (13): 3878–3894. DOI:10.1021/jm8001197

Avian influenza virus subtype H5N1 is a potential pandemic threat with human-adapted strains resistant to antiviral drugs. Although virtual screening (VS) against a crystal or relaxed receptor structure is an established method to identify potential inhibitors, the more dynamic changes within binding sites are neglected. To accommodate full receptor flexibility, we use AutoDock4 to screen the NCI diversity set against representative receptor ensembles extracted from explicitly solvated molecular dynamics simulations of the neuraminidase system. The top hits are redocked to the entire nonredundant receptor ensemble and rescored using the relaxed complex scheme (RCS). Of the 27 top hits reported, half ranked very poorly if only crystal structures are used. These compounds target the catalytic cavity as well as the newly identified 150- and 430-cavities, which exhibit dynamic properties in electrostatic surface and geometric shape. This ensemble-based VS and RCS approach may offer improvement over existing strategies for structure-based drug discovery.

Statistical Cluster Analysis of Pharmaceutical Solvents  
Xu D, Redman-Furey N. Intl. J. Pharm., 2007, 339:175-188.

High efficiency in polymorph screening and crystallization optimization can be gained by judicious selection of solvents for the study design. Examination of all 57 (classes 2 and 3) pharmaceutical solvents may enable a complete study design but is costly in terms of time and resources. Based on a 17 descriptor dataset specifically constructed for all the classes 2 and 3 pharmaceutical solvents recognized by the International Conference of Harmonization (ICH), an optimal two-stage cluster analysis was carried out together with principal component analysis as a dimensionality and colinearity reduction pre-processor. Both hierarchical average linkage cluster analysis and non-hierarchical K-means cluster analysis converged on a 20-cluster solution with strong statistical criteria support and excellent agreement in cluster memberships, which can be reasonably interpreted from a chemical perspective. This 20-cluster solution is offered as an option for design of more efficient solid state screening studies. Rather than designing a polymorph screen to include all 57 solvents, the inclusion of a single member from each of the 20 clusters would be expected to adequately represent the full range of solvent properties exhibited by the entire 57 member solvent set.

Ab Initio Study of the Torsional Motion in Tolane DOI:10.1016/j.theochem.2007.03.028
Xu D, Cooksy AL. J. Mol. Struct. THEOCHEM, 2007, 815:119-125.

Accurate prediction of the torsional barrier height in tolane is achieved by systematically extrapolating to the Dunning complete basis set limit at the MP2 level with a spin-component-scaled correction. The zero-point energy correction is calculated at the B98/cc-pVDZ level based on a benchmark test using the experimental data for benzene. The final calculated barrier height is 202 cm−1, in agreement with the observed value. The correct barrier height enables the vibrational energy levels and other spectroscopic properties to be determined accurately through numerical integration of the torsional Schrödinger equation. This study provides a nearly complete computational solution to the torsional problem in tolane and may aid the exploration of torsional motions in similar molecules.

Alkylation of Phenol with Tert-butyl Alcohol Catalyzed by Large Pore Zeolites DOI:10.1016/S0926-860X(00)00663-38
Zhang K, Zhang H, Xua G, Xiang S, Xu D, Liu S, Li H. Applied Catalysis A: General, 2001, 207:183-190.

The tert-butylation of phenol was investigated over various zeolite catalysts using tert-butyl alcohol as alkylating agent in a down-flow tubular reactor at atmospheric pressure. Zeolite HY was beneficial to the reaction. The important variables affecting the activity and selectivity of zeolite HY, such as reaction temperature, space velocity and molar ratios of tert-butyl alcohol to phenol, were studied. Zeolite HY hydrothermally treated at high temperatures (above 873 K) was unfavorable to the reaction. The activity and selectivity of zeolite HY were not sensitive to its crystallite diameter change. In the alkylation of phenol with tert-butyl alcohol over zeolite HY catalyst, the suitable reaction temperature range was from 398 to 438 K. Lower reactant molar ratios were beneficial to p-TBP and o-TBP, while higher ones were helpful to produce 2,4-DTBP. Lower space velocities (WHSV (h−1)), i.e. <1.66 (based on phenol) were beneficial to the reaction. The present study shows that zeolite HY has a potential application in the production of tert-butyl phenols with high activity and 2,4-DTBP selectivity.

Alkylation of Phenol with Tert-butyl Alcohol Catalysed by Zeolite Hβ  DOI:10.1016/S0926-860X(97)00376-1
Zhang K, Xu D, Huang C, Zhang H, Xiang S, Liu S, Li H. Applied Catalysis A: General, 1998, 166:89-95.

The catalytic properties of zeolite beta in the tert-butylation of phenol are reported. The influence of various reaction parameters such as temperature, space velocity, molar ratio of the reactants are discussed. Medium acid sites on zeolite Hβ are advantageous in producing p-TBP, and that of strong acid sites are helpful for the formation of 2,4-DTBP, while weak acid sites are effective in producing o-TBP. In order to enhance the selectivity of p-TBP, a proper reaction temperature (418 K), lower reactant ratio (molar ratio) and moderate acidities on zeolite Hβ are recommended. On the basis of the behaviour obtained in the cases mentioned, information about the reaction path models is provided.