# COMPUTATIONAL AND EXPERIMENTAL APPROACHES IN MULTI-TARGET PHARMACOLOGY

EDITED BY: Thomas J. Anastasio PUBLISHED IN: Frontiers in Pharmacology

#### *Frontiers Copyright Statement*

*© Copyright 2007-2017 Frontiers Media SA. All rights reserved. All content included on this site, such as text, graphics, logos, button icons, images, video/audio clips, downloads, data compilations and software, is the property of or is licensed to Frontiers Media SA ("Frontiers") or its licensees and/or subcontractors. The copyright in the text of individual articles is the property of their respective authors, subject to a license granted to Frontiers.*

*The compilation of articles constituting this e-book, wherever published, as well as the compilation of all other content on this site, is the exclusive property of Frontiers. For the conditions for downloading and copying of e-books from Frontiers' website, please see the Terms for Website Use. If purchasing Frontiers e-books from other websites or sources, the conditions of the website concerned apply.*

*Images and graphics not forming part of user-contributed materials may not be downloaded or copied without permission.*

*Individual articles may be downloaded and reproduced in accordance with the principles of the CC-BY licence subject to any copyright or other notices. They may not be re-sold as an e-book.*

*As author or other contributor you grant a CC-BY licence to others to reproduce your articles, including any graphics and third-party materials supplied by you, in accordance with the Conditions for Website Use and subject to any copyright notices which you include in connection with your articles and materials.*

> *All copyright, and all rights therein, are protected by national and international copyright laws.*

*The above represents a summary only. For the full conditions see the Conditions for Authors and the Conditions for Website Use.*

ISSN 1664-8714 ISBN 978-2-88945-252-1 DOI 10.3389/978-2-88945-252-1

### About Frontiers

Frontiers is more than just an open-access publisher of scholarly articles: it is a pioneering approach to the world of academia, radically improving the way scholarly research is managed. The grand vision of Frontiers is a world where all people have an equal opportunity to seek, share and generate knowledge. Frontiers provides immediate and permanent online open access to all its publications, but this alone is not enough to realize our grand goals.

### Frontiers Journal Series

The Frontiers Journal Series is a multi-tier and interdisciplinary set of open-access, online journals, promising a paradigm shift from the current review, selection and dissemination processes in academic publishing. All Frontiers journals are driven by researchers for researchers; therefore, they constitute a service to the scholarly community. At the same time, the Frontiers Journal Series operates on a revolutionary invention, the tiered publishing system, initially addressing specific communities of scholars, and gradually climbing up to broader public understanding, thus serving the interests of the lay society, too.

### Dedication to Quality

Each Frontiers article is a landmark of the highest quality, thanks to genuinely collaborative interactions between authors and review editors, who include some of the world's best academicians. Research must be certified by peers before entering a stream of knowledge that may eventually reach the public - and shape society; therefore, Frontiers only applies the most rigorous and unbiased reviews.

Frontiers revolutionizes research publishing by freely delivering the most outstanding research, evaluated with no bias from both the academic and social point of view. By applying the most advanced information technologies, Frontiers is catapulting scholarly publishing into a new generation.

### What are Frontiers Research Topics?

Frontiers Research Topics are very popular trademarks of the Frontiers Journals Series: they are collections of at least ten articles, all centered on a particular subject. With their unique mix of varied contributions from Original Research to Review Articles, Frontiers Research Topics unify the most influential researchers, the latest key findings and historical advances in a hot research area! Find out more on how to host your own Frontiers Research Topic or contribute to one as an author by contacting the Frontiers Editorial Office: researchtopics@frontiersin.org

# **COMPUTATIONAL AND EXPERIMENTAL APPROACHES IN MULTI-TARGET PHARMACOLOGY**

### Topic Editor:

**Thomas J. Anastasio,** University of Illinois at Urbana-Champaign, United States

Multi-target/multidrug therapeutics: the next frontier in pharmacology. Cover image by Thomas J. Anastasio

The next frontier in pharmacology is the development of Multi-target strategies in which pathological processes are controlled by pharmacologically manipulating them at many different points at once. Designing Multi-target strategies will require deep understanding of the complex physiology that underlies pathological processes. It will also require the development of single drugs with multiple targets, or combinations of drugs with compatible pharmacokinetics that work synergistically to maximize desirable effects while minimizing unwanted side effects. This e-Book contains ten original articles, each addressing a different aspect of this challenge. Together they open new perspectives and show the way forward in the development of Multi-target therapeutics.

**Citation:** Anastasio, T. J., ed. (2017). Computational and Experimental Approaches in Multi-target Pharmacology. Lausanne: Frontiers Media. doi: 10.3389/978-2-88945-252-1

# Table of Contents

### **1. Multi-target Treatments for Multifactorial Diseases**

*05 Editorial: Computational and Experimental Approaches in Multi-target Pharmacology*

Thomas J. Anastasio

### **2. From Medicinal Plants to Multidrug Strategies**


Li Gao, Wei-Feng Wu, Lei Dong, Gui-Ling Ren, Hai-Di Li, Qin Yang, Xiao-Feng Li, Tao Xu, Zeng Li, Bao-Ming Wu, Tao-Tao Ma, Cheng Huang, Yan Huang, Lei Zhang, Xiongwen Lv, Jun Li and Xiao-Ming Meng

### **3. Drug Combination Identification using Computational Brain Models**


Samuel A. Neymotin, Salvador Dura-Bernal, Peter Lakatos, Terence D. Sanger and William W. Lytton

### **4. The Benefits and Challenges of Multi-target Pharmacology**


*101 Developing Multi-target Therapeutics to Fine-Tune the Evolutionary Dynamics of the Cancer Ecosystem*

Lei Xie and Philip E. Bourne


# Editorial: Computational and Experimental Approaches in Multi-target Pharmacology

Thomas J. Anastasio\*

*Computational Neurobiology Laboratory, Department of Molecular and Integrative Physiology, Beckman Institute for Advanced Science and Technology, University of Illinois at Urbana-Champaign, Urbana, IL, United States*

Keywords: polypharmacy, drug combination, drug repurposing, synergy, systems biology, multifactorial process, computer modeling, high throughput screening

**Editorial on the Research Topic**

**Computational and Experimental Approaches in Multi-target Pharmacology**

### MULTI-TARGET TREATMENTS FOR MULTIFACTORIAL DISEASES

Picture yourself in the cockpit of the new BoeingTM 737 MAX airliner, or at the control console of a new American AtomicsTM nuclear reactor. You are in charge, and hundreds to thousands of lives depend on your skillful control of a very complex man-made system. Fortunately, these systems are highly automated, so you need do little more than watch a few displays. Then the airliner goes into a nosedive or the reactor overheats, and the computer fails! You need to take manual control to avoid disaster. To make matters more interesting, imagine that, in order to control that nosediving airliner or that overheating reactor, you have access not to all the controls, or even to several controls, but to only one control. Can you further image that you would succeed in averting disaster?

### Biomedical researchers of many stripes are engaged in battles against multifactorial disease processes that are fought within the dense jungles of very complex physiological systems. Most of them still seem to imagine that they will win the battle by using a single drug to alter the biological properties of a single drug target. How is that working out for them? Take Alzheimer Disease as an example. For decades the Alzheimer field has focused on a single peptide, the amyloid-β peptide, and has devoted vast resources to lowering it using drugs targeting its synthetic enzymes (Armstrong, 2014; Hardy et al., 2014). After all this effort we still lack effective means to halt the neurodegenerative processes associated with Alzheimer Disease. We can't even slow them down.

Increasingly, forward-thinking researchers are calling for the development of multitarget/multidrug treatments for Alzheimer Disease (Bajda et al., 2011; Leon and Marco-Contelles, 2011; Carmo Carreiras et al., 2013). I had my epiphany while creating a computational model of the metabolism of amyloid-β. When I read the literature on the effects of estrogen on this process, in order to connect estrogen with the other elements of my model, I found that this hormone targets not one but at least 10 different elements of the system that regulates amyloidβ (Anastasio, 2013). Hormones, naturally occurring interventional agents that have evolved over eons, achieve control of complex physiological systems by manipulating many system elements simultaneously. We should strive to do the same in identifying treatments for Alzheimer Disease and other multifactorial disorders.

Diseases having multifactorial etiologies include Alzheimer and other neurodegenerative diseases, cancer and cardiovascular disease, diabetes and obesity, and depression and schizophrenia. Multi-target treatments for some multifactorial diseases already exist, and multidrug

Edited and reviewed by: *Salvatore Salomone, University of Catania, Italy*

> \*Correspondence: *Thomas J. Anastasio tja@illinois.edu*

#### Specialty section:

*This article was submitted to Experimental Pharmacology and Drug Discovery, a section of the journal Frontiers in Pharmacology*

> Received: *14 June 2017* Accepted: *20 June 2017* Published: *30 June 2017*

#### Citation:

*Anastasio TJ (2017) Editorial: Computational and Experimental Approaches in Multi-target Pharmacology. Front. Pharmacol. 8:443. doi: 10.3389/fphar.2017.00443* regimens for AIDS, infection by drug-resistant bacteria, cancer, diabetes, and even some mood disorders are by now standard. And the hunt is on for new multi-target approaches. It is widely acknowledged that the main impediment to the design of multidrug/multi-target treatments is the failure to understand the multifactorial processes themselves. New computational models are needed that can represent the interactions among the many factors involved, and new experimental methods are needed to evaluate the validity of the models. Several recent surveys describe the current landscape (Keith et al., 2005; Boran and Iyengar, 2010; Xie et al., 2012; Reddy and Zhang, 2013; Billur Engin et al., 2014; Bulusu et al., 2016). In this Research Topic, leading experts in the area of multi-target pharmacology present their most recent new findings, new models, and new ideas, and show the way forward in the identification of new multi-target/multidrug treatments for multifactorial diseases.

### FROM MEDICINAL PLANTS TO MULTIDRUG STRATEGIES

Medicinal plants are the original multidrug medicines, and many traditional treatments involve plants that have verifiable medicinal properties. For example, Borreria verticillata has been used traditionally in Brazil to treat pain. Silva et al. demonstrate that crude extracts of this plant do indeed have antinociceptive properties, and proceed to analyze its constituents experimentally and computationally.

Medicinal plants were discovered by trail-and-error but multi-target/multidrug therapies could be designed de novo. An example of a designer drug pair is the "binary weapon" of Grixti et al. in which the tumor cell toxicity of one compound is increased through downregulation of its efflux transporter by another compound. The Kell lab provides evidence that various small molecule drugs can increase the toxicity to pancreatic cancer cells of the nucleoside analog gemcitabine. In a study that unifies the traditional and the modern, Gao et al. show how protocatechuic aldehyde, a compound isolated from the Lamiaceae root used in traditional Chinese medicine, can ameliorate some of the serious adverse side effects of the chemotherapeutic agent cisplatin.

### DRUG COMBINATION IDENTIFICATION USING COMPUTATIONAL BRAIN MODELS

Neurological and psychiatric disorders exemplify the challenge of understanding a pathophysiological process well enough to identify an effective polypharmacological treatment for it. Increasingly, computational models are being used to aid the design of effective drug combinations for the treatment of brain diseases. Geerts et al. have developed a computational model of cerebral cortex, featuring a network of many biologically realistic pyramidal neurons and interneurons. Using computational analogs of the working memory tasks that are used to assess cognitive impairment in schizophrenics, they perform in silico screens to predict novel drug combinations that would be effective in ameliorating schizophrenic symptomatology. In a similar vein, Neymotin et al. present a computational model of dystonia, a movement disorder associated with involuntary muscle contractions involving several interacting brain regions. They produced a computational model of these brain regions containing a multitude of biologically realistic model neurons, and use it to suggest new multidrug treatments.

### THE BENEFITS AND CHALLENGES OF MULTI-TARGET PHARMACOLOGY

Perhaps the most obvious way to strike multiple pharmacological targets is to administer multiple drugs, but major challenges in the design of multidrug treatments are mismatches in the pharmacokinetics of the different drugs in the combination. This issue is obviated using single compounds that can strike multiple targets, but finding or synthesizing such multi-target ligands pose challenges of their own. Talevi gives numerous examples of effective multi-target drugs and suggests new ways to identify more. Rastelli and Pinzi elaborate on the multi-target ligand theme and provide an overview of computational tools and related approaches for identification of promising candidate compounds.

Physiological processes are difficult to control not only because they are complex but because they adapt. Xie and Bourne lay out the challenges associated with the development of multitarget strategies to prevent tumor growth due to the resistance to anti-cancer drugs that tumors often develop.

The hoped for response to any drug combination is a synergistic interaction that enhances the desired effects of the individual drugs, or that causes new desired effects to emerge. But synergy in the biological context can occur in various ways and quantifying it is not always straightforward. Tang et al. outline the problems and suggest that the best way to describe synergy is to combine two well known methods. One possible benefit of a multidrug combination is reduction in individual drug dosage such that the desired effect arises synergistically from the combination while unwanted side effects due to individual drugs are minimized. The flip side is the potential drawback that unwanted side effects could be exacerbated, or new side effects could emerge from the combination. The ability to predict the possible side effects of novel compounds would be of value in the design of multidrug strategies, and Lopes et al. describe a new method for doing that.

From drug-resistant bacterial infections to neurodegeneration, the biomedical community faces treatment challenges that involve confronting, understanding, and ultimately manipulating disease processes of great complexity. The articles in this Research Topic direct us along many computational and experimental avenues that we can pursue in identifying multi-target/multidrug treatments for multifactorial disorders.

### AUTHOR CONTRIBUTIONS

TA served as Topic Editor for this Research Topic and also wrote this Editorial.

### ACKNOWLEDGMENTS

I gratefully acknowledge the contributions of the authors of the articles that appear under this Research Topic, and I

### REFERENCES


greatly appreciate the expert and, in many cases, highly engaged participation of the reviewers of the manuscripts prior to publication. I also express my thanks for all the professional help and advice I received from the Frontiers in Pharmacology staff.


**Conflict of Interest Statement:** The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2017 Anastasio. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Antinociceptive Activity of *Borreria verticillata*: *In vivo* and *In silico* Studies

Rosa H. M. Silva<sup>1</sup> \*, Nathália de Fátima M. Lima<sup>1</sup> , Alberto J. O. Lopes <sup>1</sup> , Cleydlenne C. Vasconcelos <sup>1</sup> , José W. C. de Mesquita<sup>2</sup> , Ludmilla S. S. de Mesquita<sup>2</sup> , Fernando C. V. M. Lima<sup>1</sup> , Maria N. de S. Ribeiro<sup>2</sup> , Ricardo M. Ramos <sup>3</sup> , Maria do Socorro de S. Cartágenes <sup>1</sup> \* and João B. S. Garcia<sup>4</sup> \*

<sup>1</sup> Experimental Study of Pain Laboratory, Department of Physiological Sciences, Federal University of Maranhão, São Luís, Brazil, <sup>2</sup> Laboratory of Pharmacognosy, Department of Pharmacy, Federal University of Maranhão, São Luís, Brazil, <sup>3</sup> Research Laboratory Information Systems, Department of Information, Environment, Health and Food Production, Federal Institute of Piauí, Teresina, Brazil, <sup>4</sup> Experimental Study of Pain Laboratory, Department of Pain and Palliative Care, Federal University of Maranhão, São Luís, Brazil

### *Edited by:*

Thomas J. Anastasio, University of Illinois at Urbana–Champaign, United States

#### *Reviewed by:*

Haroon Khan, Abdul Wali Khan University Mardan, Pakistan Birendra N. Mallick, Jawaharlal Nehru University, India Paola Patrignani, University of Chieti-Pescara, Italy Heinrich Korner, University of Tasmania, Australia

#### *\*Correspondence:*

Rosa H. M. Silva Lenna1911@hotmail.com Maria do Socorro de S. Cartágenes scartagenes@gmail.com João Batista Santos Garcia jbgarcia@uol.com.br

#### *Specialty section:*

This article was submitted to Experimental Pharmacology and Drug Discovery, a section of the journal Frontiers in Pharmacology

> *Received:* 13 December 2016 *Accepted:* 04 May 2017 *Published:* 22 May 2017

#### *Citation:*

Silva RHM, Lima NdFM, Lopes AJO, Vasconcelos CC, Mesquita JWCd, Mesquita LSSd, Lima FCVM, Ribeiro MNdS, Ramos RM, Cartágenes MdSdS and Garcia JBS (2017) Antinociceptive Activity of Borreria verticillata: In vivo and In silico Studies. Front. Pharmacol. 8:283. doi: 10.3389/fphar.2017.00283 Borreria verticillata (L.) G. Mey. known vassourinha has antibacterial, antimalarial, hepatoprotective, antioxidative, analgesic, and anti-inflammatory, however, its antinociceptive action requires further studies. Aim of the study evaluated the antinociceptive activity of B. verticillata hydroalcoholic extract (EHBv) and ethyl acetate fraction (FAc) by in vivo and in silico studies. In vivo assessment included the paw edema test, writhing test, formalin test and tail flick test. Wistar rats and Swiss mice were divided into 6 groups and given the following treatments oral: 0.9% NaCl control group (CTRL), 10 mg/kg memantine (MEM), 10 mg/kg indomethacin (INDO), 500 mg/kg EHBv (EHBv 500), 25 mg/kg FAc (FAc 25) and 50 mg/kg FAc (FAc 50). EHBv, FAc 25 and 50 treatments exhibited anti-edematous and peripheral antinociceptive effects. For in silico assessment, compounds identified in FAc were subjected to molecular docking with COX-2, GluN1a and GluN2B. Ursolic acid (UA) was the compound with best affinity parameters (binding energy and inhibition constant) for COX-2, GluN1a, GluN2B, and was selected for further analysis with molecular dynamics (MD) simulations. In MD simulations, UA exhibited highly frequent interactions with residues Arg120 and Glu524 in the COX-2 active site and NMDA, whereby it might prevent COX-2 and NMDA receptor activation. Treatment with UA 10 mg/Kg showed peripheral and central antinociceptive effect. The antinociceptive effect of B. verticillata might be predominantly attributed to peripheral actions, including the participation of anti-inflammatory components. Ursolic acid is the main active component and seems to be a promising source of COX-2 inhibitors and NMDA receptor antagonists.

Keywords: *Borreria verticillata,* COX-2, NMDA receptor, molecular docking, molecular dynamics simulations

## INTRODUCTION

Pain is a warning system that informs the body about the occurrence of tissue damage (Nickel et al., 2012). In the pathophysiology of pain several biological actions are involved, including activation of cyclooxygenase 2 enzyme (COX-2) and N-methyl-D-aspartate (NMDA) receptor.

COX-2 is upregulated in the central nervous system in response to inflammatory factors. It is a rate-limiting enzyme for prostanoid production during inflammation (Ricciotti and Fitzgerald, 2011). Prostaglandin E2, the main pro-inflammatory prostanoid, induces painful hypersensitivity through modulation in the nociceptive pathways, activates the periphery ionic channels such as sodium, calcium, and potentiates the central activation of NMDA and α-amino-3 hydroxy-5-methylsoxazol-4-propionic (AMPA) receptors (Chen et al., 2013). Inhibition of COX-2 enzymatic activity prevents prostanoid production, thus this enzyme is a usual target of non-steroidal anti-inflammatory drugs (NSAIDs) (Zaiss et al., 2014).

The activation of NMDA receptor requires the binding of glycine and glutamate to its subunits GluN1 and GluN2, respectively (Tajima et al., 2016). It is well-known that activation of NMDA receptors causes central sensitization, amplification of spinal nociception, increased ionic conductance and membrane depolarization (Phang and Tan, 2013). For this reason, NMDA receptor antagonists (e.g., memantine) are considered an option in the management of opioid-resistant and chronic pain (Hewitt, 2000).

Therefore, NSAIDs and NMDA receptor antagonists are used to afford pain relief. However, the use of these agents is limited by the occurrence of side effects, such as dizziness, vomiting, constipation, and gastric erosions. These problems and the impact of pain in the quality of life of patients evidence the need of novel therapeutic targets for pain management.

Medicinal plants and their derivatives represent a common alternative for the treatment of diseases (Kandimalla et al., 2016; Zaia et al., 2016). Borreria verticillata (L.) G. Mey., known in Brazil as poaia, cordão-de-frade and vassourinha (Júnior et al., 2012) is traditionally used for various therapeutic purposes including the treatment of pain and inflammatory conditions (Vieira et al., 1999; Souza et al., 2013). It has shown to possess antibacterial (Neto et al., 2002; Ogunwande et al., 2010; Balde et al., 2015), hepatoprotective (Murtala et al., 2015), antioxidant (Abdullahi-Gero et al., 2014a), anti-inflammatory and analgesic (Abdullahi-Gero et al., 2014b) activity.

New technologies have been applied to the assessment of the pharmacological properties of extracts and active principles of medicinal plants, such as molecular docking and molecular dynamic, which is a computer-based approach used to give a prediction of the ligand-receptor complex structure (Meng et al., 2011). The combination of computational technique with biological assay became an important strategy toward finding plant-based drugs (Sharma and Sarkar, 2012).

Considering the factors that contribute to the mechanisms of pain and the use of medicinal plants as multi-targets therapeutic alternatives, the aim of the present study was to assess the antinociceptive activity of the crude hydroalcoholic extract and ethyl acetate fraction of B. verticillata. Furthermore, evaluate the molecular interactions of compounds present in ethyl acetate fraction with COX-2 enzyme and NMDA receptor.

### MATERIALS AND METHODS

### Botanic Material

The aerial parts of Borreria verticillata (L.) G. Mey, Rubiaceae were collected at São José de Ribamar, Maranhão state, (2◦ 33′ 13.3′′ S 44◦ 11′ 22.8′′ W), Brazil, in July 2014. A voucher specimen was deposited at Maranhão Herbarium (MAR), of Federal University of Maranhão (UFMA), under the registration number 5151.

### Obtaining the Hydroalcoholic Extract and the Ethyl Acetate Fraction

Aerial parts of B. verticillata were dried at 38◦C in an oven with circulating air and powdered with a knife mill to obtain a moderately coarse powder (particle sizes under 710 µm and over 250 µm). The powder of B. verticillata aerial parts was macerated with 70% ethanol for 5 days (this step was repeated 3 times) obtaining a solution. The solution was filtered and concentrated to a small volume at 40◦C in a rotary evaporator under vacuum, to obtain the hydroalcoholic extract of B. verticillata (EHBv). EHBv was dissolved in methanol:water (70:30,v/v) for 60 min under mechanical agitation, and successively subjected to liquid-liquid extraction with hexane, chloroform, and ethyl acetate. The solutions were filtered and concentrated at 40◦C in a rotary evaporator under vacuum, to ethyl acetate fraction (FAc).

## Phenolic and Flavonoid Content Assessment

Total phenolic content (TPC) was determined using Folin-Ciocalteu reagent and 20% sodium carbonate. The reaction was kept in the dark for 2 h at room temperature; absorbance was read with a spectrophotometer at 760 nm (Dutra et al., 2014). The PCC was calculated based on the calibration curve plotted with gallic acid standard solutions (1.0–30.0 µg/mL) and is expressed as gallic acid equivalent (mg/mL).

Total flavonoid content (TFC) was determined using a 5% methanol solution of aluminum chloride (AlCl3). The reaction was kept in the dark for 30 min at room temperature; absorbance was read with a spectrophotometer at 425 nm (Dutra et al., 2008). The TFC was calculated based on the calibration curve plotted with quercetin standard solutions (1.0–30.0 µg/mL) and is expressed as quercetin equivalent (mg/mL).

### High-Performance Liquid Chromatography with Ultraviolet-Visible Detector (HPLC UV/Vis)

EHBv and FAc were analyzed with an HPLC device (Thermo Finnigan Surveyor) coupled to an ultraviolet-visible detector and a reversed phase ACE C-18 (250 X 4.6 mm, 5µm) column was used. The components of FAc and EHBv were separated at room temperature through gradient elution at a 1 mL/min flow rate. The mobile phases consisted of purified water with 0.1% acetic acid (A) and acetonitrile (B). The gradient used was as follows: 0–5 min, 20% B; 5–10 min, 25% B; 10–15 min, 25– 23% B; 15–20 min, 23–21% B; 20–25 min, 21–18% B; 25–30 min, 18–15% B; 30–35 min, 15–0% B. The injection volume was 5 µL, and UV-Vis detection was performed at 254 nm. The compounds were identified on the basis of co-injection with standards.

### Gas Chromatography—Mass Spectrometry (GC-MS)

FAc (10 mg) was derivatized in pyridine (300 µL) and bis- (trimethylsilyl) trifluoroacetamide with trimethylchlorosilane (BSTFA/TMCS, 100 µL) and was heated at 80◦C for 1 h. The derivatized product was analyzed with a gas chromatograph (GC-2010, Shimadzu, Japan) coupled to a mass spectrometer (GCMS-QP2010 SE, Shimadzu, Japan) with an Rtx-5MS column (30 m × 0.25 mm ix 0.25 µm, Restek, USA), helium as the carrier gas and a 1.0 mL/min flow rate. The oven temperature was first kept at 70◦C and then set to increase 4◦C/min until 310◦C. The temperature was maintained at 310◦C for 4 min. The injector temperature was set to 250◦C; the injection volume was 1.0 mL at a 1:30 ratio. The mass spectra were obtained by means of electron impact ionization (70 eV) on total ion scanning mode (40 to 1,000 m/z) with the ion source at 200◦C. The compounds were identified through comparison of the obtained mass spectra with the NIST 11 library.

## *In vivo* Biological Studies

### Animals

The present study used adult, male and female Wistar Rattus norvegicus rats with weights ranging from 200 to 300 g and adult, male, and female Mus musculus mice with weights ranging from 25 to 35 g, which were procured from the Central Vivarium (Biotério Central), Federal University of Maranhão (UFMA). Animals were provided with free access to food and water in an environment with controlled temperature and 12/12 h light/dark cycle. This study was carried out in accordance with the recommendations of IASP Guidelines for the Use of Animals in Research. The experimental protocols were approved by the UFMA Ethics in Animal Use Committee (CEUA), ruling no. 17, protocol no. 23115.013545/2015-89.

### Experimental Groups

Six experimental groups with 6 animals each were used. CTRL group was treated oral (p.o) with 0.9% NaCl (0.1 mL/kg); the INDO group was treated with indomethacin (10 mg/kg p.o.); MEM group was treated with memantine (10 mg/kg p.o.); FAc groups were treated with fraction the B. verticillata at doses of 25 mg/kg p.o (FAc 25) and 50 mg/kg p.o (FAc 50) and the EHBv group was treated with the hydroalcoholic extract of B. verticillata (500 mg/kg p.o.). NaCl 0.9% was used as the vehicle to dissolve the solutions.

After the results obtained in the in silico studies, it was observed that of the active compounds ursolic acid (UA) present in the FAc presented better results. Then 6 animals were treated orally with UA (10 mg/kg p.o) and submitted to the carrageenaninduced paw edema test and tail flick.

### Carrageenan-Induced Paw Edema

This test was performed to assess the pharmacological activities of the investigated compounds after subplantar injection of carrageenan. Mice were distributed and treated as described in the "Experimental groups" section. Sixty minutes after the onset of treatments, paw edema was induced through administration of 50 µL of 1% carrageenan via subplantar injection in the right paw; the same volume of 0.9% NaCl was injected in the left paw. The paw volume was measured with a digital plethysmometer 1, 2, 3, 4, and 5 h after induction. Edema was calculated as the difference between the right and left paw volume and is expressed as paw volume variation (ml) over time (Winter et al., 1962; Sadeghi et al., 2011).

### Writhing Test

The acetic acid-induced writhing test is described as a visceralsomatic inflammatory model used for pharmacological screening of central and peripheral antinociceptive activity. Mice were distributed and treated as described in the "Experimental groups" section. Sixty minutes after treatment onset, abdominal writhing was induced through intraperitoneal administration of 0.8% acetic acid (10 mL/kg). The number of contractions was cumulatively counted for 20 min after induction (Koster et al., 1959). The results are expressed as the average number of cumulative abdominal contractions (Shamsi and Keyhanfar, 2013; Mansouri et al., 2014).

### Formalin Test

The formalin test for nociception allows assessment of the neurogenic nociceptive mechanisms triggered by activation of nociceptive fibers and the inflammatory mechanisms activated following the release of inflammatory mediators. Mice were treated as described above; 60 min later, a subplantar injection of 20 µL of 2.5% formalin was administered in the right paw. The nociceptive response, characterized by paw licking or biting, was observed during the first 5 min to assess neurogenic mechanisms and then from minutes 15 to 30 to assess inflammatory mechanisms (Hunskaar and Hole, 1987; Nemoto et al., 2015).

### Tail Flick Test

This test was performed to assess central antinociceptive activity through the stimulation of spinal reflexes. Rats were treated as described above; 60 min later, a thermal stimulus was applied to the final third of the tail (Ugo Basile, Varese-Italy), and the latency to tail flick was measured at baseline, 30, 60, 120, and 180 min. The stimulus intensity was set to obtain 3–6 s latency times; the cutoff point was set to 10 s to avoid tail injury (D'Amour and Smith, 1941; Mansouri et al., 2014).

### *In silico* Studies

### Structure of Compounds and Receptors

The compounds identified in FAc were obtained from the PubChem Project database and were structurally plotted in 3 dimensions (3D) using GaussView 5.0.8 (Dennington et al., 2009). Geometric and vibrational properties were calculated (optimized) under vacuum by means of the density functional theory (DFT) method using functional hybrid B3LYP combined with basis 6–31 ++ G (d, p) in the Gaussian 09 program (Frisch et al., 2009).

The 3D structure of Swiss mouse COX-2 (chain A) was obtained from the Protein Data Bank (PDB, #1DDX). The 3D structures of the drugs MEM and INDO were obtained from the PubChem Project (CID 4054 and 3715, respectively). The structural model of the Rattus norvegicus NMDA receptor subunits GluN1a and GluN2B was obtained by means of homology modeling.

### Homology Modeling

Homology modeling was performed following Ramos et al. (2012) with MODELER 9v14 (Webb and Sali, 2014) and the amino acid sequences of subunits GluN1a and GluN2B (NCBI GI 645985944 and GI 645985945, respectively). As the crystallographic structure of the PDB NMDA (code 4PE5) are not complete, models were generated by homology modeling (HM-GluN1a and HM-GluN2B) using the crystallographic structure of PDB code 4TLL (GluN1/GluN2B NMDA) as template. The quality of the selected models was checked with the programs ProCheck (Laskowski et al., 1993) and Errat (Colovos and Yeates, 1993), run in the SAVES server with Z-Score (ProSA-web Protein Structure Analysis) (Table S2).

### Molecular Docking

The AutoDock 4.2 package (Morris et al., 2008) was used to prepare proteins (refined models) and ligands for docking calculations using the AutoDock Tools (ADT) module, version 1.5.6, according to Ramos et al. (2012). The affinity grid centers were defined on residue Arg120 for COX-2, Tyr513 for NMDA GluN1A and Arg487 for NMDA GluN2B. The initial complex coordinates for MD simulations were selected based on the lowest energy configuration of clusters combined with visual inspection.

### Molecular Dynamics of Complexes

The MD simulations of the complexes selected after molecular docking were performed using GROMACS 4.6.7 software (Berk et al., 2014) following Ramos et al. (2012). The ligand topologies were generated with Automated Topology Builder (ATB) and Repository version 2.1 (Malde et al., 2011). The protonation states of histidine residues were determined using the H++ online server—http://biophysics.cs.vt.edu/hppdetails. php. To enlarge the sample, 3 10-ns MD simulations were performed per complex using different atomic velocities assigned according to the Maxwell distribution. The data generated for the last 4 ns in each simulated system were used for analysis. During the production step, 123 frames were obtained at 100 ps intervals. The details of the interactions were calculated with LigPlot++ software (Laskowski et al., 1993). A minimum of 50% of contact (total of hydrophobic interactions and hydrogen bonds) in the analyzed frames was defined as a criterion of binding efficacy.

### Statistical Analysis

Mean values among experimental groups were compared through univariate analysis of variance (one-way ANOVA) followed by the Newman-Keuls test at p < 0.05. The data were analyzed in the terms of means ± standard errors using GraphPad Prism 5 software.

### RESULTS

### Phenolic and Flavonoid Content Assessment

The concentrations of total phenolic compounds and flavonoids in EHBv were 9.61 and 7.67 mg/mL, respectively.

### HPLC UV/Vis Analysis

Gallic acid, ß-sitosterol, glycyrrhetinic acid, ß-amyrin, caffeic acid, coumaric acid, and quercetin were identified in EHBv (**Table 1** and Figure S1). Gallic acid, ursolic acid, caffeic acid, and ellagic acid were identified in FAc (**Table 1** and Figure S2).

### GC-MS Analysis

Thirteen components were detected in FAc, corresponding to alcohols, sugars, fatty acids, and flavonoids. The main components found were the sugars D-psicofuranose (10.2%) and glucopyranose (33.47%) and the fatty acid 10-undecenoic acid (15.70%) (**Table 2**).

### *In vivo* Biological Studies

### Carrageenan-Induced Paw Edema Test

Subplantar injection of carrageenan induced edema on the animals' paws that lasted the full period of observation, i.e., 5 h, with the peak of edema starting 4 h after induction.

Treatments consisting of EHBv 500 and FAc (FAc 25 and FAc 50) significantly reduced carrageenan-induced edema at 3, 4, and 5 h after induction compared to control. The percentage of reduction in edema caused by the treatment was 41, 42, and 43% for EHBv 500; 48, 61, and 67% for FAc 25; and 41, 53, and 61% for FAc 50 respectively. Treatment with INDO significantly reduced edema by 72, 74, and 77%, while MEM reduced edema by 20, 18, and 14% at 3, 4, and 5 h after induction, respectively, compared to control. Comparisons between the effects of B. verticillata extract, fractions and standard drugs (indomethacin and memantine) revealed that INDO, FAc 25, and FAc 50 were the most efficacious in reducing edema. The effects of FAc at doses of 25 and 50


mg/kg were equivalent, without statistically significant difference (**Figure 1** and Table S1).

### Writhing Test

Intraperitoneal administration of acetic acid caused pain to animals during entire 20 min of assessment. Treatments with EHBv 500, FAc 25, and FAc 50 significantly reduced the number of abdominal contractions by 71, 72, and 42%, respectively, compared to control (**Figure 2**). Treatment INDO and MEM significantly reduced the number of abdominal contractions by 72.5 and 33%, respectively, compared to control. The effects of INDO, EHBv 500 and FAc 25 on writhing test, and the reductions in abdominal contractions induced by FAc 50 and MEM were equivalent.

### Formalin Test

This test assessed pain at 2 different stages. In stage 1 MEM, FAc 25, and FAc 50 treatments significantly reduced neurogenic pain by 54, 45, and 57%, respectively, compared to control, while INDO did not induce a significant reduction of neurogenic

TABLE 2 | Gas chromatography-mass spectrometry (GC-MS) analysis of the ethyl acetate fraction (FAc) from *B. verticillata*.


FIGURE 1 | Paw edema induced by subplantar administration of 1% carrageenan in mice treated orally with NaCl 0.9%, indomethacin 10 mg/kg, memantine 10 mg/Kg, EHBv 500 mg/Kg and FAc (25 and 50 mg/Kg). \*p < 0.05; \*\*p < 0.01; \*\*\*p < 0.001 vs. CTRL (ANOVA; Newman Keuls).

pain. The effects of FAc 50 and MEM were similar. In stage 2, INDO, MEM, EHBv 500, FAc 25, and FAc 50 treatments significantly reduced inflammatory pain by 82, 63, 57, 57, and 59%, respectively, compared to control (**Figure 3**). The effects of EHBv 500, FAc 25 and FAc 50 were similar to MEM and INDO treatments, without statistically significant differences.

### Tail Flick Test

Treatment with MEM significantly reduced pain, as evidenced by 15, 40, and 40% increases in the latency time at 60, 120, and 180 min, respectively, compared to control. INDO, EHBv 500, FAc 25, and FAc 50 were unable to increase the latency time of the animals during the 180 min assessment period (**Figure 4**).

### *In silico* Studies Molecular Docking

ligands are described in **Table 3**.

### All compounds identified in FAc on HPLC UV-Vis analysis were selected for molecular docking. The parameters affinity (binding energy [1Gbind] and inhibition constant [Ki]) and values of all

FIGURE 2 | Writhing induced by the intraperitoneal administration of 0.8% acetic acid (10 ml/kg) in mice treated orally with NaCl 0.9%, indomethacin 10 mg/kg, memantine 10 mg/Kg, EHBv 500 mg/Kg and FAc (25 mg/kg and 50 mg/Kg). \*\*p < 0.01; \*\*\*p < 0.001 vs. CTRL (ANOVA; Newman Keuls).

FIGURE 3 | Formalin test induced by subplantar formalin administration of 2.5% in mice treated orally with NaCl 0.9%, indomethacin 10 mg/kg, memantine 10 mg/Kg, EHBv 500 mg/Kg and FAc (25 and 50 mg/Kg). \*p < 0.05; \*\*p< 0.01; \*\*\*p < 0.001 vs. CTRL (ANOVA; Newman Keuls).

TABLE 3 | Interactions by molecular docking of indomethacin, memantine and compounds identified in FAc with COX-2, GluN1a and GluN2B.


\*∆Gbind, binding energy. \*\*Ki, inhibition constant.

In relation to COX-2, UA exhibited higher 1Gbind values and inhibition constants compared with INDO, which were −9.86 kcal/mol and 0.05 µM and −8.30 kcal/mol and 0.82 µM, respectively. The 1Gbind values corresponding to ellagic acid, caffeic acid, and gallic acid were −7.54, −6.68, and −5.88 kcal/mol, respectively. Relative to the NMDA receptor, UA exhibited the best affinity parameters with GluN1a, with values of −7.02 kcal/mol and 7.13 µM, close to the values obtained for MEM, which were −7.82 kcal/mol and 1.86 µM. With respect to GluN2B, the affinity values of UA were slightly higher compared with MEM, which were −5.69 kcal/mol and 67.35 µM vs. −5.66 kcal/mol and 71.18 µM, respectively. Relative to GluN1a and GluN2B, the 1Gbind values corresponding to ellagic acid, caffeic acid and gallic acid were −5.67, −4.57, and −4.47 kcal/mol and −5.34, −5.16, and −5.04 kcal/mol, respectively.

The interactions of UA, MEM, and INDO with amino acids identified in selected configurations obtained through molecular docking calculations are described in **Table 4**.

### Molecular Dynamics Simulations

The lowest docking-energy conformation of the cluster with lowest energy was chosen as initial structure for the molecular dynamics simulations of the COX-2 (**Figures 5A–C**) GluN1a (**Figure 6A** and Figure S3) and GluN2B (**Figure 6B** and Figure S4) (UA, INDO or MEM) complexes. Interactions ≥50% were considered to be relevant.

In relation to COX-2, UA exhibited high frequencies of interaction with Lys79, Leu80, Lys83, Pro84, Arg120, Leu123, Trp155, Ser199, Phe470, Leu472, Ser481, and Glu524 (**Figure 7A**). The highest frequencies of interaction corresponded to Glu524, with 97% of hydrogen bonds, and Arg120, with 96% of hydrophobic interactions. INDO exhibited high frequencies of interaction with Asn43, Arg44, Thr62, Phe64, Leu80, Trp122, Leu123, lle124, Asp125, Phe470, and Arg469; the most relevant interactions corresponded to Trp122 and Leu123, both with 88% of hydrophobic interactions (**Figure 7B**).

With respect to (**Figure 8A**), UA exhibited high frequencies of interaction with lle497, Pro510, Ala502, Phe507, and Glu506, with the highest frequencies being 82% with Ile497 (67% of hydrophobic interactions and 15% of hydrogen bonds) and 81% with Phe507 (24% of hydrophobic bonds and 57% of hydrogen bonds). Relative to GluN2B (**Figure 8B**), UA exhibited high frequencies of interaction with lle483, Ser488, Phe493, Asp492, and Pro496. The frequencies of interaction with Asp492 and Phe493 were both 66%, corresponding to 2 and 25% of hydrophobic bonds and 64 and 41% of hydrogen bonds, respectively. MEM had maximum contact with GluN1a (**Figure 8C**) and high frequencies of interaction with Phe436, Phe511, Lys512, Ala712, Phe736, Glu764, Met763, and Leu766, with 100% of hydrogen bonds with Lys512 and 84% of interactions with Glu764 (61% of hydrophobic interactions and 23% of hydrogen bonds). Regarding to GluN2B


TABLE 4 | Interactions of ursolic acid, memantine and indomethacin for the conformations chosen by molecular docking.

NA, Not rated.

(**Figure 8D**), MEM exhibited 50% interaction with Glu489 (30% of hydrophobic interactions and 20% of hydrogen bonds).

The frequencies of interaction of UA and MEM with GluN1a and GluN2B are described in **Figure 8**.

### Test In vivo of In silico Selected Compound

After analysis of the results obtained from the in silico tests, tests of Carrageenan-induced paw edema and Tail flick with UA (Sigma-Alderich) were performed.

Treatment with UA (10 mg/kg) significantly reduced carrageenan-induced edema by at 5 h analyzed. These reductions were 57–84% when compared to the control group. Comparisons of the effects UA, INDO, and MEM revealed that UA were most efficacious in reducing edema. Treatment with INDO and of MEM reduced significantly the paw edema in the 3, 4, and 5 h when compared with CTRL. These reductions were similar to those mentioned above (**Figure 9**).

In tail flick test, treatment with UA (10 mg/kg) significantly increased the animal latency time at 30–120 min at 28–42% when compared to control. This effect was not observed after 180 min. Treatment with MEM significantly increased latency time when compared to control. These results were similar to those mentioned above (**Figure 10**).

### DISCUSSION

The assays performed in the present study showed that B. verticillata had peripheral antinociceptive effects and antiinflammatory properties. The in silico tests indicated that UA was chief among the active components identified in FAc, exhibiting relevant interactions with amino-acid residues in COX-2 and NMDA receptor active sites.

The FAc doses (25 and 50 mg/kg) used in this study were chosen according to results obtained in previous studies developed by our research group. Our preliminary studies with EHBv demonstrated that doses 250 and 500 mg/kg had an antinociceptive effect. In order to define the doses used in FAc, we reduced the EHBv dose by 100 times, resulting in doses of 25 and 50 mg/kg.

The anti-edematous effects of EHBv 500, FAc 25, and FAc 50 were similar to INDO starting 3 h after induction. Studies have shown that the paw edema test has a maximum parameter in 3–4 h lasting up to 5 h. After this period of time, induction loses its effectiveness and due to this reason we do not continue the evaluation (Eisenberg et al., 1994; Castardo et al., 2008). Paw edema characteristically develops in 2 stages: the first stage begins with carrageenan administration and lasts up to 2 h, corresponding to the exudative phase of inflammation, which is triggered by the release of inflammatory mediators such as histamine, serotonin and bradykinin. The second stage begins 3 h after carrageenan administration and is characterized by neutrophil infiltration and is sustained by prostanoids (Wilches et al., 2015; Honmore et al., 2016). Along the inflammatory process, or when tissue damage occurs, several substances are released that promote vasodilation, plasma extravasation and cell recruitment. In addition, glutamate release by primary afferent fibers (class Aδ and C) is increased, as is the activation of glutamatergic receptors, such as NMDA, with consequent propagation of action potentials and the release of excitatory neurotransmitters in the posterior horn of the spinal cord (Miller et al., 2011). Accordingly, the anti-edematous effect of B. verticillata FAc might be attributed to COX-2 inhibition and the consequent reduction of PGE<sup>2</sup> synthesis, as it also occurs in the case of oral treatment with INDO. In addition, previous studies showed that blockade of the NMDA receptor also influences inflammation, as it reduces leukocyte migration (Bong et al., 1996) proinflammatory cytokine induction (Morel et al., 2013) and partially inhibits enzyme phospholipase A<sup>2</sup> (Buritova et al., 1996); all of which explain the anti-edematous effect of MEM.

The EHBv 500, FAc 25, and FAc 50 treatments reduced the number of abdominal contractions on the writhing test by 71, 72, and 42%, respectively. The effects of EHBv 500 and FAc 25 were equivalent to INDO and the effect of FAc 50 was equivalent to MEM.

The writhing test is a visceral inflammatory pain model used for screening compounds with peripheral analgesic activity; it involves activation of somatic and visceral receptors in addition to induction of local inflammation mediated by prostaglandins, bradykinin, tumor necrosis factor (TNF) α and interleukins (ILs) 1β and 8 (Rodrigues et al., 2012; Olonode et al., 2015). Some studies have demonstrated increased activation of peripheral receptors and elevated spinal cord glutamate concentrations

in inflammatory pain models (Santos et al., 2005; Shamsi and Keyhanfar, 2013). As the effects of EHBv 500 and FAc 50 treatments were equivalent to INDO, the results of the present study might be mainly attributed to COX-2 inhibition and the consequent reduction of PGE<sup>2</sup> synthesis. In addition, we must also bear in mind the peripheral inhibition of the action of glutamate on nociceptors induced by MEM, as the effects of this drug and FAc 50 were equivalent.

The results of the present study showed that oral treatment with MEM, FAc 25, and FAc 50 reduced the nociceptive response in both phases of the formalin test. Treatments with EHBv and INDO reduced the nociceptive response time only in the inflammatory phase. The reductions in the nociceptive response induced by B. verticillata extract and FAc were similar to MEM and INDO in the first and second phases of formalin test, respectively.

The formalin-induced nociception model involves a biphasic nociceptive response. The first phase (0–5 min) is neurogenic, when class C fibers are predominantly activated; the second phase (15–30 min) corresponds to inflammation and depends

on a combination of inflammatory mediators and peripheral and central sensitization (Hunskaar and Hole, 1987; Tjolsen et al., 1992; Liberato et al., 2003).

Peripheral NMDA receptors contribute to the triggering and maintenance of peripheral sensitization under conditions characterized by cell damage and inflammation (Christoph et al., 2005). Several studies have shown that NMDA antagonists, such as MEM, are able to reduce the nociceptive response during both the first and second phases (Davidson and Carlton, 1998; Liberato et al., 2003). McRoberts et al. (2011) found a 50% reduction of the inflammatory phase in GluN1-knockout rats, which suggests that reduced expression of the NMDA receptor decreases glutamateand substance-P-mediated synaptic signaling in the spinal cord.

NSAIDs such as INDO decrease nociception in the second phase of the nociceptive response only (Tjolsen et al., 1992). On those grounds, one might attribute the reduction of nociception in the first phase to the inhibition of the peripheral action of glutamate on nociceptors and, in the second phase, inhibition of COX-2 and PGE<sup>2</sup> synthesis.

The tail flick test was performed in the present study to investigate the central analgesic potential of EHBv and FAc. The EHBv 500, FAc 25, FAc 50, and INDO treatments were not able to promote central analgesia. In turn, the latency periods were longer (60, 120, and 180 min after treatment, respectively) than with 10 mg/kg MEM. According to Danneman et al. (1994), the tail behavioral response to nociceptive stimuli is predominantly regulated by spinal and supraspinal structures. Thus, the analgesic effect of MEM is due to the blockade of the NMDA receptor, with consequent reduction of central sensitization mediated by excitatory neurotransmitters, such as glutamate (Parsons et al., 1999; Morel et al., 2013).

The absence of the dose-dependent effect can be explained by reducing the dissolution of the higher dose and the fact that probably the largest dose does not have the same amounts of active compounds present at the lowest dose.

On the basis of in vivo tests results, a peripheral antinociceptive activity of B. verticillata can be attributed to a reduction in the activity of COX-2 and the peripheral NMDA receptor. In addition the peripheral antinociceptive activity of B. verticillata detected in the present study agrees with the findings described by Abdullahi-Gero et al. (2014b). In that study, oral and intraperitoneal treatments with the ethanolic extract of B. verticillata leaves in doses of 200–1,000 mg/kg exhibited peripheral and central analgesic and anti-inflammatory effects; furthermore, the results suggested that the 50% lethal dose (LD50) for mice and rats is ≥5,000 mg/kg.

According to the chemical analysis results, the biological activity of B. verticillata might be mainly attributed to phenolic compounds and triterpenes, since constituents of these metabolite classes were identified in EHBv and FAc. Some studies showed that these compounds have anti-inflammatory and analgesic effects, as follows: gallic acid—reductions in allodynia and anti-inflammatory effects (Angélica et al., 2013); caffeic acid—reductions in leukocyte migration and free radical and nitric oxide production (Mehrotra et al., 2011); ellagic acid—interaction with opioid receptors (Mansouri et al., 2014) and UA—inhibition of nuclear factor-kappa B (NF-kB) activity (Takada et al., 2010).

In order to determine the mechanisms and possible multitargets of the peripheral antinociceptive activity of B. verticillata, the compounds identified in the FAc were submitted to in silico studies. Gallic acid, caffeic acid, ellagic acid, and UA were subjected to molecular docking analysis, targeting COX-2, and NMDA receptor.

According to Guimarães et al. (2014), negative 1Gbind values represent favorable interactions of the ligand-receptor complex. The molecular docking interactions in **Table 3** demonstrated that UA had lower binding energies with the NMDA receptor subunits GluN1a (−7.02 kcal/mol), GluN2B (−5.69 kcal/mol) and COX-2 (−9.89 kcal/mol). Therefore, the UA was the ligand that presented better interactions with the NMDA receptor and COX-2, hence it was selected for MD simulations and in vivo tests.

UA interacted with COX-2 and both NMDA receptor subunits as well as the standard drugs INDO and MEM. From a structural point of view, the COX-2 active site consists of a lipophilic channel with a gate formed by residues Arg120, Tyr355, and Glu524 (Rowlinson et al., 2003), and its activation leads to metabolic changes in arachidonic acid (AA) based on interactions with COX-2 residues Arg120, Tyr355, Tyr385, and Ser530, resulting in prostaglandin production (Xu et al., 2014).

Several studies have confirmed the relevance of interactions involving the aforementioned amino acids, thus pointing to them as potential targets in the investigation of COX-2 inhibitors. Fenamic acid derivatives were assessed for COX-2 inhibitory actions in vitro; once the efficacy of the compounds was established, the complexes were subjected to x-ray crystallography. The results indicated interactions of fenamic acid derivatives with Arg120, Tyr355, Tyr385, Trp387, Glu524, and Leu531 (Orlando and Malkowski, 2016). In vitro

FIGURE 7 | Frequency of hydrophobic contacts (green) and hydrogen bonding (orange) of ursolic acid + COX-2 (A) and indomethacin + COX-2 (B). The frequencies correspond to the end 4 of molecular dynamics simulations.

evaluation of the COX-2 inhibitory potential of meloxicam and isoxicam followed by x-ray crystallography detected interactions of both compounds with Met113, Leu117, Arg120, Ile345, Val349, Leu352, Leu359, Phe518, Ala527, Ser530, and Leu531 (Xu et al., 2014). The interaction of COX-2 with ibuprofen, another NSAID with anti-inflammatory and analgesic activities, involved the participation of COX-2 residues Arg120, Val349, Tyr355, Trp387, Met522, Val523, Gly526, Ala527, and Ser530 (Orlando et al., 2015). The results of that study further showed that residues Arg120 and Tyr355, located in the COX-2 channel gate, are crucial for the enzyme's interaction with ibuprofen.

Constant interactions between UA and significant amino acid residues, such as Arg120, Glu524, and others closed related to residues described in previous studies, such as Trp115 (close to Met113 and Leu117, involved in meloxicam and isoxicam interactions), Leu123 (close to Arg120, involved in meloxicam,

carrageenan in mice treated orally with NaCl 0.9%, indomethacin 10 mg/kg, memantine 10 mg/Kg, and ursolic acid 10 mg/Kg. \*p< 0.05; \*\*p< 0.01; \*\*\*p < 0.001 vs. CTRL (ANOVA; Newman Keuls).

isoxicam and ibuprofen interactions), and Glu524 (close to Val523, involved in meloxicam and ibuprofen interactions, and to Gly526, involved in ibuprofen interactions).

The present study further found that the UA and COX-2 interaction profile detected upon molecular docking persisted in the MD simulations, even in the presence of temperature, pressure, water, and protein flexibility variations. In addition, the interactions of UA with amino acid residues adjacent to the COX-2 active site and the AA metabolism site, such as Tyr115, Ser119, Tyr122, Leu123, Asp125, and Pro528, might also intensify the COX-2 inhibitory potential of UA. On the basis of the referred facts: (1) UA might block the access of ligands to the COX-2 active site; (2) interactions with Arg120 and amino acid residues adjacent to the COX-2 active site might reduce arachidonic acid metabolization; and (3) interactions of COX-2 active site inhibitors with Arg120 and Glu524 have already been described, which demonstrates the relevance of these amino acid residues for the possible blocking of the COX-2 active site by UA (Figure S5).

From a structural point of view, the NMDA receptor is a heterotetramer formed by 2 subunits, GluN1 and GluN2. These subunits are modulated by domains ATD, LBD, TMD, and CTD (Karakas and Furukawa, 2014). Activation of the NMDA receptor requires simultaneous binding of glutamate to GluN2 and glycine or d-serine to GluN1 through the interactions of the latter αamino and α-carboxyl groups with LBD regions. The LBD is formed by 3 transmembrane regions (M1, M2, and M3) and 2 segments (S1 and S2), which comprise the active site and are located in subdomains D1 and D2, respectively (Traynelis et al., 2010).

According to Kaye et al. (2006), the active site of subdomain D1 is formed by amino acid residues Gln405, Phe484, Thr486, Pro516, Thr518, and Arg523, and the active site of subdomain D2 is formed by Ser688, Ser687, Val6869, Trp731, Asp732, and Ser756. Several studies with GluN1 antagonists, such as 5,7-dichlorokynurenic acid (DCKA) and 1-thioxo-1,2-dihydro-[1,2,4]triazolo[4,3-a]quinoxalin-4(5H)-one (TK40), verified interactions with Phe484, Thr518, Pro516, and Arg523 (Furukawa and Gouaux, 2003) for DCKA and with Glun405, Phe484, Pro516, Leu517, Thr518, Arg523, Ser687, Ser688, and Val689 for TK40 (Kvist et al., 2013). These interactions correspond to residues Gln383, Phe462, Pro494, Leu495, Thr496, Arg501, Ser665, Ser666, and Val667 in our structure. After 10 ns MD simulations, we found high frequencies of interactions of UA with lle497, Ala502, Phe507, Ser508, Lys509, and Pro510, located close to the sites indicated by previous authors, which points to the potential of UA to inhibit glycine binding to GluN1.

Lee et al. (2014) developed a crystallographic structure of GluN2B complexed with the partial agonist trans-1 aminocyclobuthane-1,3-dicarboxylic acid (t-ACBD). Those authors revelead that t-ACBD interacted with the amino acid residues His479 and Thr507, which are equivalent to His454 and Thr482 in the present study. Addicionally, UA presented interactions with Ile483, Ser48, and all other neighboring amino acid residues, thus demonstrating that UA might prevent binding of agonists to the NMDA receptor GluN2B subunit.

Through in silico tests was demonstrated that UA has the COX-2 enzyme and NMDA receptor as multi-targets. The interactions of UA with these two pharmacological targets suggest that there is inhibition of COX-2 activation, NMDA receptor and consequently reduction of the painful process. Compounds and multi-target drugs are currently a promising alternative for the treatment of complex and multifactorial pathologies. They present advantages over polypharmaceutical therapies such as: patient adherence to treatment, improved kinetics and therapeutic efficacy, increased therapeutic range, reduction of drug interactions, adverse and toxicological effects (Talevi, 2015; Nicolik et al., 2016). Because UA has exhibited this feature (multi-targets), it may be considered an advantage over other compounds described in the literature. Verano et al. (2013) has shown that UA 10 mg/kg (i.p) has an antinociceptive effect, so we chose this dose to test its effect orally and have confirmed the antinociceptive and anti-inflammatory effects of UA through inhibitions of the TRPV1 receptor and serotonergic synergism (5-HT); reduced IL-2, IL-6, interferon (IFN)-γ, TNF-α, reactive oxygen species (ROS), phospholipase A2 and NF-κβ release; reduced COX-2 and inducible nitric oxide synthase (iNOS) expression; (Kashyap et al., 2016) and favorable interactions for complexes formed with COX-1 and COX-2 (Magalhães et al., 2012). Zhang et al. (2013) demonstrated that UA can induce apoptosis of cancer cells by reducing COX-2 expression and Ma et al. (2014) demonstrated that UA reduced COX-2 expression in CCl4-treated animals. The results presented here evidenced that UA exerts peripheral antinociceptive effect by inhibiting COX-2 in carrageenan- induced paw edema test and analgesic central in the tail flick test.

Until now there are no studies performed in vivo or in silico to analyse the antinociceptive potential of compounds present in the aerial parts of B. verticillata. We suggest that the peripheral antinociceptive effect of this plant species is mainly due to actions of anti-inflammatory components. Possibly contributing to reducing the activity of the enzyme COX-2 and peripheral NMDA receptors, with consequent reductions in pain. In silico and in vivo studies allowed the selection and suggestion of UA as the main active compound in B. verticillata, as it exhibited desirable affinity parameters, stable interactions with COX-2 and NMDA receptor subunits GluN1a and GluN2B and presented peripheral and central analgesic effect. Thus, the UA is configured a promising molecule for the development of COX-2 inhibitors and NMDA receptor antagonists.

### AUTHOR CONTRIBUTIONS

Conceived and designed the experiments: RS, MC, JD. Methodology: RS, NL, AL, CV, FL, JD, LD, RR. Analyzed the data: RS, AL, CV, JD, LD, RR, MC,

### REFERENCES


JG. Contributed reagents/materials/analysis tools: RR, MR. Wrote the paper: RS, AL, CV, JD, LD, RS, MC, JD.

### FUNDING

The present study was funded by the Maranhão Research Foundation (FAPEMA) and the National Council for Scientific and Technological Development (CNPq).

### ACKNOWLEDGMENTS

We would like to thank the National High-Performance Processing Center (Centro Nacional de Processamento de Alto Desempenho/CENAPAD-UFC) from Federal University of Ceará (UFC), for its availability to perform the in silico assays and Prof. Dr. Eduardo B. de Almeida Júnior, the botanist responsible for the Maranhão Herbarium (MAR), for the botanical identification of the plant species used in the present study.

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fphar. 2017.00283/full#supplementary-material


NR1 ligand-binding core. EMBO J. 12, 2873–2885. doi: 10.1093/emboj/ cdg303


neurons decreases pain during phase 2 of the formalin test. Neuroscience 172, 474–482. doi: 10.1016/j.neuroscience.2010.10.045


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2017 Silva, Lima, Lopes, Vasconcelos, Mesquita, Mesquita, Lima, Ribeiro, Ramos, Cartágenes and Garcia. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Enhancing Drug Efficacy and Therapeutic Index through Cheminformatics-Based Selection of Small Molecule Binary Weapons That Improve Transporter-Mediated Targeting: A Cytotoxicity System Based on Gemcitabine

#### Justine M. Grixti 1, 2, Steve O'Hagan2, 3, 4, Philip J. Day 1, 2 and Douglas B. Kell 2, 3, 4 \*

#### Edited by:

*Thomas J. Anastasio, University of Illinois at Urbana–Champaign, USA*

#### Reviewed by:

*Peng Hsiao, Seattle Genetics Inc., USA David Dickens, University of Liverpool, UK*

#### \*Correspondence:

*Douglas B. Kell dbk@manchester.ac.uk*

#### Specialty section:

*This article was submitted to Experimental Pharmacology and Drug Discovery, a section of the journal Frontiers in Pharmacology*

> Received: *13 December 2016* Accepted: *10 March 2017* Published: *27 March 2017*

#### Citation:

*Grixti JM, O'Hagan S, Day PJ and Kell DB (2017) Enhancing Drug Efficacy and Therapeutic Index through Cheminformatics-Based Selection of Small Molecule Binary Weapons That Improve Transporter-Mediated Targeting: A Cytotoxicity System Based on Gemcitabine. Front. Pharmacol. 8:155. doi: 10.3389/fphar.2017.00155* *<sup>1</sup> Faculty of Biology, Medicine and Health, University of Manchester, Manchester, UK, <sup>2</sup> Manchester Institute of Biotechnology, University of Manchester, Manchester, UK, <sup>3</sup> School of Chemistry, University of Manchester, Manchester, UK, <sup>4</sup> Centre for Synthetic Biology of Fine and Speciality Chemicals, University of Manchester, Manchester, UK*

The transport of drug molecules is mainly determined by the distribution of influx and efflux transporters for which they are substrates. To enable tissue targeting, we sought to develop the idea that we might affect the transporter-mediated disposition of small-molecule drugs via the addition of a second small molecule that of itself had no inhibitory pharmacological effect but that influenced the expression of transporters for the primary drug. We refer to this as a "binary weapon" strategy. The experimental system tested the ability of a molecule that on its own had no cytotoxic effect to increase the toxicity of the nucleoside analog gemcitabine to Panc1 pancreatic cancer cells. An initial phenotypic screen of a 500-member polar drug (fragment) library yielded three "hits." The structures of 20 of the other 2,000 members of this library suite had a Tanimoto similarity greater than 0.7 to those of the initial hits, and each was itself a hit (the cheminformatics thus providing for a massive enrichment). We chose the top six representatives for further study. They fell into three clusters whose members bore reasonable structural similarities to each other (two were in fact isomers), lending strength to the self-consistency of both our conceptual and experimental strategies. Existing literature had suggested that indole-3-carbinol might play a similar role to that of our fragments, but in our hands it was without effect; nor was it structurally similar to any of our hits. As there was no evidence that the fragments could affect toxicity directly, we looked for effects on transporter transcript levels. In our hands, only the ENT1-3 uptake and ABCC2,3,4,5, and 10 efflux transporters displayed measurable transcripts in Panc1 cultures, along with a ribonucleoside reductase RRM1 known to affect gemcitabine toxicity. Very strikingly, the addition of gemcitabine alone increased the expression of the transcript for ABCC2 (MRP2) by more than 12-fold, and that of RRM1 by more than fourfold, and each of the fragment "hits" served to reverse this. However, an inhibitor of ABCC2 was without significant effect, implying that RRM1 was possibly the more significant player. These effects were somewhat selective for Panc cells. It seems, therefore, that while the effects we measured were here mediated more by efflux than influx transporters, and potentially by other means, the binary weapon idea is hereby fully confirmed: it is indeed possible to find molecules that manipulate the expression of transporters that are involved in the bioactivity of a pharmaceutical drug. This opens up an entirely new area, that of chemical genomics-based drug targeting.

Keywords: binary weapon, cheminformatics, gemcitabine, anticancer drugs, pancreatic cancer, drug transporters, phenotypic screening

## INTRODUCTION

In a typical small molecule drug discovery programme pipeline, candidate ("hit") compounds for treating a particular disease are selected from a large chemical library, and after various modifications (to form "leads" and variants thereof) enter "phase 1," a testing for safety at low doses in healthy volunteers. "Attrition" is a term used to describe the failure of such molecules to progress further to market, via phases 2 and 3 (small and larger clinical trials) (Kola and Landis, 2004; Empfield and Leeson, 2010; Leeson and Empfield, 2010; Leeson, 2016). Nowadays attrition occurs largely for reasons of toxicity or lack of efficacy (Kola and Landis, 2004; Arrowsmith and Miller, 2013), and runs in excess of 90% (e.g., Kola and Landis, 2004; Kell, 2013, and see for full details http://csdd.tufts.edu/files/uploads/Tufts\_CSDD\_ briefing\_on\_RD\_cost\_study\_-\_Nov\_18,\_2014..pdf), with gross pharmacokinetics and pharmacodynamics (as assessed at the whole organ level) being seen as less of an issue than it once was (Kola and Landis, 2004). The simple consequence of this level of attrition is that it costs ∼10 times more than it might, per molecule, now as much as \$2.5 Bn, to bring a drug successfully to market.

## Role of Transporters in Cellular Drug Uptake

We have argued that a lack of understanding of human metabolism and of the transporters necessary to get orally active drugs across intestinal epithelia and into target cells is one of the chief causes of attrition. By now, following a similar programme in yeast (Herrgård et al., 2008), we do have a reasonable model of the human metabolic network (Swainston et al., 2013, 2016; Thiele et al., 2013), with fully one third of the steps involving some kind of transport(er), and with uptake transporters of the SoLute Carrier families (SLCs) (Hediger et al., 2004, 2013) being woefully understudied (César-Razquin et al., 2015). In particular, although it remains underappreciated, we have rehearsed on multiple occasions the abundant evidence that the non-transporter- (i.e., bilayer-) mediated uptake of drugs through intact cell membranes is normally negligible (e.g., Dobson and Kell, 2008; Dobson et al., 2009a,b; Kell et al., 2011, 2013, 2015; Lanthaler et al., 2011; Kell, 2013, 2015a,b, 2016a,b; Kell and Goodacre, 2014; Kell and Oliver, 2014; Mendes et al., 2015; O'Hagan and Kell, 2015a; Kell, 2015a,b), a striking recent example being that of Superti-Furga and colleagues (Winter et al., 2014). This shifts the agenda to one of molecular enzymology and systems biology, in which we need to discover (i) which transporters transport which drugs (Giacomini et al., 2010; Sugiyama and Steffansen, 2013), (ii) their expression profiles in different membranes and tissues, and (iii) their kinetic properties. In other words it leads us to recognize that this is fundamentally a problem of systems pharmacology (e.g., van der Greef and Mcburney, 2005; Berger and Iyengar, 2009; van der Graaf and Benson, 2011; Antman et al., 2012; Rostami-Hodjegan, 2012; Waldman and Terzic, 2012; Zhao and Iyengar, 2012; Kell and Goodacre, 2014; Westerhoff et al., 2015; Kell, 2015a).

### Transporter-Mediated Drug Targeting

A particularly nice example of the overwhelming use of transporters for drug uptake comes from the study of Superti-Furga and colleagues (Winter et al., 2014) using haploid cells and determining that very much less than 1% of sepantronium uptake could have occurred other than via a specific SLC called SLC35F2. In a similar and complementary vein, the expression profile of specific transporters allows one to target drug substrates to the particular tissues in which the relevant transporters are most highly expressed. This has been illustrated beautifully by Pfefferkorn and colleagues for both a glucokinase activator (Pfefferkorn et al., 2012; Pfefferkorn, 2013; Sharma et al., 2015) and a "statin"-type drug (Pfefferkorn et al., 2011) that are both targeted to the liver via proteins of the Organic Anion Transport Protein (OATP/SLCO/SLC21) (Hagenbuch and Stieger, 2013) family. In this case substantial concentration ratios of e.g., hepatocyte: pancreas of 50:1 (Pfefferkorn et al., 2012) and hepatocyte:myocyte of 250,000:1 (Pfefferkorn et al., 2011) could be achieved (a finding hard to explain on the basis of any significant bilayer permeability!). Other examples of tissue-selective drug targeting include a liver-targeted stearoyl desaturase inhibitor (Oballa et al., 2011; Ramtohul et al., 2011; Liu, 2013), various other liver-targeted drugs based on OATPs (Buxhofer-Ausch et al., 2013; Tu et al., 2013), and a prostatespecific targeting of an iodide transporter for radio-iodinemediated cell killing (Kakinuma et al., 2003).

These examples show what can be achieved in terms of drug targeting if the transporter distribution happens to work to one's advantage "naturally," but cannot be exploited directly when it does not.

The glucokinase activator case is important, since if such drug molecules were allowed to enter all tissues they proved toxic (Pfefferkorn et al., 2012; Pfefferkorn, 2013). A similar and particular case of interest is that of broadly cytotoxic anticancer drugs, where we evidently need mechanisms to target them solely to the tissue of interest, and where we might then greatly improve their therapeutic index. Since the tissue-dependent expression of such transporter molecules is highly heterogeneous (see e.g., almost any dataset in the human protein atlas http:// proteinatlas.org/ Uhlén et al., 2015, including those for SLC28 http://www.proteinatlas.org/search/slc28 and SLC29 http://www. proteinatlas.org/search/slc29), it must be subject to regulation (e.g., Pennycooke et al., 2001; Del Santo et al., 2001; Fernández-Veledo et al., 2004, 2007; Plant, 2016). Thus, just as with the small-molecule-driven induction of pluripotent stem cells (Okita et al., 2007; Feng et al., 2009; Desponts and Ding, 2010; Li and Ding, 2010; Zhang, 2010; Grskovic et al., 2011; Li et al., 2012, 2014; Li X. et al., 2015; Jung et al., 2014; Kang et al., 2014), that regulation can similarly be affected by pharmacological intervention with small molecule effectors. Thus, our aim was to seek small molecules that were themselves without cytotoxic effects but that could increase the response of different target cells to anti-cancer drugs that are otherwise present at only a barely cytotoxic level, in particular by modulating the level of activities of specific uptake transporters. It differs from the use of pairs of existing drugs of known activities (e.g., Borisy et al., 2003; Lehár et al., 2007, 2008, 2009a,b; Zimmermann et al., 2007; Wright, 2016), but, interestingly, bears a clear resemblance to the overall strategy used in traditional Chinese medicine where a "shi" ("courier") herb is used to assist the delivery of the main ingredient ("Jun" or "Emperor" herb) to its site of action (Zhao et al., 2015). We refer to this combination as a "binary weapon."

## Gemcitabine and Pancreatic Cancer

The nucleoside analog gemcitabine (2',2'-difluorodeoxycytidine, Gemzar <sup>R</sup> ) (Alvarellos et al., 2014) is one of the most commonly used chemotherapeutic agents in pancreatic adenocarcinoma, the carcinoma with arguably the least favorable prognosis (5 year survival time) of any (Bhattacharjee et al., 2014; Waddell et al., 2015). Like all nucleoside inhibitors of this type, it must first be transported into the cell and then be metabolized (phosphorylated) to exert its clinical action (thereby lowering its ability to act as a substrate for efflux pumps Fukuda and Schuetz, 2012, though see below). Gemcitabine has multiple intracellular targets, and up-regulation of these targets or nucleosidemetabolizing enzymes such as ribonucleotide reductase (RRM1) may confer resistance to this drug (Bergman et al., 2002, 2005; Nakano et al., 2007; Minami et al., 2015). The main uptake transporters are considered to be ENT1 (SLC29A1) and CNT1/3 (SLC28A1/3) of the SLC28/29 families (Kong et al., 2004; Podgorska et al., 2005; Veltkamp et al., 2008; Young et al., 2008, 2013; Molina-Arcas et al., 2009; Cano-Soldado and Pastor-Anglada, 2012; Molina-Arcas and Pastor-Anglada, 2013) (**Table 1**). SLC28 transporters are sodium-dependent concentrative nucleoside transporters (Smith et al., 2007), while SLC29 are equilibrative. Notably, there is considerable evidence that the potency (cytotoxicity) of gemcitabine is strongly related to the expression level(s) of these transporters (e.g., Burke et al., 1998; Mackey et al., 1998a,b; Baldwin et al., 1999; Rauchwerger et al., 2000; Cass, 2001; Achiwa et al., 2004; Spratlin et al., 2004; Giovannetti et al., 2005, 2006, 2007; King et al., 2006; Marcé et al., 2006; Mey et al., 2006; Mini et al., 2006; Leung and Tse, 2007; Mori et al., 2007; Oguri et al., 2007; Zhang et al., 2007; Cano-Soldado et al., 2008; Molina-Arcas et al., 2008; Pérez-Torras et al., 2008; Veltkamp et al., 2008; Andersson et al., 2009; Damaraju et al., 2009; Farrell et al., 2009, 2016; Köse and Schiedel, 2009; Maréchal et al., 2009, 2012; Wong et al., 2009; Hagmann et al., 2010; Lane et al., 2010; Molina-Arcas and Pastor-Anglada, 2010; Okazaki et al., 2010; Paproski et al., 2010, 2013; Santini et al., 2010, 2011; Spratlin and Mackey, 2010; Tanaka et al., 2010; Bhutia et al., 2011; De Pas et al., 2011; Gusella et al., 2011; Komori et al., 2011; Matsumura et al., 2011; Borbath et al., 2012; Choi, 2012; Gesto et al., 2012; Kobayashi et al., 2012; Koczor et al., 2012; Morinaga et al., 2012; Murata et al., 2012; Eto et al., 2013; Nakagawa et al., 2013; Skrypek et al., 2013; Xiao et al., 2013; Chan et al., 2014; Deng et al., 2014; Greenhalf et al., 2014; Khatri et al., 2014; Koay et al., 2014; Lee et al., 2014; Lemstrová et al., 2014; Liu et al., 2014; Nordh et al., 2014; Tavano et al., 2014; Wu et al., 2014; de Sousa Cavalcante and Monteiro, 2014; Hung et al., 2015; Pastor-Anglada and Pérez-Torras, 2015; Yamada et al., 2016).

ABC-type efflux transporters are heavily involved in drug resistance in both mammals (e.g., Liu et al., 2005; Fukuda and Schuetz, 2012; Rosenberg et al., 2015; Silva et al., 2015) and microbes (e.g., Putman et al., 2000; Du et al., 2014; Prasad and Rawal, 2014; Li X.-Z. et al., 2015), and may be of value in industrial biotechnology (Kell et al., 2015). Gemcitabine may also be a substrate for certain efflux transporters such as ABCG2/BRCP (König et al., 2005; Keppler, 2011; Chen et al., 2012; Lemstrová et al., 2014), although "knockdown of ABCC3, ABCC5 or ABCC10 individually did not significantly increase gemcitabine sensitivity" (Rudin et al., 2011). Finally, gemcitabine may also be deaminated in plasma, leading to its clearance (Hodge et al., 2011).

Other small molecules known to affect the response of pancreatic cancer cells to gemcitabine include nicotine (Banerjee et al., 2013, 2014), while molecules that affect nucleoside transporter expression include bile acids (Klein et al., 2009). Finally, erlotinib, gefitinib, and vandetanib inhibit human nucleoside transporters and thereby protect cancer cells from gemcitabine cytotoxicity (Damaraju et al., 2014), while a variety of kinase inhibitors (Huang et al., 2002, 2003, 2004) and dihydropyridine-type calcium channel antagonists (Li et al., 2007) may also affect nucleoside transport.

In particular, however, and not least since the small molecule indole-3-carbinol (which is probably converted to 3,3′ -diindolylmethane Banerjee et al., 2009) had been stated to increase both ENT1 expression and the sensitivity of pancreatic carcinoma cells to gemcitabine (Wang et al., 2011), possibly acting via miRNA-21 (Giovannetti et al., 2010; Hwang et al., 2010; Melkamu et al., 2010; Paik et al., 2013), as too did the molecule "S-1" (Nakahira et al., 2008; Jordheim and Dumontet, 2013), the gemcitabine/nucleoside transporter system seemed ideal for the test of our "binary weapon" strategy. The present paper reports the results of this approach.


1 | The main human nucleoside transporters and some of their

## MATERIALS AND METHODS

### Cells and Reagents

The human pancreatic duct epithelioid carcinoma cell line, Panc1 (see Gou et al., 2007), and the human embryonic kidney cell line, HEK293 were grown in Dulbecco's modified Eagle's medium (DMEM) (Sigma). The human bone marrow neuroblastoma cell line, SH-SY5Y was grown in a 1:1 mixture of Eagle's Minimum Essential Medium (Sigma) and F12 Medium (Sigma). All cell culture media were supplemented with 10% heat-inactivated fetal bovine serum (FBS), 200 mM L-glutamine, and a 5 mL solution containing 10,000 units.mL-1 penicillin and 10 mg.mL-1 streptomycin. The immortal human pancreatic duct epithelial cell line, hPDE was grown in Keratinocyte-SFM (1X) medium (ThermoFisher), supplemented with 10 mg.mL-1 streptomycin. All four cell lines were obtained and karyotyped locally. Cells were routinely maintained at 37◦C in a humidified 5% CO2 atmosphere, in continuous exponential growth at a cell density ranging between 1 × 10<sup>5</sup> and 1 × 10<sup>6</sup> cells.mL-1, by passaging every 3 or 4 days. Cell line authenticity was confirmed through karyotype testing (University of Manchester, UK).

### Cell Growth/Viability Assay

Cells were seeded in a 96-well plate at a density of 5,000 cells/well, in triplicate, and left to attach. Gemcitabine, present at different concentrations, was added directly to the cells, and left to incubate for an additional 96 h. Cells were then subjected to the MTT Cell Proliferation Assay as per the manufacturer's instructions (Sigma). Absorbance at 570 nm was measured 3 h after the addition of 10 µL of MTT salt reagent/well.

### Maybridge Fragment Screening

Maybridge fragments (MBFs) obeying the "rule of three" (Congreve et al., 2003) were supplied at 100 mM in DMSO and were deployed into the assay plates using an ECHO contactless liquid handler (Labcyte, Inc). For screening purposes, the first 500 MBFs (Library 1) were pooled, i.e., each well in a 96-well plate had a pool of six MBFs. Cells were seeded in a 96-well plate at a density of 5,000 cells/well, in triplicate, and left to attach overnight. Following incubation, the growth medium was replaced with fresh medium containing the pooled MBFs, each fragment present at 10 µM, followed by an additional 24 h incubation. The cells were further incubated with the fragments in the presence of gemcitabine at 20 nM for 96 h. Cells were then subjected to the MTT Cell Proliferation Assay as described above.

To study the effect of each MBF on its own rather than in a pool, the candidate pooled fragments (i.e., showing activity) were de-convolved, i.e., one MBF/well, and cells were plated and treated as described above.

### Specificity Experiments

SH-SY5Y cells were seeded at a density of 12,500 cells/well, HEK293 and hPDE cells at a density of 10,000 cells/well, in triplicate in a 96-well plate, and left to attach overnight. Following incubation, the medium was replaced with fresh medium containing the MBF hits (i.e., MBF D1, B1, 10, 11, 12, and 20) at 10 µM, followed by an additional 24 h incubation period. Cells were further incubated with the fragments in the presence of gemcitabine at 100 nM for 72 h (SH-SY5Y cells) and 96 h (HEK293 and hPDE cells). Cell viability was then assessed using the MTT Cell Proliferation Assay.

### Cheminformatic Analyses

These were all performed as in our previous work of this type (O'Hagan and Kell, 2015a,b,c; O'Hagan et al., 2015; O'Hagan and Kell, 2016), using the KNIME workflow system (see e.g., Berthold et al., 2008; Mazanetz et al., 2012; Meinl et al., 2012; Warr, 2012; O'Hagan and Kell, 2015b and http://knime.org/).

### Maybridge Fragment Titration Experiments

Cells were seeded in a 96-well plate at a density of 5,000 cells/well, in triplicate, and left to attach overnight. Following incubation, the medium was replaced with fresh medium containing MBFs at different concentrations (3, 10, 30, 100, and 300 µM) followed by an additional 24 h incubation. Cells were further incubated with the fragments in the presence of gemcitabine at 100 nM for 96 h. Cells were then assessed using the MTT Cell Proliferation Assay as described earlier.

### Cell Culture Treatments for Gene Dysregulation Studies

To examine the effect of gemcitabine and MBFs, alone or in combination; on expression of the influx and efflux transporter genes and of the RRM1 gene, cells were seeded in a 6-well plate at a density of 30,000 cells/well, in duplicate, and left to attach overnight. Following incubation, in studies where the effects of the MBFs alone were studied, the medium was replaced with fresh medium containing MBFs at 10 µM, followed by further incubation for 24 h. For studies where the effects of gemcitabine alone were studied, the medium was replaced with fresh medium containing gemcitabine at 100 nM, followed by further incubation for 96 h. For studies in which cells were treated with gemcitabine in combination with the fragments, the cells were first pre-treated with MBFs at 10 µM for 24 h, followed by further incubation with the fragments at 10 µM in the presence of gemcitabine at 100 nM for 96 h. Cells were harvested using TRIzol <sup>R</sup> reagent (Life Technologies) and stored in −80◦C until use.

### Total RNA Isolation and Quantitative Real-Time Reverse Transcription Polymerase Chain Reaction (RT-qPCR)

Following treatment as described above, total cellular RNA was isolated from the cells using the RNeasy isolation kit (Qiagen) according to the manufacturer's instructions. RNA concentration was determined using a NanoDrop <sup>R</sup> Spectrophotometer (NanoDrop ND-1000, NanoDrop Technologies, Wilmington, USA). The OD\_260/280 nm ratios of all RNA samples were determined to be between 1.9 and 2.0, suggesting that all RNA samples were highly pure. RNA integrity was verified by the Agilent RNA 6000 Nano assay kit (Agilent Bioanalyser 2100, Agilent Technologies, Cheadle, UK) as described by the manufacturer. Single-strand cDNA used for RT-qPCR analyses was synthesized from purified total RNA using SuperScript <sup>R</sup> III Reverse Transcriptase (Life Technologies, Paisley, UK).

encoded by shape. Fragments were added in pools of 6. Pools in which there was a hit relative to the same control are marked in red. The line is a line of best fit.

RT-qPCR were performed using 384-well plates, with a final volume of 10 µL in each well, consisting of 4 µL of cDNA, 5 µL of 2x SYBR Green LightCycler 480\_TM PCR master mix (Roche Life Sciences), 0.8 µL of sterile distilled water, 0.1 µL each of 20 µM reverse and forward primers. Samples were performed in triplicates. In the no template controls (negative controls) 4µL of H2O were added, instead of the cDNA samples. RT-qPCR reactions were carried-out using the Roche LightCycler LC\_480 qPCR platform, where fluorescence signals were measured in

real-time. The protocol, set-up with thermal cycling conditions, consisted of one cycle at 95◦C for 10 min, followed by 45 cycles of amplification at 95◦C for 10 s, and 60◦C for 30 s. Roche LightCycler Data Analysis Software was used to determine the melt curve data as well as the quantification cycle values (Cq values). The changes in expression levels were normalized against two reference gene as determined via GeNorm (REF), and the relative mRNA levels of genes following treatment

were calculated using "The Comparative C<sup>T</sup> Method" (11C<sup>T</sup> Method).

## Design of Primers for RT-qPCR

The National Centre for Biotechnology Information (NCBI) website (http://www.ncbi.nlm.nih.gov/) was used to identify and obtain mRNA sequences. Exon boundaries were determined from the "European Molecular Biology Laboratories" website (http://www.ensembl.org). This procedure was performed until sets of primers were selected for each target gene. The final step involved checking the primers for similarity using NCBI BLAST (Basic Local Alignment Search Tool) (http://www.ncbi. nlm.nih.gov/BLAST), reducing the chance of primers binding non-specifically.

### Identification of Reference Genes for RT-qPCR Analysis

Samples were analyzed for the expression of each of eight candidate reference genes, namely: ACTB (Beta-Actin), B2M (Beta-2-microglobulin), GAPDH (glyceraldehyde-3-phosphate dehydrogenase), HMBS (hydroxymethyl-bilane synthase), HPRT1 (hypoxanthine phosphoribosyl transferase 1), RPL13A (ribosomal protein L13a), RPL32 (ribosomal protein L32), SDHA (succinate dehydrogenase complex, subunit A) as recommended by Vandesompele et al. (2002). RT-qPCR was performed as described previously, using the primers specific for each candidate reference gene. The GeNorm algorithm software package was used to determine the two most stable reference genes from the set of tested candidate genes by calculating a gene normalization factor, eliminating the least stable genes until a stability value (M) of 0.4 or less was reached (Vandesompele et al., 2002).

## RESULTS

## Effects of Gemcitabine and Drug Fragments on the Viability of Panc-1 Cells

A standard strategy is to choose a series of molecules that cover chemical space effectively, and for this we chose initially the main Maybridge drug fragment library. It consists of 500 rule-of-three-compliant (Congreve et al., 2003) polar molecules that cover chemical space widely, and where the molecular properties include molecular weight <300, number of hydrogen bond donors ≤3, number of hydrogen bond acceptors ≤3, ClogP ≤3, and in addition, the number of rotatable bonds ≤3 and the polar surface area ≤60Å<sup>2</sup> . While the use of fragments is commonplace in target-based assays, especially where structures are known (e.g., Erlanson and Hansen, 2004; Erlanson et al., 2004; Rees et al., 2004; Carr et al., 2005; Alex and Flocco, 2007; Ciulli and Abell, 2007; Jhoti, 2007; Jhoti et al., 2007; Hubbard, 2008; Fischer and Hubbard, 2009; Schulz and Hubbard, 2009; Whittaker et al., 2010; Leach and Hann, 2011; Erlanson, 2012; Caliandro et al., 2013), we here prefer the use of the rather more successful phenotypic screens (Swinney and Anthony, 2011; Swinney, 2013). Although it is hard to find published examples of phenotypic screens that used fragment-based libraries, we merely point out that 25% of successful (marketed) drugs are no larger than fragments (i.e., <300 Da) (O'Hagan and Kell, 2015c). The fragment-based approach also has the advantage of avoiding the

#### TABLE 2 | Six hits in the "binary weapon" assay given in three formats, plus indole-3-carbinol.


increasing "molecular obesity" (Hann, 2011; Meanwell, 2011) that is seen in some cases as inimical to the finding of successful drugs (Leeson and Springthorpe, 2007; Leeson and Empfield, 2010).

Panc1 cells are a pancreatic cancer cell line (e.g., Gradiz et al., 2016). **Figure 1** shows four separate experiments in which the effect of the pools of the Maybridge fragments (6 at a time) on the viability of cells was assessed in the presence and absence of 20 nM gemcitabine, pointing up three pools containing "hits" (which occurred in at least 3 experiments; there are a total of 336 experiments here). **Figure 2** shows the % viability of one set of Panc1 cells as a function of the gemcitabine concentration, as a result of which we later chose 100 nM gemcitabine to assess the efficacy of the individual fragments in increasing its toxicity. **Figure 3** shows a titration curve for three repeats with one of the "hits," the plot also serving to illustrate the variability of the toxicity of gemcitabine alone on different days. **Figures 4**, **5** show the distribution in chemical space of all 500 fragments in the first Maybridge library and three "hits" at 10 µM that lowered the viability of cells by at least 10% in the presence, but not the absence, of 100 nM gemcitabine. These were retested singly, then together pairwise, resulting in three hits, viz B1, D1, and B12. B12 seemed to interfere with the other two fragments by binding to them directly (UV evidence) and was not used further. Note that a significant issue is that although for a given batch of Panc1 cells the titration curves were reasonably reproducible, they were considerably less so between batches (for reasons that will become apparent below). This meant that each culture had to be used as its own control, as we did e.g., in **Figure 2**. Another interesting feature was that quite a significant fraction of the fragments (as in **Figure 1**, and see below) were even somewhat stimulatory to cell growth in the absence of gemcitabine.

There are four other Maybridge fragment libraries of 500 molecules each, covering broadly the same chemical space but in more detail (O'Hagan and Kell, 2015c), and we performed a cheminformatics analysis (MACCS encoding, Tanimoto similarity) to establish which other molecules might be similar, exactly as per the analyses in (O'Hagan et al., 2015). Some 20

clusters of molecules that could be observed.

molecules had a Tanimoto similarity within 0.7 of one of the three remaining hits and were tested. In this case, the starting % viability was much higher than those in **Figure 2**. All 20 of these fragments are in fact active, which shows that these molecules (**Figure 6**) exhibit a very considerable enrichment over the whole library, and illustrates the utility of the principle of molecular similarity (Gasteiger, 2003; Bender and Glen, 2004; Stumpfe and Bajorath, 2011; Maggiora et al., 2014). The figure also illustrates which of the original three hits the new hits are closest to, and encodes their S log P-values as the size of the marker. This enormous cheminformatics-based enrichment also gives considerable confidence in our strategy, despite the variability in sensitivity of the Panc1 cells to gemcitabine alone, since such a huge enrichment could not conceivable occur for molecules that were not active. Although none was quite as active as the original hits, all exhibited some kind of activity (**Figure 2**) (the starting viabilities for two different experiments in the presence of gemcitabine only were 78 and 84%). Of all of these, the seven most potent molecules exhibited activity at 3 µM. One was rather expensive and was again excluded. Thus, we had a total of 6 hits to consider [two from library 1 (B1 and D1), and a total of four from the other four libraries, referred to as fragments 10, 11, 12, and 20]. **Table 2** gives their names, SMILES encodings and 2D structures, along with that of indole-3-carboxylic acid (see later).

**Figure 7** shows a (symmetrical) heatmap (MACCS encoding) of the Tanimoto similarities of the 22 most potent molecules, where it can again be seen that the hits are in three clusters. These are B1, 10, and 11 (all are amines), D1 and 12 (carboxylic acids), and 20 (an aniline derivative—possibly to be avoided Benigni and Passerini, 2002; Benigni et al., 2009; Franke et al., 2010). One implication is that they each have different targets (probably plural) but attempts even to show additivity, let alone synergy, met with failure, possibly because the molecules were indeed rather similar to each other in terms of the larger chemical space. **Figure 8**—equivalent to **Figure 3**—shows data for two experiments with fragment 10, again illustrating the stimulation of growth by the fragment alone, and its inhibition in the presence of a relatively weakly inhibiting concentration of gemcitabine. Finally, **Figure 9** shows the Tanimoto similarities (TS, based on the MACCS encoding) between the six hits plus Indole-3-carbinol (I3C, see below). Fragments within a group showed a Tanimoto similarity of 0.75 or greater, while those between groups were less than 0.5. I3C was not really similar to any of the hits; its highest TS to any of the hits was 0.36. It is

especially gratifying to note that MBF10 and MBF11 were both selected and had a TS to each other of 1, as they are in fact structural isomers. Along with the other clusterings, this adds considerable weight to the validity of our assays.

### Effect of Indole-3-Carbinol on Gemcitabine Toxicity

Cruciferous vegetables such as Brassica spp. are considered to have certain anticancer properties (Higdon et al., 2007; Juge et al., 2007; Fujioka et al., 2016b), and small molecules derived from the hydrolysis of glucosinolates, such as sulforaphane and indole-3-carbinol (I3C), have been implicated in a variety of anticarcinogenic mechanisms (e.g., Chen et al., 2014; Fujioka et al., 2016a). I3C is a small molecule (MW 147.17, well within the range of "fragments"), and Lyn-Cook and colleagues (Lyn-Cook et al., 2010; Wang et al., 2011; Paik et al., 2013) have published that I3C can enhance the sensitivity of pancreatic cancer cells to gemcitabine, possibly via upregulation of ENT1 expression (Wang et al., 2011). It was thus of interest to compare I3C with the hits that we found. In our hands, however, I3C had no measurable effect on either the cell viability in the presence or absence of gemcitabine (nor on the expression profiles discussed below). This is entirely consistent with its low structural similarity to the other hits as indicated above.

## Effect of Fragments on the Growth of Panc1 Cells

Although this was not the main focus of the present paper, we did note (as mentioned above) that the fragments themselves could stimulate the growth of Panc1 cells relative to that of controls (as measured by OD). This is illustrated in **Figure 10** for 28 of the fragments on which we focussed. Also encoded with the structures are the number of H-bond donors and acceptors, the total polar surface area of the fragments, and (on the abscissa) the S log P-values. It is clear (i) that virtually every fragment could stimulate the growth of the cells, and (ii) that there was no particularly obvious relationship of the extent of such stimulation with any of the descriptors stated.

### Effect of Gemcitabine and Fragments on the Expression of Selected Transcripts in Panc 1 Cells

Given that there was evidence that the fragments did not affect gemcitabine uptake directly, we assumed that they must be working by influencing the activity or expression of appropriate targets (and certainly small molecules can affect transporter expression, (e.g., Mrozikiewicz et al., 2014). To this end, we

square 1, circle 2, diamond 3, triangle 4), H-bond donors (by color, blue 0, green 1, red 2, yellow 3), total polar surface area (by size of symbol, up to 63 Å<sup>2</sup> ) and S log P (on the abscissa).



*Only those transcripts detectable within 35 PCR cycles are shown. Data are given as mean* ± *standard deviation. A 2-sided T-test was performed to assess statistical significance against GEM alone, P-values being encoded as* \* < *0.05,* \*\* < *0.01,* \*\*\* < *0.001.*

TABLE 4 | Changes in the transcript level of ABCC2 and RRM1 when treated with gemcitabine and/or the other fragment hits.


*A 2-sided T-test was performed to assess statistical significance vs. GEM alone, P-values being encoded as* \* < *0.05,* \*\* < *0.01,* \*\*\* < *0.001.*

designed primers to enable PCR of transcripts relevant to gemcitabine transport and metabolism. **Table 3** shows each of those that were detectable within 35 PCR cycles when treated (i) with gemcitabine alone, (ii) with Maybridge fragment D1 alone, and (iii) with both gemcitabine and D1. Strikingly, gemcitabine increases the expression of the ABCC2 efflux transporter (MRP2) more than 12-fold, and that of RRM1 more than fourfold, while the addition of D1 largely reverses both of these effects. It would seem that these are by far the largest contributors to the efficacy of fragment D1 in enhancing the cytotoxicity of gemcitabine, and the same is true for each of the other fragments (**Table 4** and **Figure 11**). However, the ABCC2 inhibitor MK-571 (e.g., Weiss et al., 2007; Noma et al., 2008) at 20 µM had no effect on the viability of Panc1 cells treated with Gemcitabine alone (data not shown), possibly implying that RRM1 was the more significant contributor to the phenotypic changes in resistance.

### Selectivity of Fragments for Increasing Transporter Expression

Having seen that various of the fragments could increase the toxicity of gemcitabine to Panc1 cells, it was of interest to see whether this was a cell-selective phenomenon. Although time did not permit an exhaustive study, we noted that fragments 10 and 20 also had these toxicity-enhancing effect for the neuroblastoma SH-SY5Y cell line while B1, D1, 11, and 12 did not (**Figure 12**). No fragments seemed to have any such effects on the noncancerous pancreatic cell line hPDE (**Figure 13**) and HEK293 cells (**Figure 14**), implying that there is or can be at least some degree of specificity in our "binary weapon" approach. Clearly a larger-scale study (including both larger libraries and more cell lines) would be able to discover molecules with both potency and selectivity.

### DISCUSSION

In the present work, we sought to develop the idea that we might affect the transporter-mediated disposition of small-molecule drugs via the addition of a second small molecule that of itself had no inhibitory pharmacological effect but that influenced the expression of transporters for the primary drug (**Figure 15**). We refer to this as a "binary weapon" strategy. The specific phenotypic effect we sought was for a molecule that on its own had no such effect to increase the toxicity of the nucleoside analog gemcitabine to Panc1 pancreatic cancer cells (**Figures 1**–**3**).

Given the recognition (O'Hagan and Kell, 2015c) that more some 25% of marketed drugs are in fact no larger than the polar "rule-of-three"-compliant (Congreve et al., 2003) molecules used in fragment-based drug discovery, we used an initial screen of a 500-member polar drug fragment library. This yielded three "hits" (**Figures 4**, **5**). The structures of 20 of the other 2000 members of this library had a Tanimoto similarity greater than 0.7 to those of the initial hits, and each was itself a hit (**Figure 6**) (with the cheminformatics thus providing for a massive enrichment in the fraction of successful experiments). We chose the top six representatives for further study. They each bore reasonable structural similarities to each other (two were in fact isomers), lending strength to the self-consistency of both our conceptual and experimental strategies (**Figures 7**, **8**).

Existing literature had suggested that indole-3-carbinol might play a similar role to that of our fragments, but in our hands it was without effect, and nor was it structurally similar to any of our hits (**Figure 9**). We therefore discounted it.

There is an interesting issue when the phenotypic activity being measured is in fact cell death, as it is then impossible legitimately to compare bulk measurements of biochemical changes with individual-cell viabilities. This is because with bulk or ensemble measurements one does not know if say a lowering of a biochemical parameter by 50% means that all of the cells have lost half of it or half of the cells have lost all

of it (or anything in between) (Kell et al., 1991, 1998; Davey and Kell, 1996). In the event, the mechanism was very clear, however.

Because the fragments were themselves without negative effects on the cells in the absence of gemcitabine (interestingly, many of them actually stimulated cell growth, **Figure 1**, so each had to be compared to the appropriate control!), we next designed suitable primers to assess the expression levels of all the candidate transporters plus ribonucleotide reductase. In our hands, only the ENT1-3 uptake and ABCC2,3,4,5, and 10 efflux transporters displayed measurable transcripts, along with RRM1.

### REFERENCES


Very strikingly, the addition of gemcitabine alone increased the expression of the transcript for ABCC2 (MRP2) by more than 12-fold, and that of RRM1 by more than fourfold, and each of the fragment "hits" served to reverse this, at least in part (**Figure 11**). The effects on ABCC2 are thus consistent with the finding (Horiguchi et al., 2013) that it may be a major efflux pump for gemcitabine.

It seems, therefore, that while the effect was here mediated more by efflux than influx transporters, the binary weapon idea is hereby fully confirmed: our results show that it is possible to find molecules that manipulate the expression of transporters that are involved in the bioactivity of a pharmaceutical drug, and that there is a certain degree of specificity in this for pancreatic cancer cells (**Figures 12**–**14**). This could explain, at least in part, the basis for the selective toxicity of a drug that is otherwise cytotoxic generally (**Figure 15**). The next steps will involve determining much more extensively how much any such activity differs, or can be made to differ (as do most transcript levels), between different cells.

### AUTHOR CONTRIBUTIONS

DK and PD designed the study. All the experimental work was performed by JG, who was supervised by PD and DK. Some of the cheminformatics analyses were performed by SO. All authors contributed to and approved the writing of the manuscript.

### ACKNOWLEDGMENTS

We thank the Biotechnology and Biological Sciences Research Council for supporting this work (BBSRC grants BB/K019783/1 and BB/M017702/1). This is a contribution from the Manchester Centre for Synthetic Biology of Fine and Speciality Chemicals (SYNBIOCHEM).


cells, and enhances adenovirus-mediated oncolysis. Cancer Biol. Ther. 15, 1256–1267. doi: 10.4161/cbt.29690


gemcitabine in patients with pancreatic cancer. Gastroenterology 136, 187–195. doi: 10.1053/j.gastro.2008.09.067


dFdCTP metabolite disposition in patients with solid tumours. Br. J. Cancer 110, 304–312. doi: 10.1038/bjc.2013.738


with gemcitabine-cisplatin-based combination chemotherapy. BJU Int. 108, E110–E116. doi: 10.1111/j.1464-410X.2010.09932.x


regulatory subunit M1 (RRM1) expression is a powerful predictor of survival in patients with pancreatic carcinoma treated with adjuvant gemcitabine-based chemotherapy after operative resection. Surgery 153, 565–575. doi: 10.1016/j.surg.2012.10.010


barrier permeability in the rat embryo. Methods Mol. Biol. 686, 247–265. doi: 10.1007/978-1-60761-938-3\_11


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2017 Grixti, O'Hagan, Day and Kell. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

,

# Protocatechuic Aldehyde Attenuates Cisplatin-Induced Acute Kidney Injury by Suppressing Nox-Mediated Oxidative Stress and Renal Inflammation

, Gui-Ling Ren1,2,3, Hai-Di Li1,2,3, Qin Yang1,2,3

,

#### Edited by:

Li Gao1,2,3†

Xiao-Ming Meng1,2,3

, Wei-Feng Wu1,2,3†

\*

Thomas J. Anastasio, University of Illinois at Urbana–Champaign, USA

#### Reviewed by:

Mariarosaria Santillo, University of Naples Federico II, Italy Md Abdul Hye Khan, Medical College of Wisconsin, USA

#### \*Correspondence:

Xiao-Ming Meng mengxiaoming@ahmu.edu.cn †These authors have contributed equally to this work.

#### Specialty section:

This article was submitted to Experimental Pharmacology and Drug Discovery, a section of the journal Frontiers in Pharmacology

> Received: 10 August 2016 Accepted: 23 November 2016 Published: 06 December 2016

#### Citation:

Gao L, Wu W-F, Dong L, Ren G-L, Li H-D, Yang Q, Li X-F, Xu T, Li Z, Wu B-M, Ma T-T, Huang C, Huang Y, Zhang L, Lv X, Li J and Meng X-M (2016) Protocatechuic Aldehyde Attenuates Cisplatin-Induced Acute Kidney Injury by Suppressing Nox-Mediated Oxidative Stress and Renal Inflammation. Front. Pharmacol. 7:479. doi: 10.3389/fphar.2016.00479 <sup>1</sup> School of Pharmacy, Anhui Medical University, Hefei, China, <sup>2</sup> Anhui Institute of Innovative Drugs, Hefei, China, <sup>3</sup> Key Laboratory of Anti-inflammatory and Immune Medicine, Ministry of Education, Hefei, China, <sup>4</sup> Department of Pediatrics, Division of Hematology/Oncology, Aflac Cancer and Blood Disorders Center, Children's Healthcare of Atlanta, Emory University School of Medicine, Atlanta, GA, USA

, Lei Dong<sup>4</sup>

Xiao-Feng Li1,2,3, Tao Xu1,2,3, Zeng Li1,2, Bao-Ming Wu1,2,3, Tao-Tao Ma1,2,3

Cheng Huang1,2,3, Yan Huang1,2,3, Lei Zhang1,2,3, Xiongwen Lv1,2,3, Jun Li1,2,3 and

Cisplatin is a classic chemotherapeutic agent widely used to treat different types of cancers including ovarian, head and neck, testicular and uterine cervical carcinomas. However, cisplatin induces acute kidney injury by directly triggering an excessive inflammatory response, oxidative stress, and programmed cell death of renal tubular epithelial cells, all of which lead to high mortality rates in patients. In this study, we examined the protective effect of protocatechuic aldehyde (PA) in vitro in cisplatintreated tubular epithelial cells and in vivo in cisplatin nephropathy. PA is a monomer of Traditional Chinese Medicine isolated from the root of S. miltiorrhiza (Lamiaceae). Results show that PA prevented cisplatin-induced decline of renal function and histological damage, which was confirmed by attenuation of KIM1 in both mRNA and protein levels. Moreover, PA reduced renal inflammation by suppressing oxidative stress and programmed cell death in response to cisplatin, which was further evidenced by in vitro data. Of note, PA suppressed NAPDH oxidases, including Nox2 and Nox4, in a dosage-dependent manner. Moreover, silencing Nox4, but not Nox2, removed the inhibitory effect of PA on cisplatin-induced renal injury, indicating that Nox4 may play a pivotal role in mediating the protective effect of PA in cisplatin-induced acute kidney injury. Collectively, our data indicate that PA blocks cisplatin-induced acute kidney injury by suppressing Nox-mediated oxidative stress and renal inflammation without compromising anti-tumor activity of cisplatin. These findings suggest that PA and its derivatives may serve as potential protective agents for cancer patients receiving cisplatin treatment.

Keywords: protocatechuic aldehyde, acute kidney injury, Nox, oxidative stress, inflammation, necroptosis

### INTRODUCTION

fphar-07-00479 December 3, 2016 Time: 13:59 # 2

Cisplatin is widely used in the treatment of various cancers including ovarian, head and neck, testicular and uterine cervical carcinomas (Pabla and Dong, 2008; Sung et al., 2008). Although regarded as one of most effective chemotherapeutic agents by directly interfering with DNA synthesis and inducing apoptosis, adverse effects such as nephrotoxicity put this promising anticancer agent in a precarious position. In fact, approximately 30% of patients experience a marked decline in renal function after a single dose injection of cisplatin (Sung et al., 2008). To this point it is important to prevent cisplatin-induced acute kidney injury, which is growing in clinical significance. Although emerging evidence indicates several pathological mechanisms, including excessive inflammatory response, oxidative stress, apoptosis and death of renal tubular epithelial cells, clear therapeutic targets and effective therapies are still lacking (Sancho-Martinez et al., 2015; Yang et al., 2016; Zuk and Bonventre, 2016).

Traditional Chinese Medicine (TCM) may be a promising and novel avenue to effectively treat human diseases with lower toxicity. Indeed, artemisinin is an impressive example of success in the war against malaria. Previous studies indicated that several TCM monomers, such as resveratrol (Kim et al., 2011), Luteolin (Domitrovic et al., 2013), and Emodin (Liu et al., 2016) relieved cisplatin-induced acute kidney injury in animal models. Our recent findings also revealed that 18β-glycyrrhetinic acid alleviated cisplatin-induced apoptosis of renal tubular epithelial cells by targeting HDAC2/BMP-7 axis (Ma et al., 2016). However, the identification of more efficient and low toxic therapeutic agents needs more research. In our pilot study, we tested 10 potential anti-inflammatory TCM monomers including aloin, barbaloin, icariin, protocatechuic acid, protocatechuic aldehyde (PA), puerarin, sodium houttuyfonate, sophoridine, wogonin, and wogonoside, which have not been investigated in the kidney field or in cisplatin-treated tubular epithelial cells (data not shown). Results presented here show that PA, isolated from the root of S. miltiorrhiza, is one of the most powerful protective TCM monomers. It significantly suppressed cisplatin-induced injury of tubular epithelial cells and the inflammatory response. This was possibly mediated by inhibition of oxidative stress and programmed cell death. Moreover, the protective role of PA was further evidenced in vivo in cisplatin nephropathy, where it prevented decline of renal function and attenuated renal injury. More importantly, results of MMT assay in three tumor cell lines demonstrated that treatment of PA didn't alter the antitumor property of cisplatin. These findings indicate that PA may be a potential therapeutic agent for preventing cisplatin-induced acute kidney injury.

### RESULTS

### PA Ameliorated Cisplatin-Induced Death in HK2 Cells

We used an MTT assay to analyze the impact of PA on cell viability in the human tubular epithelial cell line (HK2). Results show that PA treatment began to reduce cell viability at concentrations greater than 1 µM (**Figure 1A**). Moreover, PA in concentration of 0.25, 0.5, and 1 µM significantly restored cell viability after cisplatin treatment (20 µM; **Figure 1B**). Given these results, we chose 0.25, 0.5, and 1 µM PA for subsequent experiments. We also determined whether PA limited the antitumor activity of cisplatin in three solid tumor cells lines, SMCC-7721, BEL-7402, and U87. MTT assay data show that cisplatin reduced cell viability of hepatic cancer SMCC-7721 cells, particularly at 48 h, and administration of PA didn't reduce the tumor-suppressive effect of cisplatin (**Figure 1C**). This was further supported by the findings that PA didn't protect against cisplatin-induced tumor cell death in human hepatic cancer line BEL-7402 and malignant gliomaU87 cell line.

### PA Protected against Cisplatin-Induced Cell Damage and Inflammatory Response

To assess whether PA reduces kidney damage, we examined mRNA and protein expression of kidney injury molecule-1(KIM1). Western blot and real-time PCR results show that cisplatin upregulated KIM1. This was decreased by PA treatment in a time-dependent manner in HK-2 cells (**Figures 2A,B**). Additionally, real-time PCR and ELISA analysis show that PA protected against inflammatory response as evidenced by decreased chemokine monocyte chemotactic protein (MCP-1), inflammatory cytokine (IL-8), and proinflammatory cytokine TNF-α expression levels (**Figures 2B,C**).

### PA Inhibited Cisplatin-Induced Cell Necroptosis and Apoptosis

We tested the protective effects of PA on cell death of HK2 by flow cytometric analysis of PI/AnnexinV staining. Results show that PA alleviated cisplatin-induced necroptosis and apoptosis (**Figure 3A**). Mechanistically, PA significantly reduced the key signaling molecules mediating necroptosis, including RIP1, RIP3, and phosphorylation of downstream MLML in HK2 cells (**Figure 3B**). Furthermore, cleaved-caspase-8, cleaved-caspase-3, cleaved-caspase-12, and phosphorylation of p53 were also markedly decreased in response to PA treatment (**Figure 3B**).

### PA Suppressed Cisplatin-Induced Injury In vitro via Blocking Nox-Mediated Oxidative Stress

A reactive oxygen species (ROS) assay was performed using DCF fluorescence in HK2 cells to measure the effect of PA on oxidative stress. Results show that PA largely decreased ROS in cisplatin-stimulated HK2 cells (**Figure 4A**). Moreover, results of DHE staining show that superoxide levels, the reaction product of Nox enzymes, were suppressed by PA treatment (**Figure 4B**). We know that apocynin or PA treatment decreases cisplatininduced damage of tubular epithelial cells. Interestingly, we found that after blocking Nox enzymes with apocynin, PA failed to further reduce the renal damage (**Figure 4C**). This indicates a role for the Nox family in mediating the protective effect of PA. Therefore, we analyzed the main members of the Nox family in kidney, including Nox2 and Nox4. Western blot results show that

PA treatment reduced cisplatin-induced Nox2 and Nox4 protein expression in a dose-dependent manner (**Figures 4D,F**). This was further supported by real-time PCR data showing reduced Nox2 mRNA and Nox4 mRNA (**Figures 4E,G**).

### PA Attenuated Cisplatin-Induced Injury through Nox4-Dependent Mechanisms

We determined whether Nox2 and Nox4 are important targets for PA in mediating its protective effects. We found that Nox2 and Nox4 were largely suppressed when shRNA plasmids were transfected into cells (**Figures 5A,B**). Western blot results show that PA further suppressed KIM1 levels in absence of Nox2. But, PA failed to further attenuate renal tubular cell damage when Nox4 was disrupted (**Figure 5C**). This indicates an essential role for Nox4 in mediating protective effects of PA.

### PA Attenuated Cisplatin-Induced Necrosis, Inflammation, and ROS through Nox4-Dependent Mechanisms

We analyzed the role of PA on the Nox4 pathway in more detail in vitro. We found that PA failed to reduce cleavedcaspase-3 levels in Nox4 knockdown HK2 cells (**Figure 6A**). Consistently, when Nox4 was silenced, PA didn't further suppress the production of inflammatory cytokines (**Figure 6B**). This indicates that PA protects against cisplatin-treated tubular

epithelial cells in a Nox-dependent manner. These results were further confirmed by detecting MDA level, which showed that ROS levels were consistent (**Figure 6C**).

### PA Inhibited Cisplatin-Induced Acute Kidney Injury in Mice

The protective effect of PA was also evaluated in vivo in cisplatin nephropathy. Results of periodic acid-Schiff (PAS) staining revealed that administration of PA in concentrations of 0.45, 0.9, and 1.8 mg/kg reduced tubular necrosis, dilation, and cast formation compared with model groups (**Figure 7A**). The therapeutic effects of PA were further confirmed by detection of renal function, including serum creatinine and blood urea nitrogen. Results show that PA largely prevented the decline of renal function in a dose-dependent manner (**Figures 7B,C**).

### PA Significantly Reduced Cisplatin-Induced Tubular Injury and Inflammation Response in Mice

Western blot and quantitative data show that KIM1, a key marker for tubular injury, was significantly elevated by cisplatin but reduced by PA treatment in a dose-dependent

phosphorylation of p53 in HK2 cells; Data represent the mean ± SEM for 3–4 independent experiments. <sup>∗</sup>p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001 compared to the control. #p < 0.05, ##p < 0.01, ###p < 0.001 compared to cisplatin-treated group. Cis, cisplatin; PA, protocatechuic aldehyde.

manner (**Figure 8A**). This was further evidenced by results from IHC and real-time PCR (**Figures 8C,E**). PA treatment also suppressed the increased level of Neutrophil gelatinaseassociated lipocalin (NGAL), another kidney injury marker, in urine in cisplatin nephropathy (**Figure 8B**). Additionally,

results from IHC show that PA reduced TNF-α positive signals in injured kidney. This was consistent with real-time PCR results showing reduced mRNA levels of proinflammatory cytokines and chemokines including TNF-α, IL-6, and MCP-1 (**Figures 8D,E**).

FIGURE 4 | Protocatechuic aldehyde inhibited cisplatin-induced cell oxidative stress in HK2 cells by suppressing Nox signaling. (A) DCF Assay of Reactive Oxygen Species. Results of ROS assay demonstrate that PA inhibited cisplatin-induced cell oxidative stress in HK2 cells; (B) DHE staining of intracellular ROS levels. PA attenuates cisplatin-induced ROS generation. (C) Role of Nox enzymes in the effect of PA on cisplatin-treated HK2 cells. (D,F) Western blot analysis and quantitative data of Nox2 and Nox4 in cisplatin-treated HK2 cells. Results show that administration of PA substantially suppressed protein levels of both Nox2 and Nox4 in cisplatin-stimulated HK2 cells. (E,G) Real-time PCR in HK2 cells. Results demonstrate that treatment of PA largely reduced cisplatin-induced mRNA levels of Nox2 and Nox4; Data represent the mean ± SEM for 3–4 independent experiments. <sup>∗</sup>p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001 compared to the control. #p < 0.05, ##p < 0.01, ###p < 0.001 compared to cisplatin-treated group. Cis, cisplatin; PA, protocatechuic aldehyde; APO, Apocynin.

### PA Protected against Cisplatin Nephropathy by Attenuating Oxidative Stress

We then investigated the underlying mechanisms by which PA improves renal function and limits kidney injury. We found that MDA levels were significantly increased in kidneys of cisplatin-injected mice. But, PA reduced the upregulation of MDA levels in a dose-dependent manner (**Figures 9A,B**). This is consistent with the finding that PA restored the cisplatinsuppressed level of GSH, an antioxidant index. This indicates that PA substantially decreased cisplatin-induced production of ROS. Moreover, Western blot and quantitative data show that PA reduced protein levels of Nox2 and Nox4 in a dosagedependent manner in cisplatin nephropathy. This was further confirmed by immunohistochemistry that revealed that Nox4 was downregulated in cisplatin-injured kidney treated with PA (**Figures 9C,D**). Additionally, the function of Nox4 was detected in vivo. Results of PAS staining indicate that lentivirus-mediated knockdown of Nox4 significantly attenuated cisplatin-induced kidney damage (**Figure 9E**). This demonstrates it may be a critical target to mediate protective effects of PA.

### PA Protected against Cisplatin Nephropathy by Attenuating Cell Death

Our results show that PA suppressed activation of the RIP1/RIP3/MLKL axis, which is regarded as the key pathway mediating necroptosis. We found PA also decreased levels of cleaved-caspase-3 in cisplatin nephropathy (**Figure 10**).

disrupted HK2 cells. (A) Western blot analysis and quantitative data of cleaved-caspase-3. (B) Real-time PCR analysis of inflammatory response in Nox4 disrupted HK2 cells. Results show that when Nox4 was disrupted, PA failed to further decrease mRNA levels of KIM1, MCP-1, TNF-α, and IL-8. (C) Malondialdehyde (MDA) levels in HK2 cells. Results indicate that when Nox4 was knocked down, PA failed to further decrease the MDA levels in cisplatin-stimulated HK2 cells. Data represent the mean ± SEM for 3–4 independent experiments. <sup>∗</sup>p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001 compared to the control. #p < 0.05, ##p < 0.01, ###p < 0.001 compared to Nox4 EV group. \$p < 0.05, \$\$p < 0.01, \$\$\$p < 0.001 compared to cisplatin-treated group. Cis, cisplatin; PA, protocatechuic aldehyde; EV, empty vector; KD, knockdown.

### MATERIALS AND METHODS

### Murine Model of Cisplatin-Induced AKI

Mice were obtained from Laboratory Animal Center of Anhui province. All animal procedures were approved by the Institutional Animal Experimentation Ethics Committee of Anhui Medical University. Cisplatin with a single dose at 20 mg/kg was injected intraperitoneally in 8-week-old male mice and their littermates injected with saline were set as normal control. Protocatechuic aldehyde (Dalian Meilun Biotech Co. Ltd.), concentrations of 0.45, 0.9, and 1.8 mg/kg, were given via intraperitoneal injection 6 h before cisplatin treatment and injected daily. Mice were sacrificed under anesthesia (10% chloralic hydras by intraperitoneal injection) 3 days after cisplatin injection. Samples of kidney tissues and blood were harvested for further analysis, including BUN (Nanjing Jiancheng Bioengineering Institute) and creatinine, paraffin embedding (Nanjing Jiancheng Bioengineering Institute) and molecular analysis. Paraffin sections (3–5 mm) were stained with Periodic acid Schiff (PAS) Staining kit (Fuzhou Maixin Biotech. Co., Ltd.) and analyzed via immunohistochemistry.

### Reagents and Materials

Antibodies RIP1, RIP3, KIM-1, Nox2, Nox4, TNF-α, and β-actin were purchased from Santa Cruz Biotechnology (Santa Cruz, CA, USA); Rabbit anti-P-MLKL and Anti-cleavedcaspase-3, cleaved-caspase-12 were obtained from Cell Signaling

Technology (CST, Danvers, MA, USA). Lipofectamine 3000 was purchased from Science Biotechnology (Invitrogen, Beijing, China). Protein Assay Kit was purchased from Beyotime Institute of Biotechnology (Jiangsu, China). Cell Malondialdehyde (MDA) assay kit (Colorimetric method), reduced glutathione (GSH) assay kit, creatinine (Cr) Assay kit (sarcosine oxidase) and urea (BUN) assay kit were obtained from Nanjing Jiancheng Bioengineering Institute(Nanjing, China). Reactive Oxygen Species Assay (DCF Assay) Kit and Dihydroethidium (DHE) were purchased from Beyotime Institute of Biotechnology (Jiangsu, China). ITC AnnexinV/propidiuiodide was purchased from Bestbio (Shanghai, China). Mouse NGAL ELISA kit was purchased from CUSABIO (Wuhan, China).

### Cell Culture

Solid tumor cell lines (SMCC-7721, BEL-7402, and U87), kidney tubular epithelial cells of human (HK2), and Nox2, Nox4 knockdown HK2 cells were cultured in 5% FBS-containing HyCloneTM DMEM/F12 medium at 37◦C in humidified 5% CO2. After overnight starving in DMEM/F12 medium containing 0.5% FBS, HK2 cells were pretreated with PA (Meilun Biology Technology, Dalian, China) for 6 h before being exposed to cisplatin (20 µM). Twenty-four hour later, the cells were harvested for cell viability, the indexes of kidney injury, oxidative stress, programmed cell death, and inflammatory response such as KIM1, Nox2, Nox4, cleaved-caspase-3, 8, and 12, p-p53, p53, RIP1, RIP3, phospho-MLKL, TNF-α, and MCP-1 using Western blot analysis, real-time PCR or other methods. Three to four in vitro experiments were performed independently.

### Cell Viability Assay

Cell viability is determined by MTT assay, according to a purple formazan product produced by mitochondrial dehydrogenase of viable cells. Human HK2 cells were grown in 96-well plates with treatments by a set of concentrations of PA (arranged from 0.25 to 8 µM). After 12 h, the cells were exposed to cisplatin (20 µM) for 24 h in incubator. They were harvested after addition of 5 mg/ml MTT solution to each well for 4 h. Optical density (OD) was determined in microplate reader (Multiskan MK3, Thermo, USA) at 492 nm wavelength.

### RNA Extraction and Real-Time PCR

Total RNA was isolated using the RNeasy Isolation Kit (Qiagen, Valencia, CA, USA) according to the manufacturer's

instructions (Li et al., 2012). The RNA concentration was detected by a NanoDrop 2000 Spectrophotometer (Thermo Scientific, USA). RNA, nuclease-free water and RealMasterMix (Bio-Rad, Hercules, CA, USA) were used for cDNA synthesis. Real-time PCR was performed in a total volume of 9 µl, including 2 µl cDNA solution, 4 µlBio-Rad iQ SYBR Green supermix with Opticon 2 (Bio-Rad, Hercules, CA, USA), 2.4 µl nuclease-free water,

and 0.6 µl each primer. The sequences of primers are as follows:

Human IL-8, forward 5<sup>0</sup> -AGGACAAGAGCCAGGAAGAA-3<sup>0</sup> , reverse 5<sup>0</sup> - ACTGCACCTTCACACAGAGC-3<sup>0</sup> ; Human TNF-α, forward 5<sup>0</sup> -CCCAGGGACCTCTCTCTAATCA-3<sup>0</sup> , reverse 5<sup>0</sup> - GCTACAGGCTTGTCACTCGG-3<sup>0</sup> ;

Human KIM-1, forward 5<sup>0</sup> -CTGCAGGGAGCAATAAGGAG-3<sup>0</sup> , reverse 5<sup>0</sup> -TCCAAAGGCCATCTGAAGAC-3<sup>0</sup> ; Human Nox4, forward 5<sup>0</sup> -GGATCACAGAAGGTCCCTAGCAG-3<sup>0</sup> , reverse 5<sup>0</sup> -GCGGCTACATGCACACCTGAGAA-3<sup>0</sup> ; Human Nox2, forward 50–30TTCCAGTGCGTGTTGCTCGAC, reverse 50–30GATGGCGGTGTGCAGTGCTAT; Human β-actin, forward 5<sup>0</sup> -CGCCGCCAGCTCACCATG-3<sup>0</sup> , reverse 5<sup>0</sup> -CACGATGGAGGGGAAGACGG-3<sup>0</sup> ; Mouse IL-6, forward 5<sup>0</sup> -GAGGATACCACTCCCAACAGACC-3<sup>0</sup> , reverse 5<sup>0</sup> -AAGTGCATCATCGTTGTTCATACA-3<sup>0</sup> ; Mouse TNF-α, forward 5<sup>0</sup> - CATCTTCTCAAAATTCGAGTGACAA-3<sup>0</sup> , reverse 5<sup>0</sup> -TGGGAGTAGACAAGGTACAACCC-3<sup>0</sup> ; Mouse MCP-1, forward 5<sup>0</sup> - CTTCTGGGCCTGCTGTTCA-3<sup>0</sup> , reverse 5<sup>0</sup> -CCAGCCTACTCATTGGGATCA-3<sup>0</sup> ; Mouse KIM-1, forward 5<sup>0</sup> -CAGGGAAGCCGCAGAAAA-3<sup>0</sup> , reverse 5<sup>0</sup> -GAGACACGGAAGGCAACCAC-3<sup>0</sup> ; Mouse β-actin, forward 5<sup>0</sup> - CATTGCTGACAGGATGCAGAA-3<sup>0</sup> , reverse 5<sup>0</sup> -ATGGTGCTAGGAGCCAGAGC-3<sup>0</sup>

Assays were run over 40 cycles with the following conditions: denaturation at 95◦C for 20 s, annealing at 58◦C for 20 s, and elongation at 72◦C for 20 s. β-actin was used to normalize the expression values of the other genes.

### Western Blot Analysis

Protein was isolated from pulverized tissue or cells from 6 well plates in ice-cold RIPA-Buffer (Beyotime, Jiangsu, China). BCA protein quantitative kit (Beyotime, Jiangsu, China) was

used to evaluate the protein concentration. For Western blots, total protein were loaded in 10% SDS-PAGE and transferred onto nitrocellulose membranes. After blocking, membranes were incubated with rabbit anti-Nox2, Nox4, anti-KIM-1, anti-RIP1, anti-RIP3, anti-P-MLKL, anti-cleaved-caspase-3 antibody, anti-cleaved-caspase-8, anti-cleaved-caspase-12, and mouse antiβ-actin antibody for 18 h at 4◦C, then incubated with IRDye 800 conjugated secondary antibody for 1.5 h at room temperature (1:10000, Rockland immunochemicals, Gilbertsville, PA, USA). Images were detected by Li-Cor/Odyssey infrared image system (LI-COR Biosciences, Lincoln, NE, USA) and quantified using the Image J software (NIH, Bethesda, MD, USA).

### Flow Cytometry

The extent of programmed cell death was detected by flow cytometry (BD FACSVerse, USA) using AV-FITC/PI apoptosis detection kit (Bestbio, Shanghai, China). Briefly, HK2 cells were harvested and washed twice with PBS after incubation in cell culture bottle with/without cisplatin and PA for 24 h. The cells were centrifuged at 1500 rpm/min for 5 min and stained with 5 µL Annexin V-FITC and 10 µL PI in the dark, followed by flow cytometry and quantified using FlowJo 7.6 software.

### Determination of MDA and GSH

The levels of MDA and GSH in cell or in mouse tissues were measured with a commercial kit (Jiancheng Co., Nanjing, China) according to the manufacturer's instructions. Thiobarbituric acid reacts with MDA, degradation product of lipid peroxidation in, to generate red compound which has maximum absorbance at 532 nm. This assay is commonly called thiobarbituric acid reactive substances assay (TBARS assay). Homogenate was centrifuged at 4000 rpm/min (10 min) after incubation with TBA reagent for 40 min at 100◦C. Supernatant was measured at 530 nm. 5,5-dithiobis-2-nitrobenzoic acid (DTNB) reacts with sulfhydryl compounds to generate a yellow compound, whose absorption peaks at 405 nm. GSH concentration was determined based on the absorbance of yellow compound. The homogenate was obtained and centrifuged at 3500 rpm/min for 10 min. The supernatant was reacted with DTNB and generated yellowcolored complex, which was measured at 405 nm.

### DCF Assay

DCF, the oxidation product of 2,7-dichlorodihydro-fluorescein diacetate, is a marker of cellular oxidation. Cisplatin-induced the generation of ROS as revealed by increased 2,7-dichlorodihydrofluorescein, which was measured by fluorescence microscopy with excitation of 488 nm and emission of 525 nm after cells were incubated with DCF (10 µL/L) for 20 min at 37◦C in no FBS-containing DMEM/F12 medium.

### DHE Staining

The oxidation of DHE, ethidium bound to DNA and fluoresced red, was used to estimate intracellular ROS levels. Cells were incubated with 5 µM freshly prepared DHE solution (Beyotime, Jiangsu, China) for 30 min at 37◦C and then measured under fluorescence microscopy.

## Knockdown of Nox2 and Nox4 in HK2 Cells

Nox2 and Nox4 were silenced by transfection with sequencespecific or non-targeting siRNA (GenePharm, Shanghai, China)/shRNA (GeneChem Co. Ltd., Shanghai, China) using LipofectamineTM 2000 reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions. Briefly, for each well of a 6-well plate, 5 µl siRNA/shRNA and 5 µl lipofectamine 2000 were diluted in 200 µl Opti-DMEM separately, and incubated for 5 min at room temperature in the dark. The diluted siRNA/shRNA was combined with the diluted Lipofectamine 2000 and incubated for 20 min at room temperature. Finally, the mixture with 100 ml Opti-DMEM was applied to the cells. After incubation for 12 h, the medium was replaced with fresh DMEM supplemented with 5% FBS. Cells with Nox4 shRNA screened by puromycin were cultured in incubator at 37◦C. The transfection efficiency was then evaluated using Western blotting.

### Renal Histology and Immunohistochemistry

Renal tissues were fixed in 4% PFA immediately. After dehydration, samples were embedded in paraffin. According to the manufacturer's instruction, PAS staining was performed in Paraffin sections (3–5 µm) to assess the degree of tubulointerstitial damage and examined by light microscope (Olympus, Japan) at 200 × magnification. On PAS-stained kidney sections (n = 6–8), kidney damage in the cortical proximal was scored as the approximate extent of tubules that displayed tubular necrosis, cast formation, and tubular dilation as follows: 0 = normal; 1 = 10%; 2 = 10–25%; 3 = 26–50%; 4 = 51–75%; 5 = 75–95%; 6 = more than 96%. Immunohistochemistry was performed in paraffin sections using a microwave-based antigen retrieval technique (Li et al., 2012). Sections were incubated with rabbit anti-Nox4, anti-KIM-1, anti-TNF-α, and rabbit anti-F4/80 antibody overnight at 4◦C. After incubation in secondary antibody and chromagen liquid DAB (3,30-diaminobenzidine tetrahydrochloride), the slides were counterstained with hematoxylin. The results were analyzed by Image Analysis System (AxioVision 4, Carl Zeiss, Jena, Germany).

### Statistical Analyses

Data are expressed as the mean ± SEM. Statistical significance was analyzed by two-tailed unpaired t-test or one-way analysis of variance (ANOVA), followed by Tukey post hoc tests using GraphPad Prism 5 software.

### DISCUSSION

Nephrotoxicity leads to high mortality in cisplatin-treated patients with cancer, therefore identification of preventive agents for cisplatin-induced AKI is needed for clinical treatments. Here, we identified a novel TCM monomer, protocatechuic aldehyde, that protects against cisplatin-induced injury both in vivo and in vitro via inhibiting Nox-mediated oxidative stress,

renal inflammation, and programmed cell death of renal tubular epithelial cells.

Protocatechuic aldehyde is a phenolic acid compound isolated from several types of Chinese herbs, including roots of miltiorrhiza, leaves of Stenolomachusanum (L.) Ching and Ilex chinensis Sims (Li et al., 2012). In several animal models, including cerebral ischemia model (Guo et al., 2016) and experimental model of sepsis (Xu et al., 2012). PA has shown pharmacological effects on inflammation and oxidative stress, which are highly correlated with pathogenesis of cisplatin nephropathy (Gao et al., 2011; Wei et al., 2013). For example, PA alleviated cerebral ischemia-reperfusion-induced oxidative injury via activating protein kinase Cε/Nrf2/HO-1 pathway (Guo et al., 2016). In addition, treatment of PA suppressed TNF-α-induced NF-κB phosphorylation and HMGB1 expression, attenuating inflammatory response in RAW264.7 cells (Xu et al., 2012). PA was also reported to reduce oxidative stress in sh-sy5y cells by targeting Dj-1. However, whether PA prevents cisplatin-induced AKI and the underlying mechanisms are still unknown.

In the current study, we found PA inhibited cisplatin-induced decline of renal function, renal damage, cell death and apoptosis, and inflammatory response in a dosage-dependent manner. Interestingly, accumulating evidence shows that PA possesses anti-cancer property by targeting cyclin D1-regulated cell cycle of tumor cells (Jeong and Lee, 2013; Choi et al., 2014). Here, we used three solid tumor cells lines, SMCC-7721, BEL-7402, and U87, to test the effect of PA on tumor cell viability. We found PA failed to reinforce the tumor suppressive effect of cisplatin. However, data show that PA gently promoted anti-tumor activity of cisplatin in U87 cell lines 24 h after cisplatin incubation. PA didn't reduce anti-tumor effects of cisplatin when protected against cisplatin nephropathy, indicating PA-based therapy for cisplatin-induced nephrotoxicity in cancer patients may be effective and promising.

Our results also show that PA significantly attenuated oxidative stress by reducing ROS production in cisplatinchallenged tubular epithelial cells and kidney tissues. This may be one of the most critical mechanisms by which PA attenuates cisplatin nephropathy. To date, evidence has shown that oxidative stress plays critical roles in pathogenesis of cisplatininduced nephrotoxicity (Humanes et al., 2012). Multiple sources for ROS in cells have been identified, including mitochondria, xanthine oxidase, cytochrome P-450, and uncoupled nitric oxide synthase (Humanes et al., 2012). NAPDH oxidases, a set of membrane-associated proteins using NADPH to transfer electrons across biological membranes and therefore generating ROS, are regarded as one of key sources for ROS in the kidney (Schreck and O'Connor, 2011). The Nox family consists of seven members (Noxs 1–5, Duox1 and 2). Nox2 and Nox4 are highly expressed within the kidney and Nox4 is known as the predominant form, which plays important roles in renal oxidative stress and kidney injury (Gill and Wilcox, 2006; Sedeek et al., 2013; Kim et al., 2016; Oh et al., 2016). Our data show that lentivirus-mediated knockdown of Nox4 in vivo significantly attenuated cisplatin nephropathy. In the present study, both in vivo and in vitro studies show that PA blocked cisplatininduced ROS generation, which was evidenced by results of GSH assay, MDA assay, and ROS assay. Of note, cisplatin significantly increased protein levels of Nox4 both in cisplatin-challenged HK2 cells and cisplatin nephropathy, which was largely blocked by PA treatment. More importantly, in vitro data indicated that silencing Nox4, partly blocked the inhibitory effects of PA on cisplatin-upregulated levels of KIM1, cleaved-caspase-3 and production of inflammatory factors. This demonstrates that Nox4 may be the primary target in mediating the anti-oxidative stress and protective role of PA in cisplatin-induced nephrotoxicity. Additionally, it is noteworthy that Nox2 was significantly induced by cisplatin, but suppressed by PA treatment, both in vivo and in vitro; however, results generated from cisplatin-stimulated Nox2 knockdown HK2 cells showed the less important role of Nox2 in mediating the effects of PA compared with Nox4, as a critical enzyme in the injury of kidney and other organs, (Gill and Wilcox, 2006). Whether Nox2 is involved in PA treatment needs to be further validated.

We also found that PA diminished cisplatin-induced programmed cell death, especially necroptosis. As the best characterized regulated necrosis, necroptosis is known to play a critical role in cisplatin-induced AKI (Xu et al., 2015; Linkermann, 2016). Administration of Nec-1, an inhibitor for necroptosis, alleviated acute kidney injury induced by cisplatin (Tristão et al., 2012), cyclosporin A (Ouyang et al., 2012), and ischemia-reperfusion injury (Zhang et al., 2013). Emerging evidence shows that RIP1, RIP3, and downstream MLKL serve as predominant regulators in necroptosis while genetic deletion or pharmalogical inhibition of these key genes significantly reduced cisplatin-induced AKI (Linkermann et al., 2013, 2014; Xu et al., 2015). Compared with apoptosis, necroptosis plays more pivotal roles in the induction of inflammatory responses. Cell membrane collapse induces the release of damage-associated molecular patterns (DAMPs) including high-mobility group box 1 (HMGB1), heat-shock proteins, uric acid, and IL-33, which may interact with receptors (such as Toll-like receptors) and initiate downstream signaling pathways to enhance renal inflammation in a positive-feedback loop (Scaffidi et al., 2002; Vilaysane et al., 2010; He et al., 2011). However, the mechanisms by which PA blocks necroptosis is still to be determined. PA may direct target key mediators in necroptosis-regulated pathways or through indirect mechanisms in which Nox-dependent oxidative stressare invovled.

Collectively, as shown in **Figure 11** our study demonstrated that PA substantially alleviated the decline of renal function and renal damage while preventing renal oxidative stress, necroptosis, and consequent inflammatory response. This may be correlated with the inhibitory effect of PA on Nox4. Considering the finding that PA significantly protects against cisplatin-induced acute kidney injury without compromising the anti-cancer properties of cisplatin, it should be further explored as a preventive agent for cisplatin-treated cancer patients.

### AUTHOR CONTRIBUTIONS

X-MM and JL designed experiments and took part in the critical revision of manuscript; LG and W-FW carried out experiments and participated in drafting of manuscript; G-LR and H-DL

provided a series of experimental instructions and help; QY and X-FL conducted primary screening test of relevant drugs; TX, ZL, and B-MW analyzed experimental results; T-TM, CH, and LD analyzed and interpretation of data; YH, LZ, and XL assisted with experiments on animals.

### FUNDING

This study was supported by the National Natural Science Foundation of China (No. 81300580, No.

### REFERENCES


81570623) and by Science and Technological Fund of Anhui Province for Outstanding Youth of China (Grant number: 1608085J07) and the Grants for Scientific Research of BSKY (XJ201307) from Anhui Medical University.

### ACKNOWLEDGMENTS

We thank Dr. Austin Cape at ASJ Editors for careful reading and feedback.



**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2016 Gao, Wu, Dong, Ren, Li, Yang, Li, Xu, Li, Wu, Ma, Huang, Huang, Zhang, Lv, Li and Meng. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

## Assessing the synergy between cholinomimetics and memantine as augmentation therapy in cognitive impairment in schizophrenia. A virtual human patient trial using quantitative systems pharmacology

#### Hugo Geerts 1, 2 \*, Patrick Roberts <sup>3</sup> and Athan Spiros <sup>1</sup>

#### Edited by:

*Thomas J. Anastasio, University of Illinois at Urbana-Champaign, USA*

#### Reviewed by:

*Alexander Dityatev, German Center for Neurodegenerative Diseases (DZNE), Germany Filippo Caraci, University of Catania, Italy*

#### \*Correspondence:

*Hugo Geerts, In Silico Biosciences, 686 Westwind Dr, Berwyn, PA 19312, USA hugo-geerts@ in-silico-biosciences.com*

#### Specialty section:

*This article was submitted to Experimental Pharmacology and Drug Discovery, a section of the journal Frontiers in Pharmacology*

> Received: *23 June 2015* Accepted: *31 August 2015* Published: *22 September 2015*

#### Citation:

*Geerts H, Roberts P and Spiros A (2015) Assessing the synergy between cholinomimetics and memantine as augmentation therapy in cognitive impairment in schizophrenia. A virtual human patient trial using quantitative systems pharmacology. Front. Pharmacol. 6:198. doi: 10.3389/fphar.2015.00198* *1 In Silico Biosciences, Berwyn, PA, USA, <sup>2</sup> Perelman School of Medicine, University of Pennsylvania, Philadelphia, PA, USA, <sup>3</sup> Department of Veterinary and Comparative Anatomy, Pharmacology and Physiology, Washington State University, Pullman, WA, USA*

While many drug discovery research programs aim to develop highly selective clinical candidates, their clinical success is limited because of the complex non-linear interactions of human brain neuronal circuits. Therefore, a rational approach for identifying appropriate synergistic multipharmacology and validating optimal target combinations is desperately needed. A mechanism-based Quantitative Systems Pharmacology (QSP) computer-based modeling platform that combines biophysically realistic preclinical neurophysiology and neuropharmacology with clinical information is a possible solution. This paper reports the application of such a model for Cognitive Impairment In Schizophrenia (CIAS), where the cholinomimetics galantamine and donepezil are combined with memantine and with different antipsychotics and smoking in a virtual human patient experiment. The results suggest that cholinomimetics added to antipsychotics have a modest effect on cognition in CIAS in non-smoking patients with haloperidol and risperidone and to a lesser extent with olanzapine and aripiprazole. Smoking reduces the effect of cholinomimetics with aripiprazole and olanzapine, but enhances the effect in haloperidol and risperidone. Adding memantine to antipsychotics improves cognition except with quetiapine, an effect enhanced with smoking. Combining cholinomimetics, antipsychotics and memantine in general shows an additive effect, except for a negative interaction with aripiprazole and quetiapine and a synergistic effect with olanzapine and haloperidol in non-smokers and haloperidol in smokers. The complex interaction of cholinomimetics with memantine, antipsychotics and smoking can be quantitatively studied using mechanism-based advanced computer modeling. QSP modeling of virtual human patients can possibly generate useful insights on the non-linear interactions of multipharmacology drugs and support complex CNS R&D projects in cognition in search of synergistic polypharmacy.

Keywords: cognition, polypharmacy, antipsychotics, cholinomimetic, schizophrenia

## Introduction

While polypharmacy is more of a rule than an exception in reallife clinical treatment, preclinical animal models are ill-equipped to address the issue of comedication because of fundamental species-specific differences in drug metabolism, the impact of human-specific genotypes, the different pharmacology of the drugs for human vs. rodent targets and incomplete pathology (for a review see, Geerts, 2009).

Cognitive Impairment in Schizophrenia (CIAS) is a major unmet medical need for this patient population; while psychosis can be readily managed by the current antipsychotic drug armentarium, cognitive and negative symptoms are hampering patients to return to a more normal professional life (Kitchen et al., 2012). This has prompted the major stakeholders from industry, regulatory agencies and academia to develop a regulatory path for cognitive enhancement, resulting in the development of the Matrics battery (Green and Nuechterlein, 2004). Over the last 15 years, many novel highly selective drugs have been tested for cognitive enhancement as augmentation therapy without much success (Dunlop and Brandon, 2015).

Symptomatic treatment has been successful in Alzheimer's disease with cholinomimetics and memantine. Galantamine is an acetylcholinesterase inhibitor (AChE-I) with allosteric potentiating effects on nicotinic receptors (nAChR) currently approved for Alzheimer's Disease (Tariot et al., 2000) and has been tested for cognitive improvement in schizophrenia (Norén et al., 2006; Deutsch et al., 2008; Dyer et al., 2008; Sacco et al., 2008). Donepezil, a pure AChE-I has been studied for cognitive improvement in schizophrenia with mixed results (Akhondzadeh et al., 2008; Keefe et al., 2008; Gauthier and Molinuevo, 2013). On the other hand, memantine is a weak NMDA-antagonist approved for moderate-to-severe AD (Reisberg et al., 2003) with preferential affinity against the excitatory-inhibitory synapses in cortical networks (Kotermanski and Johnson, 2009). The effect of these compounds in CIAS has yielded equivocal results. Possible reasons for the lack of clear results include nontrivial pharmacodynamic relationships between the investigative drugs and the different baseline antipsychotics. For instance, some antipsychotics such as olanzapine and clozapine affect the muscarinic cholinergic receptors; indeed olanzapine is a documented antagonist of M1, M2, M3, M4, and M<sup>5</sup> mAChR with respective affinities of 26, 18, 52, 17, and 20 nM (Bymaster et al., 1996). Therefore, it is to be expected that the dose-response of cholinomimetics in CIAS might be highly dependent upon the different antipsychotics.

To complicate the situation even further, anticholinergics are often used in treating the symptoms of extrapyramidal motor symptoms and can affect the cognitive function substantially (Ogino et al., 2014). In addition, a disproportionate fraction of schizophrenia patients tend to smoke (Dalack et al., 1998), adding more complexity to the pharmacodynamic interaction of cholinomimetics at the level of nicotinic receptors.

A recent paper suggested that the combination of memantine and AChE inhibitors, in particular galantamine, would show a bigger effect in CIAS because of their complementary pharmacology on pyramidal cells and interneurons (Koola et al., 2014).

Testing such a combination therapy is likely beyond the capability of preclinical animal testing due to the complexity and combinatorial challenges of the trial design. To explore the possible clinical applications of this combination therapy, other approaches need to be explored. In this paper we propose an advanced version of a computer-based Quantitative Systems Pharmacology (QSP) platform, a mechanism-based computer model of the relevant humanized cortical networks that has been developed for clinical readouts in psychiatry and neurology, and calibrated with group average clinical data. Such an approach is similar to Physiology-Based PharmacoDynamic (PBPD) Modeling in line with new terminology around Physiology-Based Pharmacokinetic Modeling (PBPK) and is gaining traction in pharmaceutical research and development. The platform has been able to blindly, prospectively and correctly predict an unexpected clinical outcome in schizophrenia and Alzheimer's disease (AD) (Geerts et al., 2012; Nicholas et al., 2013; Liu et al., 2014) and has been calibrated for clinical cognitive outcomes in conditions of chronic schizophrenia (Geerts et al., 2013).

Testing the platform in a number of practical clinical situations with known outcomes, such as the augmentation therapy of cholinomimetics and memantine on antipsychotics is mandatory to help constrain the platform. Once calibrated and constrained, with human clinical data, this QSP approach can then be used in rationally designed multi-target drug discovery programs.

### Methods

We use a previously described (Geerts et al., 2013) biophysically realistic and mechanism-based QSP platform to simulate the impact of augmentation therapy with cognitive enhancers in virtual schizophrenia patients. This QSP platform has been calibrated against observed clinical effects on the N-back working memory test with various therapeutic interventions in diverse patient populations and recapitulates the negative pharmacodynamic clinical effect of risperidone augmentation on clozapine. Basically, the platform consists of a receptor competition model that allows accurate quantification of drug target exposure, a biophysically realistic neuronal network that captures the microarchitecture of a cortical column and a calibration module that relates computer model outcome to

**Abbreviations:** 5HTTLPR, Serotonin transporter promotor: long vs. short isoform; ACh, acetylcholine; AChE, acetylcholinesterase; AChE-I, acetylcholinesterase inhibitors; CIAS, Cognitive Impairment in Schizophrenia; COMT, Catechol-O-Methyl Transferase; CYTP450, Cytochrome P450 enzyme; DA, Dopamine; dlPFC, dorso-lateral Prefrontal Cortex; EPS, Extra-pyramidal symptoms; GABA, gamma-amino butyric acid; GPCR, G-protein coupled receptors; mAChR, muscarinic acetylcholine receptor; MATRICS, Measurement And Treatment Research to Improve Cognition in; Schizophrenia; nAChR, nicotinic acetylcholine receptor; NMDA, N-methyl-D-aspartate (glutamate receptor subtype); PBPD, Physiology-based pharmacodynamic modeling; PBPK, Physiology-based pharmacokinetic modeling; PDSP, Psycho-active Drug Screening Program; PET, Positron Emission Tomography; PFC, Prefrontal Cortex; PK/PK, pharmacokinetic-pharmacokinetic; QSP, Quantitative Systems Pharmacology; RBANS, Repeatable Battery for the Assessment of Neuropsychological Status; R&D, (Pharmaceutical) Research and Development.

actual clinical results. The platform includes the neurophysiology of over 30 CNS targets, ranging from catecholamine GPCR over various glutamate, GABA receptors and ligand and voltage-gated ion channels to enzymes such as Catechol-O-methyl transferase (COMT) and PDE10 and neurotransmitter transporters.

### Defining Target Exposure in Quantitative Systems Pharmacology Model

The receptor competition model (Spiros et al., 2010; Spiros and Geerts, 2012) calculates the degree of activation of various postsynaptic receptors (dopamine, serotonin, norepinephrine, and cholinergic neurotransmitters) in the presence of antipsychotics. The affinity of the parent molecule and its major metabolite for both pre- and post-synaptic receptors, derived from the Psychoactive Drug Screening Program (Besnard et al., 2012) is used to calculate the competition with endogenous neurotransmitters. The presynaptic autoreceptor neurophysiology properties are calculated from preclinical data using fast-cyclic voltammetry constrained with clinical imaging data (Nicholas et al., 2013). The functional intrasynaptic concentration of the specific antipsychotic is determined from calculating the concentration that corresponds to the clinically observed displacement of a radio-active D2R specific PET tracer, such as raclopride (Spiros et al., 2012). We assume regular clinical doses 400 and 600 mg, i.e., 6 mg risperidone, 10 mg haloperidol, 15 mg olanzapine, 200 quetiapine, and 30 mg aripiprazole.

### Quantitative Systems Pharmacology Model for Cognitive Impairment in Schizophrenia

The QSP model consists of a network of 80 four-compartment pyramidal and 40 two-compartment interneurons with the effects of dopaminergic, serotonergic, noradrenergic and cholinergic modulation (including a spatio-temporal receptor state model for allosteric modulation of the different nicotinic ACh receptors) and has been described in detail elsewhere (Geerts et al., 2013). A subset of pyramidal cells is stimulated for a very short period at a certain time point, reflecting a sensory stimulus typical of a working memory paradigm. The actual membrane potential of each compartment can be calculated from the actual conductances, which are dependent upon the activation level of various G-protein coupled receptors. The resulting state diagram shows a synchronized firing of the cells during 5–10 s after they have been stimulated for a short period (100 ms). While this computational neuroscience model has been designed using in vivo electrophysiological single-unit recordings in non-human primates (Williams and Goldman-Rakic, 1995) performing a working memory task and therefore probably only reflects the maintenance phase, the outcome could be generalized to the strength of a memory trace (Roberts et al., 2012; Geerts et al., 2013). We have shown previously that the duration of this synchronized firing correlates well with actual 2-Back working memory task in a variety of experimental interventions in humans (Geerts et al., 2013).

Schizophrenia pathology is implemented using insights from human neuroimaging, genetic and neuropathology data and includes a hypodopaminergic cortical D1R tone (Durstewitz and Seamans, 2008), NMDA-R hypofunction (Coyle, 2006) documented by a hypocortical-hyperstriatal imbalance in metabolic imaging (Meyer-Lindenberg et al., 2002), a GABA deficit (Volk and Lewis, 2002) applied here to the network interneurons, and a noisier background signal (Winterer et al., 2000), resulting in a clinical cognitive deficit which is dependent upon the cognitive domain, but on average is 1.5 standard deviations lower than healthy controls (Saykin et al., 1994). The pathology in the computer model leads to a similar deficit between a "healthy environment" and the schizophrenia condition.

### Implementation of Pharmacology for Cognitive Enhancers

Donepezil is an AChE-inhibitor with a K<sup>i</sup> of 20 nM while galantamine inhibits AChE-I with a much lower affinity of 800 nM and in addition weakly and allosteric potentiates α<sup>7</sup> and α4β<sup>2</sup> nAChR (Woodruff-Pak et al., 2002). Imaging studies with <sup>11</sup>C-PMP have suggested that 10 mg donepezil and 24 mg galantamine lead to brain AChE-inhibition levels of 35% (Shinotoh et al., 2001; Darreh-Shori et al., 2008). These clinically observed inhibition levels can be used to calculate the daily dose to affect 50% brain AChE-inhibition, which corresponds to 18.5 mg for donepezil and 44.5 mg for galantamine, resulting in inhibition levels of 20% for 5 mg donepezil, 15% for 8 mg galantamine and 24% for 16 mg galantamine. ACh half-life, T, in the cholinergic receptor competition model is then calculated as T0/(1-Enzyme inhibition), with T<sup>0</sup> being the half-life in untreated patients. The AchE is one of the fastest enzymes in the human body (Iwanaga et al., 1994), leading to a half-life in the untreated situation of 5 ms. This leads to ACh half-lives of 6.9 and 7.7 ms for donepezil at 5 and 10 mg and to half-lives of 5.9, 6.8, and 7.7 ms for galantamine at 8, 16, and 24 mg.

In addition, galantamine has a small allosteric potentiating effect on nAChR (Woodruff-Pak et al., 2002), which we implemented as a 5, 10, or 15% (respectively for 8, 16, and 24 mg) relative increase in both α<sup>7</sup> nAChR and α4β<sup>2</sup> nAChR activation levels.

### Implementation of Smoking

As a disproportionally large fraction of schizophrenia patients smoke (Dalack et al., 1998), we implement the effect of nicotine on both α4β<sup>2</sup> nAChR and α<sup>7</sup> nAChR. Nicotine has a much higher affinity for α4β<sup>2</sup> nAChR than for a7 nAChR and imaging studies with the PET radiotracer 18F-2-Fluoro-A85380 showed an almost complete saturation of α4β<sup>2</sup> nAChR in smokers (Brody et al., 2006). We assume an increase in α4β<sup>2</sup> nAChR activation of 20% as the receptors are already naturally active. However, this level of α4β<sup>2</sup> nAChR activation, together with the continuous nicotine exposure likely overall leads to receptor desensitization (Grady et al., 2012). Because α4β<sup>2</sup> nAChR regulates GABA release (McClure-Begley et al., 2009; Zappettini et al., 2011) we implement the desenitization induced by the smoking condition as a two-fold decrease in GABA conductances, leading to a greater firing of the network. Given the relative much lower affinity of nicotine for the α<sup>7</sup> nAChR (20,000 nM vs. 100 nM) (Buisson et al., 1996), we assume smoking does not affect α<sup>7</sup> nAChR. Note that the amount of ACh bound to α4β<sup>2</sup> (and of α7) nAChR is further determined by the galantamine or donepezil

mediated AChE inhibition in addition to inhibition of the presynaptic M2 mAchR autoreceptor by specific antipsychotics, such as olanzapine. In our model, this is illustrated by the fact that binding of ACh to α4β<sup>2</sup> nAChR ranges from 24% (non-smoking patient on haloperidol) to 62% (olanzapine in smoking patients on 24 mg galantamine).

### Implementation of Memantine Pharmacology

Memantine is a relatively weak NMDA-R inhibitor that has a larger affinity for the NMDA-NR2C/2D subunit (Kotermanski and Johnson, 2009) in physiological conditions. Based upon the observation that the NR2C/2D subunits are preferentially located on inhibitory interneurons (Monyer et al., 1994) in rats, memantine's pharmacology is implemented using a two-fold greater inhibition of the NMDAR on interneurons as compared to the NMDAR on pyramidal cells.

Data further suggest that the functional memantine concentration in the human brain is relatively small; in our earlier paper on the cognitive model (Roberts et al., 2012) for Alzheimer's Disease, a 1% decrease in gNMDA on interneurons resulted in a selective positive impact on moderate to severe AD cases but not in the situation of mild cases, suggesting that such an inhibition corresponded to the clinical dose of 20 mg.

### Implementation of Pharmacological Profile of Antipsychotics

The affinity parameters for each individual drug and neurotransmitter for human receptors were derived from the standardized PDSP database (http://pdsp.med.unc.edu/ indexR.html) (Besnard et al., 2012). Importantly, the active moiety of antipsychotics, taking into account the pharmacology of metabolites was used (see **Figure 1**). We took great care in determining the functional intrasynaptic concentration of the various antipsychotics using published <sup>11</sup>C-raclopride displacements observed with specific antipsychotic dose combinations using the receptor competition model, described above.

### Results

### Augmentation Therapy with Memantine

We studied the effect of increasing memantine doses on the performance of the in silico network for CIAS in the presence of antipsychotics. **Figures 2A,B** shows the effect of memantine on the estimated 2-back working memory outcomes, respectively in the absence and presence of nicotine. In the absence of smoking, with the exception of quetiapine, all drugs improve cognitive readout with increasing memantine doses with the greatest effect observed for aripiprazole (from 69 accuracy to 81% in a 2-back test). For smoking conditions, the increase in frequency of accurate responses as a function of the optimal memantine dose is amplified for Risperidone (maximal increase from 5 to 14%), Haloperidol (maximal effect increases from 5 to 12%) and olanzapine (maximal effect from 6 to 11%), but not for aripiprazole and quetiapine. Note that this maximal effect happens at memantine doses of 40 mg. At clinically relevant doses of 20 mg, the effect is about half as much.

### Augmentation Therapy with AChE-I

We then simulated the effect of augmentation strategy with ACh inhibitors added to antipsychotics on cognitive outcome (**Figure 3**), both in the absence and presence of nicotine. In non-smokers, AChE-I dose-dependently improved cognitive outcomes with risperidone, aripiprazole (donepezil only) and haloperidol, with only the 24 mg galantamine showing a robust improvement of >10% in correct responses. In the presence of olanzapine and aripiprazole with galantamine, a tendency was observed for an inverse U-shape dose-response. In the presence of quetiapine, increased AChE-I worsened responses.

Smoking tended to increase the procognitive effect as a stand alone (i.e., without cholinomimetics) and in the presence of risperidone, haloperidol and aripiprazole (donepezil only). For instance the fraction of correct responses went up from 66 to 71% for risperidone with 24 mg galantamine. However, smoking also tended to diminish the cognitive response when augmented with quetiapine and olanzapine. When aripiprazole is augmented with galantamine, smoking shifted the peak response of the inverse U-shaped curve.

### Augmentation Therapy with a Combination of Memantine and AChE-I

This section deals with the simulation of combination therapy of memantine with cholinomimetics on cognitive outcomes. Note that we simulate the outcome of up to five agents in the same patient, i.e., the parent molecule and active metabolite of an antipsychotic, an AChE-I such as donepezil and galantamine, memantine and nicotine. This is a situation that is often encountered in clinical practice.

From the data a complex picture emerges, ranging from a negative effect (adding AChE-I lowers the effect of memantine) with quetiapine and aripiprazole, to an additive effect with risperidone (see **Table 1**). Synergism is clearly observed with olanzapine and in some cases with haloperidol in the nonsmoking case and with haloperidol in the smoking condition.

**Figure 4** shows a dose-response of memantine on cognitive effects in the presence of olanzapine, and with donepezil or galantamine. Memantine increases cognitive effects dosedependently in the absence of the AChE-I and the slope (in % correct on the 2-back WM test for 0–40 mg of memantine) of 0.08 is increased to 0.23 and 0.34 when adding donepezil 5 and 10 mg, respectively. When adding galantamine, the slope is increased to 0.13, 0.30, and 0.16 for 8, 16, and 24 mg, respectively. This suggests a synergistic effect with increasing concentrations of donepezil and galantamine, although at the highest galantamine dose of 24 mg the effect is somewhat attenuated. This synergistic effect with olanzapine however disappears in the smoking conditions.

An interesting pattern emerges with regard to different doses of quetiapine. In the non-smoking condition, memantine has a clear beneficial effect on cognition in the absence of any AchE-I with a tendency for greater effect at greater quetiapine doses. Adding donepezil or galantamine reverses this effect at all quetiapine (200–600 mg daily) doses into a negative interaction, i.e., cognitive performance drops with increasing memantine dose. The negative interaction decreases slightly as the quetiapine

clear that many antipsychotics have complex pharmacologies leading to non-linear interactions in the human networks. (A) pharmacology of active moiety of aripiprazole, haldoperidol, and olanzapine, (B) pharmacology of risperidone and quetiapine.

dose increases. This is likely because adding AchE-I to quetiapine substantially improves cognition; so that further block of the NMDA receptors with memantine reduces GABA tone and drives the network to fire at a very high frequency determined only by the refractory period of the pyramidal neurons, reducing the variability of the interspike distribution and the information bandwidth.

In smoking conditions, the slope of cognitive improvement with memantine without AchE-I is much smaller for all quetiapine doses, because the baseline performance of smoking and quetiapine is already higher than in the non-smoking condition. Adding AchE-I in the smoking condition to quetiapine and memantine has a relative smaller negative effect than in the non-smoking condition, although the absolute values of the slopes are similar. This illustrates the complex non-linear pharmacodynamic interactions between different comedications.

**Table 1** shows the slopes of the memantine dose-response with donepezil and galantamine in the presence of specific antipsychotics and smoking/no-smoking conditions. As mentioned above, there is a synergistic effect in non-smokers on olanzapine and in smokers treated with haloperidol for both donepezil and galantamine.

### Discussion

This study addresses the question of a suggested synergistic interaction between AChE-I such as donepezil and galantamine and memantine on cognitive readouts in schizophrenia (Koola et al., 2014) by simulating real-life treatment combinations, including the effect of smoking. We approached this through a mechanism-based QSP computer modeling approach where the pathology of CIAS is combined with the pharmacology and target exposure of the respective drug combinations.

### Combination of AChE-I with Antipsychotics

The combination of AChE-I with antipsychotics shows a number of interesting and unexpected outcomes. The 24 mg galantamine has a higher response than any of the two doses of donepezil Geerts et al. Schizophrenia polypharmacy virtual patient trial

TABLE 1 | Slopes of the memantine dose-response for non-smokers and smokers, calculated from the trendline of the dose-responses in the different conditions and comparing the effect of adding 5 and 10 mg donepezil and 8, 16, and 24 mg galantamine to memantine in the presence of five different antipsychotics (RIS, risperidone; QUE, quetiapine; OLA, olanzapine; ARI, aripiprazole, HAL, haloperidol).


*The cases with synergy are noted in yellow, while the situation where addition of AChE-I to memantine worsens the dose-response is noted in red. The other cases suggest a pure additive effect. From the data a complex picture emerges, ranging from a negative effect (i.e., lower slopes) in the presence of quetiapine and aripiprazole (in both nonsmoking and smoking conditions) to an additive effect with olanzapine and quetiapine combined with donepezil in smoking conditions. Synergism is observed with olanzapine in the non-smoking case with a tendency for synergism in haloperidol non-smokers and with haloperidol in the smoking condition.*

for risperidone and haloperidol, an effect that is sustained in smokers. Interestingly galantamine shows an inverse U-shape dose-response with olanzapine in non-smokers; the lowest dose of 8 mg had a higher effect than the highest dose. Smoking also drives the responses of both donepezil and galantamine into an inverse dose-response when combined with aripiprazole or olanzapine. One could argue that the interaction of olanzapine with the presynaptic M<sup>2</sup> mAChR (a Ki of 18 nM) has a disproportionately larger impact given the affinity of ACh for the M<sup>2</sup> mAChR (which is in the 300 nM range). Blocking the presynaptic autoreceptor could further increase the release of Ach beyond the added effect of AChE inhibition which would then drive the postsynaptic nAChR into desensitization. The level of free ACh in schizophrenia patients is relatively normal, in contrast to the reduced free ACh in Alzheimer's patients, suggesting that the nAChR work on a very different baseline.

Aripiprazole and quetiapine both have strong 5-HT1A agonism and lack the inhibition of the presynaptic 5-HT1BR that is a hallmark of risperidone; these two properties might account for a relatively good baseline performance as compared to risperidone. The 5-HT1A pathway has been documented to play a role in antipsychotic response (Takekita et al., 2015), at least in negative symptoms. Therefore, the higher ACh tone resulting from inhibition of the AChE in the presence of quetiapine and aripiprazole has less dynamical range to improve.

A recent meta-analysis indeed suggest that antipsychotics do indeed have different effects on cognition (Désaméricq et al., 2014), with quetiapine having the greatest beneficial impact on global cognitive score, attention and speed of processing.

The effect of adding AChE-I to antipsychotics on cognition has been studied extensively in clinical trials with effect sizes in the range 0.4–0.6 (Ribeiz et al., 2010; Choi et al., 2013). Studies with donepezil have yielded controversial and mixed results, from positive results (Zhu et al., 2014) to negative results (Kohler et al., 2007; Keefe et al., 2008) with similar negative (Dyer et al., 2008; Lindenmayer and Khan, 2011) and positive results (Schubert et al., 2006; Buchanan et al., 2008) or no effect (Lee et al., 2007) for galantamine leading to the overall perception that these drugs only work marginally or not at all. Our modeling suggests that the nature of the antipsychotic and the condition of smoking does matter. In non-smokers, donepezil and galantamine work best with risperidone, haloperidol, aripiprazole, and olanzapine, but not with quetiapine. In smokers, both galantamine and donepezil enhance cognition with risperidone and haloperidol but have an inverse U-shape dose-response in olanzapine (lower AchE-I doses work best) while there is no effect in aripiprazole and quetiapine. Interestingly high-dose galantamine is inferior to placebo for a number of cognitive readouts in schizophrenia when added to antipsychotics (Dyer et al., 2008). Possible explanations for these observations include the interaction of olanzapine with the muscarinic receptors, especially the presynaptic M<sup>2</sup> mAChR autoreceptor. This would affect the amount of presynaptically released ACh and therfore interfere in a complex way with the increased half-life of ACh with AChE-I and the allosteric modulatory effect of galantamine. It is of interest to note that the trials with donepezil that showed some efficacy for donepezil, had risperidone and olanzapine as baseline medication (Akhondzadeh et al., 2008; Zhu et al., 2014) with a majority of patients on the lower donepezil dose of 5 mg.

Two trials with galantamine (Schubert et al., 2006; Lindenmayer and Khan, 2011) added to risperidone showed opposite effects. It is worthwhile to examine these two clinical trials in more detail. The positive Schubert trial was a short duration, 8 week trial (n = 16), almost all smokers with an average risperidone dose of 5.4 mg (with patients on risperidone for almost 3 years) where anticholinergics were excluded and showed a clear clinical and statistical benefit in delayed memory and attention on the RBANS scale. The negative Lindenmayer trial was a 52-week study of long-acting injection Risperidone (25 and 50 mg) in 32 patients in which galantamine was titrated up to 24 mg; no data on smoking were available. Anticholinergics were allowed for treatment of EPS side-effects, but no data are available on the frequency of this comedication.

Possible clinical trial design differences in these two studies leading to the opposite outcomes include the fraction of smokers (smoking tends to amplify galantamine's effect on cognition) and the somewhat lower dose of 25 mg Ris Consta in the Lindenmayer study as compared to 6 mg oral risperidone. Sensitivity analyses in the model show that lower risperidone doses tend to perform better in the cognitive model, leaving somewhat less dynamic

galantamine and donepezil show a dose-dependent improvement in the presence of risperidone, aripiprazole, and haloperidol, although the effect sizes differ. However, there is no improvement for quetiapine, probably due to the high baseline, and a more complex dose-response for olanzapine. Smoking tends to slightly enhance the responses of cognition in augmentation therapy with risperidone and haloperidol. However, smoking also tends to slightly decrease the responses of cognition in augmentation therapy with quetiapine and olanzapine. This is probably due to the many non-linear interactions between cholinergic modulation and the complex pharmacodynamics of antipsychotics.

range for additional effect of pro-cognitive enhancers; this is partially explained by the anti-cognitive effect of the dosedependent presynaptic 5-HT1B autoreceptor block that affects free 5-HT levels. Attention deficit has indeed been shown to be correlated with the amount of D2R blockade by risperidone (Uchida et al., 2009), at least in late-life schizophrenia; therefore higher doses of risperidone leading to a lower baseline provide a greater dynamic range for pro-cognitive effects of galantamine.

In general, however, the information available from peerreviewed articles usually does not have the granularity needed to identify the different comedications for individual patients in the different treatment arms. Exploration of other databases such as ADNI, where the comedications are given at the level of individual subjects or the database from electronic health records from the South London and Maudsley National Health Systems registry (Kadra et al., 2015) are a possibility to test the QSP model outcomes to real-life situations.

### Combination of Memantine with Antipsychotics

The simulations in this paper also suggest that memantine as augmentation therapy has a very modest effect on cognition that is enhanced in the smoking conditions. Again the best response is observed when adding memantine to risperidone and haloperidol, but at concentrations that are about twice as high as currently used in Alzheimer's disease. A meta-analysis of clinical trials (Kishi and Iwata, 2013) indeed suggests a modest effect of memantine on cognition in schizophrenia.

It is worthwhile to expand upon the unexpected clinical procognitive finding of a (weak) NMDA-antagonist. Memantine's interaction under in vivo conditions with Mg present, shows a higher affinity for the NMDA-NR2C subunit on the excitatoryinhibitory glutamatergic synapses as compared to the NMDA-NR2A/B subunits that are present on excitatory-excitatory synapses (Kotermanski et al., 2009). Smoking tends to desensitize the α4β<sup>2</sup> nAChR more than the α<sup>7</sup> nAChR; this would lead to a lower GABA tone as one of α4β<sup>2</sup> nAChR mediated processes regulates the GABA release (McClure-Begley et al., 2009; Zappettini et al., 2011). Memantine, through its antagonism on the excitatory-inhibitory glutamatergic synapse, also lowers the GABA-tone, opening the possibility for an additive or synergistic mechanism by addressing one of the major pathological hallmarks of the schizophrenia condition (Volk and Lewis, 2002).

### Combination of Memantine and AChE-I with Antipsychotics

When adding AChE-inhibitors to antipsychotics and memantine, the effect greatly depends upon both the nature of the antipsychotic and the smoking condition. Overall, the effect is additive, with the exception of aripiprazole and quetiapine for both smokers and non-smokers and for olanzapine in the smoking condition, where addition of AChE-I decreases the effect of memantine.

A synergistic effect is observed for olanzapine in the nonsmoking condition and for haloperidol in both conditions. The observation that adding AChE-I to memantine can improve cognitive readouts is in line with clinical data in Alzheimer's patients (Gauthier and Molinuevo, 2013). The inhibition on the AChE leads to slightly shorter half-life for AChE in the case of galantamine (5.9 ms for 8 mg galantamine vs. 6.5 ms for 5 mg donepezil). The observed lack of differentiation between donepezil and galantamine could be due to the fact that the impact of galantamine's allosteric modulation on nAChR has a quite limited effect (especially given the low efficacy) and is barely capable to compensate for the somewhat lower level of AChE-I. This is in line with clinical observations in the Alzheimer's field that suggest no detectable difference in treatment with either donepezil or galantamine (Tan et al., 2014).

The observation that olanzapine favors a synergistic effect between memantine and AChE-I points to an important role for the GABA interneurons in cortical circuits. Olanzapine has a substantial antagonism at the 5-HT3R that regulates GABA interneuron firing (Puig, 2004), a property it shares with clozapine and as such can compensate for the observed GABA dysfunction in schizophrenia (Volk and Lewis, 2002). Indeed, as mentioned above, memantine preferentially interacts with the NMDA receptor subunit on the excitatory-inhibitory synapses, increased levels of ACh through blocking of the AChE can activate α4β<sup>2</sup> nAChR that regulate GABA release and galantamine has an additional allosteric potentiating effect at the same α4β<sup>2</sup> nAChR. The fact that the synergism disappears in smoking olanzapine patients suggests that this interaction is highly non-linear in nature.

Smoking by itself tends to improve cognitive outcome, probably through its effect on α4β<sup>2</sup> nAChR GABA tone, perhaps underscoring its capacity for self-medication. Smoking can also enhance the effect of either memantine or AChE-I in a number of situations, but with both memantine and AChE-I added to antipsychotics the effect is hard to predict due to a number of non-linear interactions. This further underscores the importance of the excitation-inhibition balance in these cortical networks.

It is of interest to elaborate on the negative pharmacodynamic interactions of certain conditions. Increasing evidence suggest that information in the human brain is not encoded in simple firing rates; examples include oscillatory behavior of local field potential in the subthalamic nucleus of basal ganglia that code for motor symptoms (Little et al., 2013; Beudel et al., 2015). Our model readout also takes into account the interspike variability of the action potential train(Geerts et al., 2013) which might explain some of the non-monotonic dose-responses observed with quetiapine and aripiprazole, where increased firing rate does not lead to enhanced cognitive performance.

### Limitations of the Model

The network is calibrated on a working memory task, represents only the maintenance phase and therefore probably does not capture the intricacies of the different cognitive tasks. However, we would argue that the network also captures the strength of a memory trace representation that is a necessary step in a number of other cognitive tasks. For instance, in an episodic memory task, an existing memory trace needs to be retrieved from its memory bank and kept for a certain time in memory, although at a different time scale, to perform calculations and to compare it to a novel sensory stimulus.

A major issue is the choice of the biological processes and the changes associated with the pharmacological treatment. For instance, the allosteric effect of galantamine on nAChR has been documented to be as high as 70% (Samochocki et al., 2000) but a later study puts the effect more at 40% and lower (Samochocki et al., 2003). Furthermore, we have somewhat arbitrarily set the effect of smoking as a 20% increase in activation. The effect of such differences can be studied using sensitivity analyses.

Another major limitation is the "blind" pharmacology of antipsychotic drugs; unlike animal models, the computer model is bound by the available knowledge of the pharmacology of the drugs and does not take into account the possible effect of drugs on targets that are undocumented. However, the functional concentration of antipsychotics as measured by PET raclopride radiotracer in clinically relevant situations is in the low nM range, except for quetiapine. This effectively reduces the probability of off-target effects that might play a role at higher doses.

This simulation generates average values for "generic" patients on a fixed dose of an antipsychotic and, memantine and/or an AChE-I. In real-life practice, doses are often titrated to the best trade-off between efficacy and side-effects. In addition, the variability in clinical outcome might be affected by CYTP450 interaction between the different drugs where inhibition or stimulation of a CYTP450 enzyme by one drug might affect the levels of other drugs (so called pharmacokineticpharmacokinetic or PK/PK interactions). Finally, variability might be a consequence of different genotypes, for instance the COMTVal158Met and the 5-HTTLPR s/l that all can affect cognitive outcomes. In principle, the QSP platform is able to simulate virtual individual patients complete with PK/PK interactions, different individual doses of the different drugs and any combination of the common genotypes. As an example, we have implemented the COMTVal158Met genotype in the QSP platform using human imaging studies in unmedicated healthy volunteers as a different half-life of dopamine and norepinephrine in cortical synapses (Spiros and Geerts, 2012).

The major message from this simulation exercise is the unexpected impact of different antipsychotics each with their own fingerprint of pharmacological activity on the dose-response of the same drug-drug combination and/or combined with smoking. The non-linear interactions of drug pharmacodynamics plays an important role in real-life drug treatment, as polypharmacy is more the rule than the exception. This is particularly important when considering the design of clinical trials; careful consideration of inclusion/exclusion of comedications can substantially reduce the patient variability in the different treatment arms and increase the probability that a clinical signal can be detected. Failure to take these differential interactions into account might lead to reduction of the clinical signal as treatment arms become inadvertently populated with drugs that do not work or work adversely. We suspect that this might be one of the causes of clinical trial failures of drugs tested with augmentation therapy in Cognitive Impairment with Schizophrenia. Obviously other processes contribute to the variability in clinical trial outcomes such as different genotypes, different pathological baselines and PK variability. In this regard, we would argue that the concept of chlorpromazine equivalents (Beckmann and Laux, 1990) where the differentiation between antipsychotics is solely driven by their dose and their corresponding D2R occupancy, although very handy and appealing because of its simplicity, might need a substantial correction.

### References


An additional collorary is that new drugs with a specific selective pharmacology aimed at augmentation therapy could only be combined with specific antipsychotics, because of negative pharmacodynamic interaction of other antipsychotics with specific drug pharmacology of the new investigative compound. In addition, smoking through its effect on nicotinic receptors can significantly affect cognitive outcome. The results of this paper suggest that the complex poly-pharmacy profile of marketed antipsychotics can lead to non-linear pharmacodynamics interactions beyond their simple D2R occupancy and that this can significantly impact clinical outcome.

galantamine in schizophrenia: a sustained alpha7 nicotinic agonist strategy. Clin. Neuropharmacol. 31, 34–39. doi: 10.1097/wnf.0b013e31806462ba


with schizophrenia and schizoaffective disorder: a review and meta-analysis of the literature. CNS Drugs 24, 303–317. doi: 10.2165/11530260-000000000- 00000


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2015 Geerts, Roberts and Spiros. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

## Multitarget Multiscale Simulation for Pharmacological Treatment of Dystonia in Motor Cortex

Samuel A. Neymotin1, 2 \*, Salvador Dura-Bernal <sup>1</sup> , Peter Lakatos <sup>3</sup> , Terence D. Sanger 4, 5 and William W. Lytton1, 6, 7, 8

*<sup>1</sup> Department Physiology and Pharmacology, SUNY Downstate Medical Center, State University of New York, Brooklyn, NY, USA, <sup>2</sup> Department Neuroscience, Yale University School of Medicine, New Haven, CT, USA, <sup>3</sup> Nathan S. Kline Institute for Psychiatric Research, Orangeburg, NY, USA, <sup>4</sup> Department Biomedical Engineering, University of Southern California, Los Angeles, CA, USA, <sup>5</sup> Division Neurology, Child Neurology and Movement Disorders, Children's Hospital Los Angeles, Los Angeles, CA, USA, <sup>6</sup> Department Neurology, SUNY Downstate Medical Center, Brooklyn, NY, USA, <sup>7</sup> Department Neurology, Kings County Hospital Center, Brooklyn, NY, USA, <sup>8</sup> The Robert F. Furchgott Center for Neural and Behavioral Science, Brooklyn, NY, US*

#### Edited by:

*Hugo Geerts, In Silico Biosciences, USA*

#### Reviewed by:

*Gaute T. Einevoll, Norwegian University of Life Sciences, Norway Christian Meisel, University Hospital Carl Gustav Carus Dresden, Germany James B. Aimone, Sandia National Laboratories, USA*

> \*Correspondence: *Samuel A. Neymotin samn@neurosim.downstate.edu*

#### Specialty section:

*This article was submitted to Experimental Pharmacology and Drug Discovery, a section of the journal Frontiers in Pharmacology*

> Received: *17 February 2016* Accepted: *30 May 2016* Published: *14 June 2016*

#### Citation:

*Neymotin SA, Dura-Bernal S, Lakatos P, Sanger TD and Lytton WW (2016) Multitarget Multiscale Simulation for Pharmacological Treatment of Dystonia in Motor Cortex. Front. Pharmacol. 7:157. doi: 10.3389/fphar.2016.00157* A large number of physiomic pathologies can produce hyperexcitability in cortex. Depending on severity, cortical hyperexcitability may manifest clinically as a hyperkinetic movement disorder or as epilpesy. We focus here on dystonia, a movement disorder that produces involuntary muscle contractions and involves pathology in multiple brain areas including basal ganglia, thalamus, cerebellum, and sensory and motor cortices. Most research in dystonia has focused on basal ganglia, while much pharmacological treatment is provided directly at muscles to prevent contraction. Motor cortex is another potential target for therapy that exhibits pathological dynamics in dystonia, including heightened activity and altered beta oscillations. We developed a multiscale model of primary motor cortex, ranging from molecular, up to cellular, and network levels, containing 1715 compartmental model neurons with multiple ion channels and intracellular molecular dynamics. We wired the model based on electrophysiological data obtained from mouse motor cortex circuit mapping experiments. We used the model to reproduce patterns of heightened activity seen in dystonia by applying independent random variations in parameters to identify pathological parameter sets. These models demonstrated degeneracy, meaning that there were many ways of obtaining the pathological syndrome. There was no single parameter alteration which would consistently distinguish pathological from physiological dynamics. At higher dimensions in parameter space, we were able to use support vector machines to distinguish the two patterns in different regions of space and thereby trace multitarget routes from dystonic to physiological dynamics. These results suggest the use of *in silico* models for discovery of multitarget drug cocktails.

Keywords: dystonia, multiscale modeling, computer simulation, motor cortex, beta oscillations, corticospinal neurons, multitarget pharmacology, support vector machines

### 1. INTRODUCTION

A large number of physiomic pathologies can produce hyperexcitability in cortex. In motor cortex, this hyperexcitability will manifest as alterations in movement and muscle tone. At the most extreme, hyperexcitability leads to a seizure with uncontrolled movement, as seen in epilepsia partialis continuans. Lesser hyperexcitability produces a variety of hyperactive movement disorders, including tics, chorea, tremor, etc, whose pathophysiology is not restricted to cortex, but involves multiple brain areas including basal ganglia, thalamus, cerebellum, and others. We focus here on dystonia, a movement disorder that produces prolonged involuntary muscle contractions (Neychev et al., 2008; Crowell et al., 2012).

The large variety of dystonias of different etiologies may present with involvement of one or several parts of the body. Pediatric causes of dystonia include cerebral palsy and are generally distinct from adult-onset cases. Common adult dystonias are torticollis, causing involuntary head turning, and movement-overuse dystonias such as writers cramp. Despite these differences, dystonias in different patient populations are primarily treated with the same therapies. While most research in dystonia has focused on basal ganglia, much pharmacological treatment is provided directly at muscles. Similarly, we propose that treatment could be targeted elsewhere in the motor pathway, here focusing on motor cortex as a potential target for therapy.

As with many other movement disorders, the dystonias generally lack a reliable biomarker and are diagnosed by semiology, the assessment of signs and symptoms. However, all dystonias feature excessive muscle activation that is associated with hyperactivity in multiple motor areas associated with movement activation. Electrophysiological studies of dystonia patients confirms a pattern of hyperactivation in cortex. Healthy individuals show low beta oscillations (∼15–20 Hz) in motor cortical local field potential (LFP). This beta is reduced in amplitude and synchrony during movement (Jasper and Penfield, 1949; Pfurtscheller and Aranibar, 1979; Crone et al., 1998; Miller et al., 2007). In dystonia patients, motor cortex shows increases in neuronal activity levels (Nobrega et al., 1995; Pratt et al., 1995), with relatively high beta amplitude and high functional connectivity at the beta frequency (Schnitzler and Gross, 2005; Jin et al., 2011). There is also excessive neural synchrony both at rest and in certain phases of movement (Toro et al., 1994; Kristeva et al., 2005; Mallet et al., 2008; Crowell et al., 2012).

Some dystonias, in common with several other movement disorders, are thought to have their origin in the basal ganglia. Other dystonias, such as those associated with cerebral palsy and with movement overuse, probably have a strong cortical component. In all cases, however, the interconnections of brain motor systems makes it clear that multiple brain areas will be "in the loop" of abnormal activity. Following some primary insult or insults to a brain area, a secondarily-involved brain area will contribute further to the disorder by reacting to the alterations in input activity through its own homeostatic responses. In some cases these homeostatic changes may be compensatory so as to reduce the severity of the symptoms. However, in other cases, plasticity may actually exacerbate the abnormal movements (Sanger et al., 2003; Neychev et al., 2008; Casellato et al., 2014).

There are at least two, and perhaps more, cerebello-thalamostriato-cortical loops that play a role in movement disorders. There may also be additional contributions from still longer loops involving recurrent connections from spinal cord or from muscle spindles. One or more of these sites may have associated pathology. Regardless of the locus of primary pathology, multiple sites are potential targets where therapy could interrupt pathophysiological dynamics. Currently, brain pharmacotherapy often fails and patients are treated with botulinum toxin to partially paralyze muscles by blocking nicotinic cholinergic transmission at the affected muscle. Another treatment is deep brain stimulation using implanted electrodes. In this paper, we take two or three steps back from the level of muscle treatment by proposing interventions at the level of motor cortex.

Complex multifocal diseases may require complex multitarget treatments (Viayna et al., 2013). In the context of brain disease, multitarget therapy can hit multiple brain regions or multiple receptors in a region or both. High-level models that include many brain areas can assist in understanding how different brain areas contribute to a disorder (Sanger and Merzenich, 2000; Sanger, 2003; Hendrix and Vitek, 2012; Kerr et al., 2013). However, these models typically lack biological detail, making them unsuitable for assessing the impact of specific pharmacological manipulations. Detailed models are not yet elaborated to the point of handling multiple brain areas but do provide the details needed to assess pharmacological intervention more directly.

Single agent treatments for disease are traditionally tested in vitro or in vivo. As noted above, single agent treatments for dystonia have not had much success. There is, however, the potential for success with multitarget drug cocktails that could target multiple locations in the brain, or multiple drug receptor targets at a single location, or both (Delnooz and van de Warrenburg, 2012). Due to combinatorial explosion, evaluating combinations of drugs in different dosages in this way can not be readily done in tissue and is most feasible in silico (Viceconti et al., 2008; Kohl and Noble, 2009; Lytton et al., 2014; Action, 2016; Viceconti et al., 2016). In this study, we use our detailed multiscale model of primary motor cortex to assess potential multitarget pharmacological therapies for treatment of dystonia. The model contains 6 cortical layers with multiple classes of excitatory and inhibitory neurons, using wiring based on mouse microconnectomic data (Shipp, 2005; Weiler et al., 2008; Kiritani et al., 2012; Hooks et al., 2013). Excitatory neurons contain intracellular molecular mechanisms that contribute to persistent activity and hyperexcitability (Neymotin et al., 2016). These mechanisms include endoplasmic reticulum associated calcium stores released by activation of IP3Rs, and ryanodine receptors, both with affinity for caffeine, an agent that can exacerbate dystonia symptoms (Richter and Hamann, 2001). Plasma membrane calcium, sodium, and potassium channels also contribute to cellular excitability.

Since our model does not include spinal cord and muscle, we defined dystonia pathology as a state of cortical hyperactivation characterized by increased beta oscillations with excessive and hypersynchronous firing in layer 5 corticospinal neurons. These layer 5 neurons project downward to brainstem and spinal cord, and their sustained firing would lead to the increased muscle contractions of dystonia. We distinguished the hyperexcitability of dystonia from the still greater hyperexcitability of a seizure by excluding simulations that showed higher levels of activity with higher frequency oscillation and a strong tendency to "latch-up" through multicell depolarization blockade (Lytton and Omurtag, 2007). Classification in 11-dimensional space demonstrated that we could identify different regions in parameter space for these different states—baseline normal, dystonia, epileptiform—and predict pharmacological combinations that would lead from pathology back to the physiological activity state. As in our previous investigations of epilepsy (Lytton and Omurtag, 2007), we found multiple parameter combinations that were consistent with the pathological state, as well as multiple parameter combinations to produce our baseline physiological state. Such parameter degeneracy is typical of complex neural systems (Edelman and Gally, 2001; Golowasch et al., 2002).

### 2. MATERIALS AND METHODS

Network simulations consisted of 1715 reduced compartmental cell models with single compartments for inhibitory cells and five compartments for pyramidal cells, arrayed by layer with connectivity taken from experimental results on motor cortex (Weiler et al., 2008; **Figures 1A,B**). Parallel-conductance electrophysiological simulation in the pyramidal cells was complemented by chemophysiological simulation focused on Ca2<sup>+</sup> handling, based on our prior models (Neymotin et al., 2015, 2016; **Figure 1C**).

Simulations were run in the NEURON (version 7.4) simulation environment (Carnevale and Hines, 2006) utilizing the reaction-diffusion (RxD) Python module (McDougal et al., 2013a,b) and NMODL (Hines and Carnevale, 2000). Two seconds of simulation time took ∼3 min using 24 nodes on a Linux cluster with parallel NEURON, run with a fixed time-step of 0.1 ms. The full model is available on ModelDB (https://senselab.med.yale.edu/ModelDB/ShowModel. cshtml?model=189154).

We briefly describe the scales of the multiscale model from smaller to larger in the following sections (**Table 1**). For more details, readers are referred to our previous papers (Neymotin et al., 2015, 2016).

### 2.1. Intracellular Molecular Scale

Our Ca2<sup>+</sup> dynamics (**Figure 1C**), are based on (Neymotin et al., 2016). We modeled a one-dimensional RxD system of intracellular neuronal Ca2<sup>+</sup> signaling in all compartments of neocortical pyramidal (PYR) neurons. Within each compartment, we modeled cytosolic and endoplasmic reticulum (ER) sub-compartments by using a fractional volume for each.

IP<sup>3</sup> was produced through a reaction sequence initiated by glutamate binding to the metabotropic glutamate receptor (mGluR), based on a reaction scheme developed by Ashhad and Narayanan (2013) (ModelDB #150551). IP<sup>3</sup> diffused outward from the synapse location and decayed following first-order kinetics (τIP<sup>3</sup> of 1 s). Baseline mGluR synaptic weight was normalized to represent the increase in the amount of glutamate bound to mGluR. Extracellular glutamate did not diffuse but was represented by a local Glu value that was incremented in response to an event delivered due to a presynaptic spike. Glu showed bind/unbind kinetics on mGluR and was eliminated by first-order degradation (lower left of **Figure 1C**).

The ER Ca2<sup>+</sup> model involves IP<sup>3</sup> receptors (IP3Rs), ryanodine receptors (RYR) (Sneyd et al., 2003), SERCA pumps, and a Ca2<sup>+</sup> leak. IP3R dynamics involved a slow Ca2<sup>+</sup> inactivation binding site state (De Young and Keizer, 1992; Li and Rinzel, 1994). The SERCA pump is a pump rather than a channel and so is modeled with Hill-type dynamics. Calcium buffering followed Ca + B 5 −−−− ↽−−−−⇀ 9.5·10−<sup>4</sup> CaB where B is diffusible buffer with diffusion coefficients D = 0.043 µm<sup>2</sup> /ms for both B and CaB, about half the rate of Ca2+diffusion (Anwar et al., 2012). Calcium extrusion across the plasma membrane was modeled by first-order decay with τex = 5 ms.

### 2.2. Synapses

AMPA/NMDA synapses were modeled by standard NEURON double-exponential mechanisms (**Table 2**). All excitatory projections were mixed AMPA (rise,decay τ : 0.05, 5.3 ms) and NMDA (rise,decay τ : 15, 150 ms). NMDARs were scaled by 1/(1 + 0.28 · Mg · exp(−0.062 · V)); Mg = 1mM (Jahr and Stevens, 1990). 13% of INMDA was carried by Ca2+(Spruston et al., 1995). AMPA and NMDA receptors had reversal potentials of 0 mV.

Inhibitory synapses were mediated by GABA<sup>A</sup> and GABA<sup>B</sup> receptors. GABA<sup>A</sup> synapses were modeled with a doubleexponential mechanism. The GABA<sup>B</sup> synapse had second messenger connectivity to a G protein-coupled inwardlyrectifying potassium channel (GIRK). LTS cells connected to apical dendrites of PYR cells using GABA<sup>A</sup> receptors (GABAAR; rise,decay τ : 0.2, 20 ms) and GABA<sup>B</sup> receptors (GABABR) and onto somata of FS and other LTS cells with GABA<sup>A</sup> Rs (rise,decay τ : 20, 40 ms). GABAARs had reversal potentials of −80 mV, and GABABRs −95 mV. GABABRs provide longer-lasting activation compared to GABAARs.

### 2.3. Cell Scale

The network consisted of pyramidal cells (PYR; 3 apical dendrite compartments, 1 basal dendrite compartment, 1 somatic compartment), fast spiking soma-targeting interneurons (FS; one compartment) , and dendrite-targeting low-threshold spiking interneurons (LTS; one compartment; Wang and Buzsaki, 1996; Wang, 2002; Monyer and Markram, 2004; Bartos et al., 2007; Neymotin et al., 2011a,b; **Tables 3**, **4**). Reaction-diffusion mechanisms (Ca2+,IP3,buffer) were restricted to the PYR cells in this network. Properties of pyramidal neurons in the different layers were identical except for apical dendrite length which is longer in deep pyramidal neurons than in superficial (Hay et al., 2011; Castro-Alamancos, 2013): 900 µm in L5-6; 450 µm in L2/3 and L4.

Voltage-gated ionic current models were based on prior models of our own and others (McCormick and Huguenard,

lines proportional to synaptic weights. E cells synapse with AMPAR/NMDARs; I cells synapse with GABAAR / GABABRs. Circles with self-connects denotes recurrent connectivity. (A) Excitatory connections. E5P projects to spinal cord (not modeled). (B) Inhibitory connections. (C) Chemical signaling in pyramidal cells showing fluxes (black arrows) and second- (and third- etc) messenger modulation (red back-beginning arrows). We distinguish membrane-associated ionotropic and metabotropic receptors and ion channels involved in reaction schemes in red (in reality, it is likely that almost every membrane-bound protein is modulated). External events are represented by yellow lightning bolts—there is no extracellular diffusion; the only extracellular reaction is glutamate binding, unbinding, and degradation on mGluR1 after an event. Ca2<sup>+</sup> is shown redundantly in blue—note that there is only one Ca2<sup>+</sup> pool for extracellular, 1 pool for cytoplasmic, and 1 pool for ER (PLC, phospholipase C; DAG, diacyl-glycerol; cAMP, cyclic adenosine monophosphate; PIP2 , phosphatidylinositol 4,5-bisphosphate). Adapted from Figure 1 of Neymotin et al. (2016).

1992; Migliore et al., 2004; Stacey et al., 2009; Neymotin et al., 2011b,a, 2013). Voltage sensitive channels generally followed the Hodgkin-Huxley parameterization, whereby x˙ = (x<sup>∞</sup> − x)/τ<sup>x</sup> (x = m for activation particle and h for inactivation particle). Steady-state x<sup>∞</sup> and time constant τ<sup>x</sup> are either related to channel opening α(V) and closing kinetics β(V) as x<sup>∞</sup> = α/(α + β), τ<sup>x</sup> = 1/(α + β), or are directly parameterized: x∞(V), τx(V). Kinetics for channels were scaled by Q<sup>10</sup> from an experimental temperature (where available) to simulation temperature of 37◦C. Q<sup>10</sup> = 3 was used when no experimental value was available. All cells contained leak current, transient sodium current INa, and delayed rectifier current IK−DR, to allow for action potential generation. Additionally, PYR cells contained in all compartments IK−A, IK−<sup>M</sup> providing firing-rate adaptation (McCormick et al., 1993). Pyramidal cells also had Ih, voltagegated calcium channels (VGCCs) in all compartments: IL, IT, I<sup>N</sup> (Kay and Wong, 1987; McCormick and Huguenard, 1992; Safiulina et al., 2010; Neymotin et al., 2015), and SK and BK

#### Neymotin et al. Multiscale Modeling for Dystonia Therapies

#### TABLE 1 | Summary of model.


*E (I) denote excitatory (inhibitory) neurons. No plasticity modeled (Table format based on Nordlie et al., 2009).*



#### TABLE 3 | Summary of neuron models.


*E (I) denote excitatory (inhibitory) neurons. Reaction-diffusion (RxD) mechanisms/ compartments described more fully in intracellular scale.*

calcium-activated potassium currents (IKCa). LTS cells contained IL, non-Ca2+-dependent Ih, SK, and Ca2<sup>+</sup> decay.

HCN channels in different cell types have somewhat different voltage dependence and different kinetics (Hagiwara and Irisawa, 1989; Schwindt et al., 1992; Chen et al., 2001; Wang et al., 2002; Robinson and Siegelbaum, 2003). The hyperpolarizationactivated HCN current I<sup>h</sup> used in pyramidal cells was modeled with second messenger and calcium dependence taken from Winograd et al. (2008) (ModelDB #113997), and modified to provide the faster voltage-sensitivity time constants found in cortex (Harnett et al., 2015), and provides PYR cells longerlasting firing activity via augmentation of the HCN channel conductance. The mechanism for Ca2<sup>+</sup> regulation of HCN channels in PYR cells in Winograd et al. (2008) is modeled empirically in order to produce the relationship between cytosolic Ca2<sup>+</sup> levels and I<sup>h</sup> activation without simulating the details of Ca2<sup>+</sup> effects on adenyl cyclase (see schematic for HCN chan in **Figure 1C**).

g¯<sup>h</sup> was 0.0025 S/cm<sup>2</sup> in PYR soma, basal dendrites and exponentially-increasing in apical dendrites with distance from soma with 325 µm space constant, hence e-fold augmented at 325 microns as described by Kole et al. (2006). Apical dendrite IK−DR, IK−A, IK−<sup>M</sup> density also increased with the same length constant, based on data showing HCN and Kv channel colocalization (Harnett et al., 2015, 2013).

### 2.4. Network Scale

The network consisted of 1715 cells (**Table 4**). The network contained 157,507 synapses for an overall connection density of ∼5% (see **Table 6**). PYR cells synapsed onto each-other's dendrites. PYR-to-PYR synaptic locations on the dendrite were randomized between basal and apical compartments (Markram et al., 1997). PYR cells synapsed onto somata of FS and LTS cells (single-compartment models).

Neuronal populations were arranged by cortical layer based on our prior models (Neymotin et al., 2011a,c; Chadderdon et al., 2014; Neymotin et al., 2016), with additional data from direct measurements from mouse motor cortex (Shipp, 2005; Weiler et al., 2008; Kiritani et al., 2012; Hooks et al., 2013), and recent experiments which demonstrate a thin layer 4 in mouse motor cortex (Yamawaki et al., 2014). The network consisted of 13 populations of 3 excitatory cell types, intratelencephalic (IT), pyramidal-tract (PT), and corticothalamic (CT), and 2 inhibitory cell types, low-threshold spiking (LTS) and fastspiking (FS). These were distributed across cortical layers 2/3, 4, 5a, 5b, and 6 (Harris and Shepherd, 2015), with cell numbers for each population based on estimated cell densities and volume (**Table 4**). The volume of each population was calculated assuming a horizontal area (x and z dimensions) of 120 × 120 µm, and a layer-dependent cortical depth range (y dimension).

Connection probabilities pij (**Tables 5**, **6**) of presynaptic excitatory populations were dependent on pre- and pothstsynaptic type and layer. For presynaptic inhibitory populations, connection probabilities inversely scaled based on distance pij = p¯ij · exp(− p (dx<sup>2</sup> + dz<sup>2</sup> )/15), in x, z plane perpendicular to the y-direction of layering. Connection probabilities and weights are based on data from rodent motor cortex mapping (Weiler et al., 2008; Lefort et al., 2009; Anderson et al., 2010; Fino and Yuste, 2011; Apicella et al., 2012; Kiritani et al., 2012). Individual neurons were placed randomly with uniform distribution. Weights from E cells displayed in **Table 6** are for the AMPA synapses, with colocalized NMDA weights at 10% of AMPA weights. Synaptic delays were randomized between 1.8 and 5 ms with additional delay based on distance.


TABLE 4 | Network Population, including normalized and nominal cortical depth range (ynormRange, yRange, neuron density, and number of cells).

*PYR, pyramidal; IT, intratelencephalic; PT, pyramidal tract; CT, corticothalamic; FS, fast spiking, LTS, low-threshold spiking.*

#### TABLE 5 | Summary of network connectivity rules.


*pij denotes probability of connection between type i and j; wij denotes weight. Parameters by pre- and post-synaptic type in* Table 6*.*

Background activity was simulated by excitatory and inhibitory synaptic inputs following a Poisson process, sent to all cells, representing ongoing drive from other cortical areas and other inputs. These inputs were selected to maintain low-frequency firing of neurons within the model, which would not fire otherwise, due to small network size and the requirement for multiple synaptic inputs to trigger a postsynaptic spike (Neymotin et al., 2011a). The strength of these background inputs was not based on the full source of inputs from all upstream brain areas, since these inputs are not completely mapped.

### 2.5. Simulation Variations

We ran over 5800 simulations, randomly varying each of the following parameters using an independent normal distribution: 1. E neuron mGluR density (mGluR); 2. E neuron ER RYR density (RYR); 3. E and I neuron HCN channel density; 4. E and I neuron fast Na<sup>+</sup> channel density (Na<sup>f</sup> ); 5. E neuron Kdr channel density; 6. E neuron K<sup>a</sup> channel density; 7. E neuron K<sup>D</sup> channel density; 8. E neuron K<sup>M</sup> channel density; 9. E neuron SK calcium-activated potassium channel density; 10. E neuron BK calcium-activated potassium channel density; 11. E and LTS neuron voltage-gated calcium channel (VGCC) density.

Means and standard deviations differed for the different parameters and were selected to allow each simulation to maintain activity. A subset of the simulations was used for the analyses described (**Table 7**).

We ran simulations with initial calcium concentration in the ER set to 1.25 mM (Bygrave and Benedetti, 1996), to mimic the effects of ER calcium priming via prior excitatory synaptic stimulation (Ross et al., 2005; Hong and Ross, 2007; Fitzpatrick et al., 2009; Neymotin et al., 2016).

We categorized the simulations into distinct groups by noting major differences in activity across parameter sets (**Table 8**). From the full set of 5867 simulations, 1505 did not display any firing due to random variations in ion channel densities which led to low neuronal activity (**Table 7**). The remaining 4341 simulations were Active due to higher neuronal activity, e.g., partially caused by the higher average Na<sup>f</sup> density in these simulations. Of these 4341 Active simulations, 1077 exhibited epileptic latch-up dynamics—periods of intense activity which led to depolarization blockade (Na<sup>+</sup> channel inactivation; Lytton and Omurtag, 2007). These periods where neurons did not fire lasted 200–300 ms (gaps between E5P spikes: E5P gap in **Table 8**). We categorized the top and bottom 2nd percentile of the Active/non-latch-up simulations ranked by E5P firing rate into dystonia (n = 65) and physiological (n = 65) sub-sets. We used E5P firing rate as a criterion for dystonia classification because E5P neurons project downward to brainstem spinal cord, and sustained overactive E5P firing can lead to the tonic muscle contractions symptomatic of dystonia.

### 2.6. Data Analysis

We formed multiunit activity (MUA) time-series, which count the number of spikes in each bin (10 ms) for a given population. To calculate neuronal population rhythms, we took the power spectral density (PSD) of the mean-subtracted MUA time-series; we then calculated the peak frequencies and amplitudes in the PSD. We used the average Kendall's τ nonparametric rank-correlation coefficient (Kendall, 1938; Knight,

#### TABLE 6 | Network Connectivity Parameters.


*p*¯*ij and wij are distance-independent probability of connections from Pre to Post neuronal types and synaptic weights, respectively.*

1966) between pairs of neuron binned spike train time-series for calculating the synchronization of populations of neurons (denoted population-synchrony). Kendall's τ non-parametric rank correlation, defined as:

$$
\sigma = \frac{n\_c - n\_d}{\frac{1}{2}n(n-1)},
$$

is used with these data. Kendall's τ is a normalized difference between concordant (nc) and discordant pairs (nd); ties are taken into account by the normalizing term, <sup>1</sup> 2 n(n − 1) , which represents the total number of ordered pairs in the

time-series. We used the Python scikit-learn library (Pedregosa et al., 2011) for performing principal component analysis (PCA) and support-vector machine (SVM) classification (Cortes and Vapnik, 1995; Orrù et al., 2012). Dystonia and physiological simulation classes were characterized on the basis of layer 5 corticospinal pyramidal neuron (E5P) firing rates. The clearest examples of both classes (bottom/top 2nd percentiles as physiological/dystonia classes) were used for the majority of the analyses described in the Results (**Figures 3**–**8**). The NuSVC variant of SVMs was used to classify physiological and dystonia simulations and to find which simulation parameters contributed the most to classification accuracy. SVM inputs


TABLE 7 | Parameter ranges (average ± standard deviation) for all simulations (n = 5867), active simulations (n = 4341), latch-up simulations (n = 1077), active/non-Latch-up simulations (n = 3264), physiological simulations (n = 65), and dystonia simulations (n = 65).

*Plasma membrane ion channel conductance density values are in S/cm*<sup>2</sup> *. mGluR and RYR density are in arbitrary units used to scale channel conductance.*

TABLE 8 | Dynamic measures (average ± standard deviation) for All simulations (n = 5867), Active simulations (n = 4341), Latch-up simulations (n = 1077), Active/Non-Latch-up (n = 3264), physiological simulations (n = 65), and dystonia simulations (n = 65).


*E5P gap measures number of 300 ms gaps between individual E5P neuron firing times; E5P MUA amplitude and E5P MUA beta amplitude in arbitrary units; E5P FV sim measures similarity between E5P population firing rate vectors using average pairwise Pearson correlation.*

were vectors consisting of normalized parameter values. Each of these input vectors was labeled into either of two distinct binary classes: physiological (0) or dystonia (1). SVM parameters, including kernel type (linear, polynomial, radial-basis function), γ , tolerance, ν, and polynomial degree were selected using a grid search with N-fold cross validation run 10 times for each combination of parameters. SVM classification accuracy surpassed the accuracy of other machine learning methods, including logistic regression (not shown). Figures were drawn with Matplotlib (Hunter, 2007).

### Neymotin et al. Multiscale Modeling for Dystonia Therapies

## 3. RESULTS

### 3.1. Simulation Overview

We ran over 5800 network simulations, randomizing 11 ion channel/receptor densities independently. A typical 2 s simulation took ∼3 min using 24 cores on Linux with parallel NEURON. After running simulations, we calculated neuronal population firing rates, synchronization, and power spectra.

### 3.2. Characterization of Dystonia Pathophysiology

Simulations were grouped into physiological and pathological based on differences in firing patterns (**Table 8**, **Figure 2**). 1505 of 5867 simulations produced no activity. The remaining simulations were characterized as physiological or pathological. Pathological simulations showed increased activity. High activity alternating with latch-up condition was defined as an epileptiform simulation with periods of >200 ms of depolarization blockade across multiple cells (1077 simulations). 1077 simulations were classified as epileptiform based on activity latch-up resulting in sustained periods. The different classes of simulations formed distinct clusters in multiple slices of excitatory corticospinal (ESP) activity feature-space (**Figure 2**). Physiological simulations showed E5P rates ≤2 Hz with low to intermediate E5P firing vector (FV) similarity. Dystonia simulations primarily occupied the upper-right quadrant of the scatterplot in **Figure 2A**, but displayed either high or low FV similarity which overlapped with the range of values displayed by the physiological simulations. Epileptiform simulations had intermediate average E5P rates due to high activity alternating with periods of quiescence caused by depolarization blockade. Across simulation types, higher E5P firing increased the excitatory drive to I5 neurons, causing increased I5 neuron firing (**Table 8**). Higher I5 and E5P neuron firing then caused higher E5P synchronization via recurrent E5P excitation and feedback inhibition (**Figure 2B**). Stronger E5P and I5 interactions then increased beta rhythm amplitude (**Figure 2B**), however with substantial variability. Peak oscillatory frequency was held relatively stable across simulations (**Table 8**). Physiological and epileptiform simulations had lower overall E5P synchrony and beta power compared to the dystonia simulations, which occupied the upper-right quadrant of **Figure 2B**.

E5P FV similarity showed temporal recurrences which further distinguished the three simulation types (**Figure 2C**). The physiological simulation showed intermediate self-similarity (0.17) due to sparse firing of different subsets of pyramidal cells at different times. In contrast, the dystonia simulation firing patterns showed strong self-similarity (0.56) and recurrence over time (recurring orange/red blobs in **Figure 2C**), indicating stereotyped dynamics. The example epileptiform simulation showed relatively weak self-similarity (0.16) due to its two distinct firing patterns: high E5P synchrony alternating with E5P silence produced by depolarization blockade. Epileptiform and dystonia simulations showed a brief period of high similarity when the epileptiform simulation showed strong oscillations during the initial period. There was weak similarity between epileptiform and physiological (0.12) and dystonia and physiological (0.22) simulations, indicating that both pathological dynamics were distinct from the physiological.

E5P neurons in a representative physiological model fired sparsely with low synchrony (population-synchrony = 0.01; **Figures 3A,D**; Supplementary Figure 1 has all physiological rasters), with multiple downstream effects. Low excitatory drive from E5P to I5 and I5L neurons caused them to fire slowly. This low L5 inhibition allowed E5a neurons to fire quickly. The weak E5P and L5 interneuron interactions produced only weak beta rhythms which were confined to layer 5 (**Figure 4A**).

In a representative dystonia simulation, E5P neurons had sustained, synchronous, rapid firing (**Figures 3B,D**; Supplementary Figure 2 shows all dystonia simulation rasters). This promoted strong, continuous layer 5 interneuron activation. The L5 interneurons then suppressed E5a intratelencephalic

FIGURE 2 | Distinct dynamics in in physiological, dystonia, and latch-up simulations. (A) Average E5P firing rate vector (FV) similarity vs. average E5P firing rate for individual simulations. (B) E5P MUA Beta oscillation amplitude vs. E5P synchrony for individual simulations. (A,B) [light blue: physiological, purple: dystonia, orange: epileptiform, black: remaining Active simulations, large circles represent simulations shown in (C) and Figure 3]. (C) Pearson correlations between all pairs of E5P FVs. Solid black lines demarcate FVs from example physiological, dystonia, and epileptiform simulations. All FVs used 50 ms intervals, forming 40 FVs per 2 s of simulation.

neurons, which fired at reduced rates. In contrast, E5b firing increased with the faster E5P firing, due to excitation spreading in the network. The relatively high recurrent connectivity (18% density) and strong synaptic weights between E5P neurons allowed the E5P neurons to remain activated despite strong feedback inhibition. The strong feedback inhibition also activated the E5P HCN channels, which produced rebound excitation. The strong E5P activation coupled with the feedback inhibition also enabled E5P neurons to synchronize (population-synchrony = 0.83; vertical stripes in **Figure 3B**) at a strong beta rhythm (∼20 Hz; **Figure 4B**). These synchronous beta rhythms also spread to other populations and layers (E2, I5, I5L, E5b, and E6).

Epileptiform simulation also displayed strong intermittent beta oscillations and strong synchrony (population-synchrony = 0.05; **Figures 3C**, **4C**), but this activity alternated with lengthy periods (200–300 ms) where E neurons were not firing due to depolarization blockade. Even with these periods of

depolarization blockade, most E neurons fired at higher average rate than in the physiological simulations (**Figure 3D**). Such increased synchrony with high excitatory cell activity is seen in epilepsy patients (Meisel et al., 2015). In contrast to the dystonia simulations, the synchronous periods of epileptiform oscillations were largely confined to layer 5 and did not spread to other layers.

### 3.3. Need for Multitarget Approach

No individual parameter determined physiological vs. dystoniadynamical-condition in the network (**Figure 5**). Therefore, no single parameter adjustment would routinely provide an effective "treatment" that would reliably restore physiological activity in most pathological models. We therefore went on to explore whether multitarget manipulation would be able to find such treatment routes.

Although no single parameter could predict physiological vs. pathological dynamics, the outliers of certain individual parameters were predictive. At the pathological margin, simulations had parameters which are expected to produce high activity: high Na<sup>+</sup> or Ca2<sup>+</sup> channels promoting inward currents, high HCN channel densities providing high resting membrane potential (RMP), and low K<sup>+</sup> channel densities again producing depolarization and reduced repolarization with spiking.

Further evidence for lack of predictability of dynamics based on parameters, comes from viewing the parameters in all 11 dimensions organized into 2 classes by dynamics. The parameter space showed substantial heterogeneity in the patterns producing pathology (**Figure 6A**), with weak intra-class clustering (**Figure 6B**). Correlation between parameter vectors of each simulation averaged 0.06 for physiological simulations, 0.07 for pathological simulations, with weak -0.05 anticorrelation between pathological and physiological simulations. The low correlations in both the physiological simulations (lower-left corner of **Figure 6B**) and the pathological simulations (upperright corner of **Figure 6B**) demonstrate that there is widespread degeneracy in the parameter sets that produce either the physiological or pathological states. Some of this degeneracy is unsurprising: for example K<sup>+</sup> channels with similar time courses of activation can substitute for one another to some extent. Other degeneracy is more complex and involves higher-order dynamic compensation.

### 3.4. High Dimensional Separation of Physiological and Pathological Parameters

Because of the difficulty of separating pathological from physiological with these high dimensional parameter sets, we used a SVM classification to create a separation (termed a maximum margin hyperplane) separating parameter sets producing physiological dynamics from parameter sets producing pathological dynamics. We started by training SVMs using only two parameters in combination (**Figure 7**). In order to test the efficiency of this separation, we separated out our target sets (physiological vs. pathological) into two subsets of each to serve as training and testing sets to evaluate the adequacy of the separation. By trying various random training and testing sets we got a mean and standard error for each case. Many two-parameter predictions were below chance (0.5) indicating that the SVM could not separate physiological from pathological based on that parameter pair. Two-parameter SVMs could accurately classify when the parameter pair included Na<sup>f</sup> density—the strongest predictor of excitability. The best prediction came with high Na<sup>f</sup> and low Kdr. Logistic regression methods were also tried to perform this two-dimensional separation but did not perform as well as SVM.

Going beyond 2 parameters, SVM classification accuracy increased regularly with the number of parameters used

(**Figure 8**), suggesting that a multi-target drug approach beyond two targets might produce greater effect. Moving to higher and higher dimensional spaces, we checked all possible parameter combinations at each dimensionality. In **Figure 8**, we report the parameter combination that was most predictive—e.g., at 6 dimensions we report just one of the 462 combinations of six from 11 parameters. Looking at the red blocks below, we can identify that the six dimensions that provide best prediction are Na<sup>f</sup> , four of the K<sup>+</sup> channels, and VGCC. Predictability increases up through six parameters, then plateaus, and then falls off due to the extreme sparseness of data. This sparseness was due to the so-called curse of dimensionality: given a constant number of data points n, the density falls off #bin-fold with each increase in dimension, where #bin is the binning of the space in one

dimension. Because of this, any high-dimensional method will tend to underestimate predictive strength given a limited amount of data (Bishop, 2006; Noble, 2006).

This multi-target SVM approach revealed the parameters that had the highest contribution to producing or preventing dystonia. Na<sup>f</sup> density was the most predictive parameter across all numbers of parameters used (horizontal red stripe at top of **Figure 8B**), as had been also shown using 2 dimensions alone (**Figure 7**). Again confirming the 2-dimensional result, the next most predictive parameters was Kdr. Following these came Ka, Kd, BK, SK, and VGCC densities which also significantly contributed to accurate predictions, due to their strong influence on E neuron excitability. mGLUR, RYR, and K<sup>m</sup> densities showed lesser contributions.

Increasing the percentile cutoffs for categorizing physiological from pathological simulations from the 2nd percentile to 7th percentile decreased prediction accuracy but still demonstrated the value of multitarget changes (**Figure 9**). The left column shows the same result as **Figure 8**: accuracy increased (colormap) as one goes from fewer to more parameters (bottom to top). By including more exemplars on both the physiological and pathological sides, we moved away from the best exemplars

and obtained less distinction between the two parameter sets. However, at all percentiles, there was an initial increase in classification accuracy with continued increase up to or beyond 3 parameters. This increase then declined as the number of parameters increased further due to the aforementioned sparseness at high dimensionality.

### 4. DISCUSSION

We developed a multiscale model of primary motor cortex to explore multitarget pharmacological therapies for dystonia. We searched parameter space—channel and receptor densities to create a set of models to contrast dystonia dynamics with physiological dynamics. Dystonia simulations displayed high excitability and synchrony in layer 5 corticospinal neurons (E5P), and strong beta oscillations which spread between cortical layers (**Figures 3B**, **4B**). Dystonia simulations could be distinguished from epileptiform simulations primarily by the presence of periods of latch-up with depolarization blockade in the epileptiform simulations. Physiological simulations had low excitability, asynchronous firing, and weak beta oscillations (**Figures 3C**, **4C**). Attempts to use high-dimensional visualization techniques to find potential therapeutic directions in the parameter space were limited by the solution degeneracy in the 11-dimensional parameter space with scattered parameter vectors with low correlation (**Figure 6**). We therefore turned to a SVM classification to identify hyperplanes in high-dimensional space that would separate the two populations. As expected, the major spike generating channels, Na<sup>f</sup> and Kdr were the primary determinants of excitability, followed by additional potassium channels and calcium channels. We did not assess pharmacological effects on synapses, which would be useful to do in the future.

### 4.1. Biological Degeneracy and Personalized Therapy

Degeneracy of mechanism is a major theme in biology (Edelman and Gally, 2001), meaning that there are many different ways that a biological system can produce a particular shape in the case of an immunoglobulin, or a particular dynamics in the case of a neural system. Such degeneracy has been shown directly in the stomatogastric ganglion of lobster, where a particular cell type produces its stereotyped dynamics using many different combinations of ion channel densities (Golowasch et al., 2002). Associated with this degeneracy is failure of averaging averaging across parameter sets that produce the dynamics gives a set of parameter values that do not produce the same dynamics.

In the context of brain physiology, this means that we can expect that individuals differ in the details of how their

motor cortex produces oscillations and contributes to movement. Similarly, we can expect that individuals differ in the details of their pathology. From a pharmacological perspective this argues that we may see greater benefit from personalized medicine identifying the high-dimensional complex of pathological parameters in a particular patient in order to treat them with their own individualized cocktail of multitarget drugs to produce a dynamics that falls somewhere in the physiological regime. To this might also be added complementary individualized, perhaps multi-locus, brain stimulation (Kerr et al., 2012; Song et al., 2013; Chadderdon et al., 2014; Hiscott, 2014; Nelson and Tepe, 2014; Dura-Bernal et al., 2016). Such a personalized approach would require much more intensive, and more costly, diagnostic procedures of a type that is currently only used by epilepsy surgery centers, which typically require invasive methods to perform their diagnostic tests.

Due to the degeneracy, parameter averaging failed in our dataset—using the average of all parameters sets that produce pathological simulations does not give a pathological simulation. However, the ability of the SVM method to separate pathological from physiological populations in high dimensional parameter space does suggest that there may be some benefit to pushing all patients in that direction through a multitarget pharamacological cocktail. In future studies, we plan to test this explicitly in the simulations in order to determine what percentage improve, what percentage worsen and what percentage remain essentially unchanged with such an average treatment. This assessment will require a larger number of simulated patients than we have thus far accumulated.

### 4.2. Multilocus, Multitarget, Multiscale Approaches for Treating Dystonia

In general, single target pharmacology has not been effective in dystonia (Fahn, 1987). As in other complex diseases, many of the treatments for dystonia have highly variable effectiveness and must be used at high doses that produce side-effects (Jankovic, 2006). For these reasons, botulinum toxin, targeting the final endpoint —the muscle movement—is commonly used as a treatment (Jankovic, 2006; Sanger et al., 2007; Bragg and Sharma, 2014). Deep-brain stimulation, an invasive procedure, is also used to partially restore normal brain dynamics (Tarsy, 2007; Johnson et al., 2008; Air et al., 2011; Bhanpuri et al., 2014).

Multilocus, multitarget approaches may be particularly useful in movement disorders because movement produces coordination by utilizing coordination among multiple brain areas including the basal ganglia, thalamus, cerebellum, sensory, and motor cortices (Neychev et al., 2008; Crowell et al., 2012; Delnooz and van de Warrenburg, 2012). Pathology within any one region, or disturbances in communication between any of the regions can potentially lead to disorders. To begin to address these multiple challenges, we focused our computer modeling here on a multiscale model of motor cortex and multitarget pharmacology based in this one area. In the future, this model will be expanded to encompass more areas and will include synaptic receptor targets in each area.

## AUTHOR CONTRIBUTIONS

All authors listed, have made substantial, direct and intellectual contribution to the work, and approved it for publication.

### ACKNOWLEDGMENTS

The authors would like to thank Ben Suter and Gordon MG Shepherd (Northwestern University) for help with the model; Tom Morse (Yale) for ModelDB support; R.A.N. for helpful suggestions. The authors declare no competing financial interests. Research supported by NIH grant R01 MH086638, NIH grant U01 EB017695, NIH grant R01 NS064046, NIH grant R01 DC012947. The NIH had no role in study design; in the collection, analysis and interpretation of data; in the writing of the report; and in the decision to submit the article for publication.

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fphar. 2016.00157

### REFERENCES


Kendall, M. (1938). A new measure of rank correlation. Biometrika 30, 81–93.


causing hypertonia in childhood. Pediatrics 111, e89–e97. doi: 10.1542/peds. 111.1.e89


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2016 Neymotin, Dura-Bernal, Lakatos, Sanger and Lytton. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# **Multi-target pharmacology: possibilities and limitations of the "skeleton key approach" from a medicinal chemist perspective**

*Alan Talevi\**

*Medicinal Chemistry, Department of Biological Sciences, Faculty of Exact Sciences, National University of La Plata, La Plata, Argentina*

#### *Edited by:*

*Thomas J. Anastasio, University of Illinois at Urbana-Champaign, USA*

#### *Reviewed by:*

*Joel D. A. Tyndall, University of Otago, New Zealand Maria G. Morgese, University of Foggia, Italy*

#### *\*Correspondence:*

*Alan Talevi, Medicinal Chemistry, Department of Biological Sciences, Faculty of Exact Sciences, National University of La Plata, 47 and 115 Street, La Plata B1900AVV, Argentina atalevi@biol.unlp.edu.ar, alantalevi@gmail.com*

#### *Specialty section:*

*This article was submitted to Experimental Pharmacology and Drug Discovery, a section of the journal Frontiers in Pharmacology*

> *Received: 20 June 2015 Accepted: 04 September 2015 Published: 22 September 2015*

#### *Citation:*

*Talevi A (2015) Multi-target pharmacology: possibilities and limitations of the "skeleton key approach" from a medicinal chemist perspective. Front. Pharmacol. 6:205. doi: 10.3389/fphar.2015.00205* Multi-target drugs have raised considerable interest in the last decade owing to their advantages in the treatment of complex diseases and health conditions linked to drug resistance issues. Prospective drug repositioning to treat comorbid conditions is an additional, overlooked application of multi-target ligands. While medicinal chemists usually rely on some version of the lock and key paradigm to design novel therapeutics, modern pharmacology recognizes that the mid- and long-term effects of a given drug on a biological system may depend not only on the specific ligand-target recognition events but also on the influence of the repeated administration of a drug on the cell gene signature. The design of multi-target agents usually imposes challenging restrictions on the topology or flexibility of the candidate drugs, which are briefly discussed in the present article. Finally, computational strategies to approach the identification of novel multi-target agents are overviewed.

**Keywords: multi-target agents, lock and key paradigm, gene profile, drug resistance, drug repositioning, drug design, designed multiple ligands**

### **Introduction**

Multi-target drugs (or multi-functional drugs or network therapeutics) have attracted considerable attention in the last decade, as potential therapeutic solutions to diseases of complex etiology (Talevi et al., 2012; Koerberle and Werz, 2014; Zheng et al., 2014) and health conditions linked to drugresistance issues (Talevi and Bruno-Blanch, 2013; Li et al., 2014). According to the "one drug, one target" paradigm, highly potent and specific (single-target) treatments would be better tolerated due to absence of off-target side-effects. However, poor correlation between *in vitro* drug effects and *in vivo* efficacy is often found with target-driven approximations (Kell, 2013; Margineanu, 2014). While target-first strategies might prove useful to approach single gene disorders, disease is often a multifactorial condition involving a combination of constitutive and/or environmental factors. Owing to compensatory mechanisms and redundant functions, biological systems are resilient to single-point perturbations (Hopkins, 2008). Under such perspective, disease often results from the breakdown of robust physiological systems due to multiple genetic and/or environmental factors, leading to the establishment of robust disease conditions (Yildrim et al., 2007). Thus, complex disorders are more likely to be healed or alleviated though simultaneous modulation of multiple targets.

Though this strategy has only been purposely applied in the last 10 to 15 years, many of the previously known therapeutic agents are in fact multi-target ligands (Yildrim et al., 2007),

which is especially true for those drugs that were discovered by serendipity, phenotypic screening or traditional medicine. Note that in all these cases, the knowledge on the pharmacological effect precedes the knowledge of the mode of action. Aspirin itself has been shown to act through a diversity of molecular mechanisms besides cyclooxygenase inhibition (Koerberle and Werz, 2014). Some therapeutic categories, e.g., mood disorder medications, are particularly abundant on classical examples of multi-target drugs (Roth et al., 2004). So actually, multi-target drugs have long been known and effectively used in the clinical practice but have majorly been found serendipitously or through phenotypic screening. What are the possibilities and limitations of tailored multi-target drugs?

## **Revisiting and Squeezing the Classical Lock and Key Paradigm**

Medicinal chemists usually resort to the traditional lock and key model to describe the interaction between a ligand and its molecular target (or an updated version of this paradigm that contemplates the ligand and target flexibility, such as the handin-glove analogy). The general idea is that the ligand (*the key*) and the target (*the lock*) should have complementary features to efficiently interact and trigger some biological response (*open the lock*). Frequently, different ligands can elicit a qualitatively similar response at a certain target. For different keys to activate the same lock alike they must share some common, essential arrangement of features (*the blade of the key*), which will be termed the *pharmacophore* (from the Greek, *what carries the medicine*). The remaining part of the key (*the bow*) may be indeed important, but less subject to structural restrictions (**Figure 1**).

A multi-target ligand might be conceived as a *skeleton or master key* capable of unlocking several locks. While selective non-selectivity might be of benefit, promiscuity (non-selective non-selectivity) might in contrast raise severe safety concerns and should be avoided. Why may a non-promiscuous ligand activate different targets? There are many possible answers to this question. First, it is frequent for a given ligand to act on several isoforms of the same protein. For instance, xilocaine (lidocaine) can produce anesthetic, antiarrhythmic and anticonvulsant effects by blocking the peripheral nervous system, heart and central nervous system sodium channels (Catterall, 2000). Alternatively, different members of a given biochemical pathway might share, to some extent, ligand specificity due to co-evolution. Finally, a ligand might display affinity to two or more unrelated targets by combining different pharmacophores in the same molecule (Morphy et al., 2004). Frequently, such combination of pharmacophores leads to molecules that are either enthalpically or entropically unfavorable, which conspires against the design of multi-target drugs, as will be later discussed in the correspondent section. This is metaphorically represented in **Figure 1**, through the awkward design of key number 3.

The contribution of biotechnology, however, has made very clear that the lock and key analogy can fall short to explain the effects of a drug on a biological system, particularly when medium- and long-term drug exposure (multiple-dose regimens) is required. After sustained exposure to a chemical agent the gene signature of a cell varies: some genes are upregulated while others are downregulated (e.g., owing to activation of nuclear receptors, compensatory mechanisms, etc.). Whereas in the past attention was directed to the direct interactions between the drug and its molecular target/s, now it is known that a more holistic perspective is needed to fully characterize the action of a drug on a biological system. For example, it has been reported that chronic administration of valproic acid and carbamazepine downregulates cytosolic phospholipase A2 and/or cyclooxygenase (with the consequent reduction of proinflammatory cytokines; Bosetti et al., 2003; Gherlardoni et al., 2004), an effect that may be involved in the effectiveness of these agents in epilepsy and bipolar disorder. The need for such holistic view is unequivocally expressed in the Connectivity Map, a publicly available resource meant to connect disease and small molecules through geneprofiles (Qu and Rajpal, 2012). The Connectivity Map stores gene expression profiles derived from the treatment of human cells cultured with a large number of drugs; when a disease signature is used as a query, it is expected that those drugs related to the disease by opposite expression changes (inverse similarity) will be potential treatments.

### **Possible Applications of Multi-Target Ligands**

Three main applications of multi-target agents in a therapeutic can be envisioned.

### **Complex Disorders**

Complex disorders are multi-factorial health conditions triggered by a number of intrinsic and/or environmental factors acting together on an organism. Among them we may mention mood disorders, neurodegenerative diseases, chronic inflammation or cancer. Despite the advances on the comprehension of the biological basis of these conditions and the huge investments made by the pharmaceutical sector, pharmaceutical solutions remain elusive. Although in some cases such disorders can be or are approached through combined therapies, multi-target ligands would present clear advantages, among them more predictive pharmacokinetics, better patient compliance, and reduced risk of drug interactions. There are several reviews available covering the potential of the multi-target approach in cancer (Petrelli and Giordano, 2008; Petrelli and Valabrega, 2009), Alzheimer's disease (Bajda et al., 2011; Dias and Viegas, 2014; Zheng et al., 2014), Parkinson's disease (Youdim et al., 2014), inflammation (Hwang et al., 2013), depression and other psychiatric disorders (Wong et al., 2008; Milan, 2014).

### **Drug Resistance**

Simultaneously impacting different targets could also be advantageous to approach individuals expressing intrinsic or induced variability in drug response due to modifications in key disease-relevant biological pathways and activation of compensatory mechanisms (Zimmermann et al., 2007; Xie et al., 2012). Apart from the obvious applications in the field of antimicrobial chemotherapy (it is less probable to develop resistance linked to single-point mutations against multi-target than single-target agents) this strategy could also be pertinent to treat non-infectious conditions characterized by high incidence of the drug resistance phenomena, e.g., epilepsy (Bianchi et al., 2009; Margineanu, 2014). One third of the epileptic patients suffer from refractory epilepsy. One of the prevalent hypotheses to explain refractory epilepsy cases proposes that at least part of the non-responsive patients might express variations in molecular targets of antiepileptic drugs (Talevi and Bruno-Blanch, 2013). Isobolographic studies in animal models and clinical experience suggest that combination of drugs with different mechanisms tends to be beneficial (Kwan and Brodie, 2006; Kaminski et al., 2009; Lee and Dworetzky, 2010; Brodie et al., 2011). On the other hand, while there exists consensus regarding the utility of single-target drugs for the treatment of some specific epilepsy types or syndromes, broad spectrum antiepileptic drugs such as valproic acid are among the most used antiepileptic agents and might be valuable in those cases where, at the onset of epilepsy, diagnosis of the specific syndrome is elusive (Bourgeois, 2007; Lagae, 2009; Löscher et al., 2013; Margineanu, 2014).

### **Prospective Drug Repositioning**

Drug repositioning (i.e., finding a second or further medical use for already known therapeutics, including approved, discontinued, shelved, and experimental drugs) has attracted enormous interest within the academic and pharmaceutical sectors during the last 10 years (Ashburn and Thor, 2004; Novac, 2013). Most of the successful drug repositioning cases have been found by serendipity or through exploitation of the original action mechanism of a drug for new indications (*on target* repositioning). Multi-target agents are natural candidates for more innovative, *off-target* drug repositioning. Computational approaches to drug repositioning have so far focused on what we will call *retrospective drug repositioning*: screening known drugs collections/libraries to find novel indications for already known therapeutic agents. *Prospective drug repositioning*, in contrast, would explore drug repositioning possibilities much earlier in the drug discovery process. While some pharmaceutical companies now consider exploring repositioning alternatives for drugs in the pipeline, the approach could be taken much further, by designing multi-purpose drugs to treat different conditions; prominently, frequently co-morbid disorders (e.g., diabetes and cardiac disease; anxiety and peptic ulcer disease, epilepsy, and depression) or, alternatively, underlying pathologies plus disease symptoms. The case of amiodarone and related compounds and Chagas disease can be illustrative. Chagas disease is a tropical parasitic disease historically endemic to Latin America. The late phase of the disease is characterized by life-threatening heart disorder in around one third of the patients. Amiodarone is a class III antiarrhythmic agent that shares many characteristics of other electrophysiological anti-arrhythmic drugs, including inhibition of sodium and potassium channels and L-type calcium channels. Interestingly, some studies showed that patients with chagasic cardiomyopathy treated with amiodarone had a more rapid recovery when compared with other patients treated with class I and class IV antiarrhythmics. This fact suggested that other mode of action could be in play. It was later demonstrated that amiodarone was able to act directly on the parasite survival, affecting the growth of *Trypanosoma cruzi* extracellular epimastigotes and *T. cruzi* amastigotes (that is, amiodarone could act on the underlying pathology). The mechanism of action of the drug was elucidated, showing that this drug directly disrupts the intracellular calcium regulation of the parasite (Benaim et al., 2006). Similar results were later observed with dronedarone (Benaim et al., 2012). Still, this example is another case of retrospective drug repositioning, since the new medical use emerged from clinical observations. A future challenge is to define whether this kind of indication expansion oriented to the treatment of co-morbid conditions could be anticipated through rational approaches at early stages of the drug development process, thus helping to provide evidence on possible advantages of new treatments compared to the existent ones, and additional criteria to decide which drug candidates should be prioritized to clinical trials and to conveniently choose the clinical endpoints of the trial that will be used to test superiority or non-superiority of the treatments under comparison. Computational networkbased approximations could prove valuable to unveil hidden connections between diseases and assist these types of initiatives.

### **Some Considerations Related to the Design and Screening of Multi-Target Agents**

Development of tailored multi-target agents with affinity to unrelated or weakly related drug targets relies mainly in two approaches (Morphy et al., 2004; Ma et al., 2010): the methodical combination of pharmacophores from selective, single-target ligands (a fragment-based approach) and; the screening of compound collections by simultaneous application of multiple computational models (or a single, multi-tasking computational model) to identify compounds with a suitable combination of activities. In the first approximation, the distinct pharmacophores are joined together by a cleavable or stable linker or, alternatively, they are overlapped by taking advantage of structural commonalities (Morphy et al., 2004). The use of linkers often leads to compounds with unfavorable biopharmaceutic or pharmacokinetic profile (e.g., compounds that violate more than two of the Lipinski's rules). Although the use of cleavable linkers might be advantageous, it also limits some of the merits of the multi-target approach in comparison with combination therapies (simplified pharmacokinetics, reduced chance of drug interactions). Moreover, the fragment-based approach could lead to poor ligand efficiency metrics (Hopkins et al., 2014), which refer to the binding efficiency *per atom*. It might be speculated that, since only a part of the molecule can interact with each of the proposed targets, the other part can become an obstacle for the binding event, reducing the binding efficiency because of enthalpic and/or entropic reasons, which is represented through the awkward topology of key number 3 (**Figure 1**). Therefore, the overlapping or merging approach (searching partially or highly integrated pharmacophores in a small molecule) seems more attractive from a biopharmaceutical viewpoint. Including some degree of flexibility in the molecule may help the common and non-common pharmacophoric features to accommodate to the correspondent binding sites of the different intended targets; however, the degree of flexibility should be carefully tuned so that an excess of flexibility does not conspire against the binding affinity (owing to unfavorable entropic loss associated to the binding event) or the bioavailability of the drug (it should be remembered that many druglikeness rules limit the number of flexible bonds in the molecule). An illustrative example of some of these principles is provided by the recent research from Jayaraman et al. (2013). These authors applied the fragment-based approach in the design of phytochemical-antibiotic conjugates conceived as multivalent inhibitors of *Pseudomonas aeruginosa* DNA gyrase subunit B (GyrB)/topoisomerase IV subunit B, dihydrofolate reductase (DHFR) and dihydropteroate synthase (DHPS). Departing from previously identified pharmacophores for inhibitors of *E. coli* GyrB and DHFR, the authors derived a common pharmacophoric model for multi-inhibition of such enzymes. Remarkably, they decided on using simple phenols (gallic acid and protocatechuic acid, simpler structural analogs of the bivalent natural product epigallocatechin gallate) conjugated through a non-cleavable linker to sulfamethoxazole and sulfadiazine (which inhibit DHPS; **Figure 2**). The decision of using simple phytochemicals as departure points resulted in four drug-like compounds with acceptable computed biopharmaceutical properties, which was checked through different drug-likeness rules (Lipinski and Veber rules) and by predicting the solubility and the percentage of absorption for the designed drug candidates. Two of the candidates displayed no violation of the rule of five and Veber rules, while the remaining two showed only one violation of Lipinski rules and marginal violation of Veber rules.

Regarding the screening approximation, one should bear in mind that the hit rate in the screening campaign is expected to be lower than the ones obtained when looking for single-target drug candidates (Talevi et al., 2012): each model used in the *in silico* screening process functions as a structural restriction that filters out all the molecules that do not gather the model requisites; thus, the more models used, the less probable it is to find chemical compounds accomplishing all the models structural constraints. For example, Nair et al. (2013) have recently performed a virtual screening campaign to identify multi-target inhibitors of DAP-kinases (a family of pro-apoptotic proteins also involve in autophagy, which are proposed as a promising target for therapeutic intervention of brain ischemia and neurodegenerative diseases). Among DAP-kinases, DRP1 has been reported to be the upstream protein of all the DAP-kinases as it is involved in the activation of other members of the family. However, modulation of DRP1 is not enough to attenuate the cell death pathways activated by DAP-kinases, owing to the existence of alternative activating sources. Searching for multi-target agents, the authors have explored a combined database of 391 known ligands of one of three members of the DAP-kinases family: DAPk1, DRP1, and ZIPk. This library was compiled from the Protein Data Bank

and ChEMBL, and it was sequentially screened through three pharmacophore hypothesis of DAPk1, DRP1, and ZIPk, in that order. Screening using the first hypothesis (DAPk1) resulted in 196 hits. Further, screening of these hits by DRP1 pharmacophore resulted in 56 hits, which contained pharmacophore features of both DAPk1 and DRP1 ligands. The 56 ligand hits were then screened by the ZIPk pharmacophore, retrieving only four ligands gathering the pharmacophoric features of all three DAP-kinases. The limited number of hits obtained when using sequential *in silico* filters/models to select multi-target agents might be compensated by the huge, ever expanding available chemical universe. Multitasking QSAR approximations (Zanni et al., 2014; Speck-Planche and Cordeiro, 2015) could prove as a valuable tool to implement this strategy.

## **Target Selection**

So far the advantages and challenges posed by the multi-target approach have been discussed. A critical question remains, however, to be made: if we are to design multi-target agents, how shall we choose our molecular targets? Obviously, a drug target needs to have the potential to be disease modifying. Secondly, if we are fighting against an infection or a deregulated cell (e.g., in cancer) the drug must display some degree of selectivity, e.g., the drug target must be exclusively or preferentially expressed in the infectious agent or in the cancerous cell, targeted proteins in a pathogen should not have homologous proteins in the host or homologous proteins in the host should be sufficiently different from those in the pathogen, etc. Furthermore, the Medicinal Chemistry community has long accepted that not all the proteins are equally "druggable," i.e., likely to be moderated by small molecules. A number of approaches to assess druggability have been proposed in the specialized literature, from "guilt by association" approximations (a protein is predicted to be druggable if it belongs to a protein family for which at least one member of the family is targeted by a drug) to methods based on binding site prediction, among others (Keller et al., 2006; Cheng et al., 2007). But still the previous are just general considerations valid for both single- and multi-target approximations. When aiming at multiple targets, the choice of the targets and the pursued type of inhibition depend on several factors, among them the nature of the disease (infectious disease? complex disorder?) and/or the possible mechanisms of drug resistance (adaptive mechanisms? target amplification or mutation?). A relevant issue that deserves attention is whether it is preferable to directly block the selected targets or to modulate them (e.g., through weak partial inhibitions). Under our modern paradigm, built on a systems biology perspective, it is understood that, in general, we are not targeting isolated proteins but pathways instead. We might target different signaling pathways (parallel targeting), which may be valuable to block escape routes, adaptive resistance mechanism and compensatory homeostatic responses; alternatively, vertical targeting (attacking the same pathway at different nodes) might prove useful against other types of resistance (e.g., target mutations; Shahbazian et al., 2012). When trying to kill pathogens or malignant cells, attacking hubs (highly connected nodes in a biochemical network) might be the strategy of choice; on the other hand, if the treatment objective is to restore a perturbed network to a healthy state, using low affinity multi-target ligands to modulate multiple non-crucial nodes neighboring key nodes ligands might be advantageous in order to avoid sever side-effects (that might be otherwise expected if blocking a key node with a crucial physiological function; Csernely et al., 2013). Metabolic control analysis constitutes a useful frame to evaluate the importance and relative contribution of individual metabolic steps in the overall functioning of a particular system and, subsequently, to identify optimal targets (Hornberg et al., 2007).

### **Conclusion**

Multi-target agents are a promising strategy to face complex, multifactor disorders and drug resistance issues. Additionally, they can prove valuable in prospective drug repositioning oriented to the treatment of comorbid conditions or both the underlying pathology and its symptoms, an overlooked application to the moment. Compared to combination therapies, they present several advantages, including more predictable pharmacokinetics, lower probabilities of drug interactions and higher patient compliance.

### **References**


We have highlighted some difficulties related to the search of tailored multi-target drugs (e.g., enthalpic and entropic considerations and potential bioavailability issues, limited number of hits when sequentially screening a virtual library). Besides the classical key and lock paradigm to approach the multi-target strategy, the effect of the drug on cell gene signatures should also be considered, especially when looking for middleand long-term treatments, which is often the case for complex disorders. Finally, network analysis might provide clues to help target selection, which is highly dependent on the nature of the treated disorder and the known mechanisms of resistance.

### **Acknowledgments**

The author would like to thank National University of La Plata (X729, X730), ANPCyT (PICT 2013-0520, PICT 2013-3175, PPL 2011/0003) and CONICET (PIP l2201201 00324 CO) for funding. The author is a member of CONICET.

expression and cyclooxygenase activity in rat brain. *Biol. Psychiatry* 56, 248–254. doi: 10.1016/j.biopsych.2004.05.012


approaches to different hypothesis," in *Pharmacoresistance in Epilepsy: From Genes and Molecules to Promising Therapies*, eds L. Rocha and E. A. Cavalheiro (New York: Springer), 207–224.


**Conflict of Interest Statement:** The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2015 Talevi. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Computational polypharmacology comes of age

Giulio Rastelli\* and Luca Pinzi

Molecular Modelling and Drug Design Lab, Life Sciences Department, University of Modena and Reggio Emilia, Modena, Italy

Keywords: polypharmacology, multitarget ligands, drug discovery, drug design, molecular modeling

In the last years, the "one target, one drug" paradigm that has traditionally dominated drug discovery has been deeply challenged by the evidence that small molecules interact simultaneously with multiple targets, a phenomenon known as polypharmacology. Today, polypharmacology is recognized as a new valuable opportunity for drug discovery and development. It is now well established that drug molecules typically bind to several targets, and that their efficacy and safety is mostly dependent on their polypharmacological profile (Jalencas and Mestres, 2012; Peters, 2013; Anighoro et al., 2014). Indeed, one of the most common reasons for terminating a drug discovery program has been promiscuity or lack of selectivity of the developed compounds. This leads to important considerations regarding the polypharmacology inherent in chemical structures and its possible exploitation for drug discovery. First, side effects caused by drug binding to unwanted off-targets (adverse polypharmacology) should be identified as early as possible in the drug discovery pipeline. Second, potential synergistic effects arising from hitting multiple targets (beneficial polypharmacology) should be taken into consideration and thoroughly incorporated in the drug design strategy. Third, polypharmacological approaches have the potential to redirect stalled drug discovery projects and to reposition valuable hits or leads (drug repositioning). Finally, prediction of polypharmacological profiles can be used to uncover new macromolecular targets for already known or new developing drugs (target identification and deconvolution). In all these areas, computational polypharmacology is gaining a foothold in drug discovery, as witnessed by the increasing number of publications reporting theoretical approaches and methods specifically put forward to address these needs.

State-of-the-art computational approaches offer the possibility to predict the activity profile of ligands to a set of targets, thereby anticipating potential selectivity issues or discovering desired multitarget activities early in the iterative design and optimization steps typical of a preclinical drug discovery project. These approaches stem from 2D or 3D shape and chemical similarity, pharmacophore analyses, target and binding site similarity assessment, docking methods, bioinformatics, graph theory and modeling, machine-learning algorithms, and chemogenomics (**Figure 1**). Broadly, these can be classified into statistical data analysis and bioinformatics, ligandbased, and structure-based approaches, all of which are well documented in the literature (Csermely et al., 2005; Boran and Iyengar, 2010; Bottegoni et al., 2012; Anighoro et al., 2014; Reddy et al., 2014). One should note that ligand-based and structure-based strategies have specific advantages and limitations. Structure-based methods use the information derived from knowledge of the 3D structure of proteins. These methods are applicable to identify ligands for a specific target or set of targets of interest, for example by performing de-novo design or virtual screening of large libraries of small molecules. In addition, they can be used to assess binding site structural similarity and to profile protein-ligand interactions among sets of targets. Their application is obviously limited to proteins with known crystal structure or to homology models derived from highly homologous crystal structure templates. Moreover, structure-based results are influenced by differences in conformations of binding site residues, which are generally difficult to predict. Ligand-based approaches do not require crystal structures of the target proteins but rely on

#### Edited by:

Thomas J. Anastasio, University of Illinois at Urbana-Champaign, USA

Reviewed by: Alessandro Giuliani, Istituto Superiore di Sanità, Italy

> \*Correspondence: Giulio Rastelli, giulio.rastelli@unimore.it

#### Specialty section:

This article was submitted to Experimental Pharmacology and Drug Discovery, a section of the journal Frontiers in Pharmacology

> Received: 17 June 2015 Accepted: 14 July 2015 Published: 28 July 2015

#### Citation:

Rastelli G and Pinzi L (2015) Computational polypharmacology comes of age. Front. Pharmacol. 6:157. doi: 10.3389/fphar.2015.00157

prior knowledge of biologically active ligands, therefore their use is limited to targets for which ligands are known. Worth of note is that in ligand-based methods the derived information is necessarily dependent on the chemical structures of the classes of compounds that have been thus far developed. As a consequence, predicting polypharmacological profiles of ligands that are too dissimilar to already synthesized classes of compounds would be impossible. Overall, ligand-based and structure-based methods appear to be applicable in conjunction to provide more robust results (Anighoro et al., 2015). Such combination offers the possibility to take advantage of the peculiar features and strengths of each approach toward the obtainment of possible candidates for polypharmacology, and appears to be one promising way to go in future investigations. For example, one risk of predicting polypharmacology by using only chemical similarity principles is that inactive compounds can exhibit high similarity with active molecules if they derive from a slight modification of an active compound at some key position crucial for its interaction with the target. In this case such similarity would lead to false positives. Likewise, false negatives can be expected considering that not all of active compounds have been identified for a given target. In these cases, structure-based methods can help overcome these potential pitfalls by estimating the steric and electrostatic complementarity of ligands with the target binding sites. For example, structure-based docking screenings of compounds that passed the desired chemical similarity filters may be independently performed on two or more biological targets of interest, and multi-target hits may be identified from compounds located at the top of all ranked lists. Finally, analysis of drug targets and drug-target associations using a network approach may provide useful information to highlight particularly interesting target combinations or chemical modulators able to perturb the network at specific nodes of disease-specific critical pathways (Csermely et al., 2013). In this context, partial inhibition of a small number of targets can be more efficient than the complete inhibition of a single target, especially for complex and multifactorial diseases (Csermely et al., 2005). This information can be used by ligand-based or structure-based methods to direct the design and screening of new drugs toward the desired set of multiple drug targets.

Polypharmacology has been mainly recognized within members of the kinome and GPCRs families (Knight et al., 2010; Jacobson et al., 2014). This is not surprising, considering that binding sites within members of conserved and evolutionarily related targets are generally conserved and thus prone to multitarget inhibition. However, one should note that the recognized specificity of a ligand or a series of ligands depends heavily on how hard on- and off-targets have been investigated, and this is surely the case for kinases and GPCRs, which have been extensively explored. We are far from having the capacity to perform an exhaustive biological profiling of ligands that enter into drug discovery pipelines, but we can expect that the more testing will be performed, the more off-targets and multitarget activities will be seen also for targets genetically and structurally unrelated to the primary intended target. In this respect, constant improvement and implementation of compound and bioactivity data deposited in publicly available databases will provide access to an increasing number of high confidence bioactivity annotations for larger sets of chemicals and therapeutic targets (Hu and Bajorath, 2013). Overall, this information will be very useful to computationally design multitarget ligands. In parallel, improvement in hardware and software performance is making it possible to handle an enormous amount of data, thus enabling the generation and analysis of big data for polypharmacology in a very cost- and time-effective way.

The rational design of molecules interacting with more than one biological target becomes most challenging when these targets are only distantly related or unrelated, i.e., when they belong to different protein families. For example, the selectivity of particularly interesting kinase inhibitors are usually profiled against a large panel of kinases of the kinome, but they are rarely screened against targets of other families due to limited capacity of experimental in vitro testing. Considering that local binding site similarities may be more important than global structural similarity to determine polypharmacological activities, especially when ligands are able to interact with key residues of more than one target, this remains a critical point for the development of multi-targeted drugs (Salentin et al., 2014). Therefore, assessing local binding site similarities and comparing protein-ligand interaction profiles, especially for distantly related sets of targets and chemical classes of ligands, will be crucial for predicting polypharmacology (both beneficial and harmful). Significant improvements are also needed on how to select the most relevant set of therapeutically important targets for a given disease, a question that can benefit of the recent progresses of proteomics and clinical molecular investigations on patients and disease states.

Progresses in modeling protein-ligand interactions and in quantitatively predicting free energies of binding of ligands to target proteins will definitely contribute to successful design of molecules with the desired polypharmacological profile. The ongoing advances in docking methods such as improvement of scoring functions and better treatment of receptor flexibility are playing an important role to meet this goal. Importantly, several free-energy based approaches, with different theoretical backgrounds and at different levels of approximation, have been proposed to rescore docking results in order to increase the accuracy of binding affinity predictions (Parenti and Rastelli, 2012). Considering that the affinity of a ligand for a target protein reflects the 1G of binding, any further improvement in our ability to accurately predict binding free energies will be important to design multi-target drug candidates.

Combining computational design and chemical synthesis of libraries of multi-target ligands provides another means to more effectively obtain bioactive compounds with the desired on- and off-target binding. For example, Reutlinger et al. very recently described the development and application of a computational molecular de novo method for designing combinatorial libraries that exhibit an accurately predicted bioactivity profile, obtaining nanomolar multitarget ligands modulating the dopamine D4 and sigma-1 receptors (Reutlinger et al., 2014). In another study, Rodrigues et al. showed that the combination of machinelearning methods with automated chemical synthesis and fast bioassay turnover enabled the generation of small molecules with the desired polypharmacology (Rodrigues et al., 2015). These investigations suggest that a combination of the two approaches may be suitable for rapidly obtaining hits and leads with the desired target engagement.

Finally, a thorough understanding of drug-target network relationships and target-disease associations is key not only to provide more effective and safer drugs, but also to uncover specific target combinations that may provide synergistic effects and/or benefits for mitigating or bypassing drug resistance. In other words, selecting the "right" combination of targets for a specific disease will probably be a major key to success, and this should be given full consideration by focusing computational experiments on target combinations suggested by clinical and/or molecular biology investigations. So far, given the high number of cellular targets and our limited ability to understand their interplay in disease states, most biologically active small molecules are likely to bind several targets and/or to activate or suppress alternate pathways or targets. Network models are providing useful information to analyze the interconnection of pathways and targets relevant to human diseases, and their relation with chemical compound networks (Schadt et al., 2009). However, the more the pathways and mechanisms of disease (especially multifactorial and complex diseases) will be understood at the molecular level, the more the polypharmacological networks can be exploited with computational methods to obtain safer and potent drugs able to modulate the desired on- and off-target activities. The recent successes in de novo predicting drug polypharmacology and the raising number of computational strategies and frameworks developed at this purpose testify that computational polypharmacology has come of age and will play an increasingly important role in drug discovery. The combination of different approaches and expertise (experimental and computational) will likely be key to success.

### References


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2015 Rastelli and Pinzi. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Developing multi-target therapeutics to fine-tune the evolutionary dynamics of the cancer ecosystem

Lei Xie1, 2 \* and Philip E. Bourne<sup>3</sup>

*<sup>1</sup> Department of Computer Science, Hunter College, The City University of New York, New York, NY, USA, <sup>2</sup> The Graduate Center, The City University of New York, New York, NY, USA, <sup>3</sup> Office of the Director, National Institutes of Health, Bethesda, MD, USA*

Keywords: multi-target drug, cancer evolution, non-linear dynamic system, polypharmacology, cell-cell communication

Multi-target therapies, either in combination or in sequential order, have been advocated to combat intrinsic and acquired resistance to anti-cancer drugs (Holohan et al., 2013; Yardley, 2013). However, the effectiveness of multi-target anti-cancer therapy in the clinic is limited. The selection of cancer cells obeys Darwin's law of evolution. Under the pressure of drug perturbation, the cancer cell can adapt versatile molecular and cellular mechanisms for survival, and often evolves into more aggressive or metastasis phenotypes (Holohan et al., 2013). At the molecular level, acquired mutations resulting from drug treatment may modify drug metabolism (e.g., increasing efflux, decreasing uptake, and enhancing detoxification etc.) and alter drug-target interactions. At the cellular level, multiple pathways support the survival of cancer cells. The inhibition of one pathway may result in the activation of an alternative pathway. Although novel approaches to optimizing combination therapies have been proposed to defer these drug resistance mechanisms (Crystal et al., 2014; Wang et al., 2015), intra-tumor heterogeneity that have been observed ubiquitously may make the drug combination fail (McGranahan and Swanton, 2015). Polygenic drug-resistance mechanisms are present in sub-clones prior to the initiation of therapy (Bozic et al., 2013). If the therapy cannot target all sub-clones that drive the cancer progress in a fast-killing mode, it would prompt the rapid growth of sub-clones that are not sensitive to the treatment (Gatenby et al., 2009). Unfortunately, the number of driver mutations in advanced tumors is substantial (Gerlinger et al., 2014). It could be an impossible mission to target all driver mutations. The existence of cancer stem cells adds another dimension of complexity. The conventional single or combinational anti-cancer drug is incapable of killing cancer stem cells. When a cancer cell is killed by chemotherapy, it could send signals to stimulate the proliferation of the cancer stem cell, leading to the repopulation of the drug-resistance tumor (Kurtova et al., 2015). Thus, new strategies are needed to combat anti-cancer drug resistance with the goal to improve the effectiveness of anti-cancer therapy.

Cancer cells originate from the host's normal cells, but eventually turn into a new "pathogen" species. To eliminate anti-cancer drug resistance, we may borrow an idea from anti-bacterial drug discovery. Anti-virulence has emerged as a novel concept in addressing the challenge of antibiotic resistance (Rasko and Sperandio, 2010). Instead of killing bacteria, the anti-virulence drug interferes with the bacterial virulence and/or cell-cell communication, disrupts pathogenhuman interactions, or enhances the host's inner immunity. The rationale is that the bacterium is less likely to evolve into a drug-resistant strain when facing less evolutionary pressure. As the cancer adapts similar mechanisms to the bacterium in acquiring drug resistance, the "anti-virulence" strategy could be applied as an anti-cancer therapy. Recent successes in anti-cancer immune therapy open a new door to exploring the "anti-virulence" strategy in cancer treatment (Johnson et al., 2015). If the cancer can be controlled as a less aggressive or non-metastasis type, it may be possible to cure the cancer by boosting anticancer immunity (Carmi et al., 2015). On the contrary,

#### Edited by:

*Thomas J. Anastasio, University of Illinois at Urbana-Champaign, USA*

#### Reviewed by:

*Masaaki Murakami, Hokkaido University, Japan Jeff M. P. Holly, University of Bristol, UK*

> \*Correspondence: *Lei Xie, lei.xie@hunter.cuny.edu*

#### Specialty section:

*This article was submitted to Experimental Pharmacology and Drug Discovery, a section of the journal Frontiers in Pharmacology*

> Received: *18 June 2015* Accepted: *07 September 2015* Published: *24 September 2015*

#### Citation:

*Xie L and Bourne PE (2015) Developing multi-target therapeutics to fine-tune the evolutionary dynamics of the cancer ecosystem. Front. Pharmacol. 6:209. doi: 10.3389/fphar.2015.00209* the chemotherapy may stimulate the production of immunosuppressive molecules (Shalapour et al., 2015). As a result, the patient's anti-cancer immune response is inhibited. Along these lines, adaptive therapy has been proposed to control the tumor growth by permitting the survival of drug sensitive cells. In this way, the growth of drug-resistant clones could be surpassed (Gatenby et al., 2009). In some cases, the cancer can be treated as a chronic disease when transferring the cancer cell into a quiescent state (Aguirre-Ghiso, 2006). For example, chronic myeloid leukemia (CML), chronic lymphocytic leukemia (CLL), and low grade non-Hodgkin's lymphoma (NHL) are slow-growing cancers. Patients can live with them for many years.

The various heterogeneous types of cancer cells form an ecosystem, cooperating and competing with each other for nutrients and spaces from the harsh environment. For example, sustained angiogenesis, one of the hallmarks of cancer, relies on the cooperation of co-existing cell lineages (Floor et al., 2012). The cooperating sub-clones either bear complementary traits or play a different role of producer or consumer of "public goods" such as diffusible growth factors (Korolev et al., 2014). Based on theoretical, experimental, and clinical results of ecology, microbiology, and cancer research, it has been proposed that tuning the population dynamics of cancer cells can be a powerful strategy in developing an anti-cancer therapy (Korolev et al., 2014). By either changing the tumor microenvironment, or confusing cancer cell-cell communications, the whole cancer ecosystem can be controlled, even eliminated.

To determine the evolutionary dynamics of the cancer ecosystem, and the drug targets that can modulate its evolutionary trajectory, we need a deep understanding of not only the drug response and resulting evolution of individual cell types but also the emergent properties of the whole system under a diverse genetic and environmental background, which is more than a simple summarization of the behavior of all cells. Recent advances in cancer biology, single cell technologies, next-generation sequencing, and systems biology provide great opportunities to dissect the evolution of the cancer ecosystem at multiple scales. Multiple molecular components such as integrin and cadherin and pathways (e.g., Rho GTPase), which are responsible for the cell-cell communications and the cellenvironment interactions, have been revealed (Brücher and Jamall, 2014). They represent potential drug targets to perturb cancer cell-cell interactions and to constraint tumor growth. The Cancer Genome Atlas has identified millions of somatic mutations (Ledford, 2015). Correlated with generic variations, drug response phenomics data are available at molecular, cellular, tissue, and organism levels (Zbuk and Eng, 2007). The systematic integration of these data may allow us to predict drugresponse phenotypes for intervention against multiple targets and pathways. Single cell sequencing has revealed the clonal evolution of breast cancer (Wang et al., 2014) and childhood acute lymphoblastic leukemia (Gawad et al., 2014), and drug resistance dynamics (Lee et al., 2014). These studies will provide critical information to predict the cancer evolutionary trajectory, leading to the development of pre-emptive treatment strategies that are in contrast to current reactive clinical approaches. In spite of this progress, one of fundamental challenges remaining is to bridge genetic and molecular mechanisms of single cell-cell interactions to the ecological dynamics of the cancer population.

Multi-scale modeling and simulation may play a key role in predicting the evolutionary dynamics of the cancer ecosystem, and identify anti-cancer therapeutic targets for pre-emptive treatment. It is possible to reconstruct context-specific whole cell models by integrating multiple omics data (Karr et al., 2012). Subsequently, their cellular functions can be simulated and predicted at different evolutionary stages under a framework of constraint-based modeling (Bordbar et al., 2014). Using a single cell or sub-clone as the building block, the cancer ecosystem can be modeled as a dynamic cell-cell interaction network, in which the node is a cell, and the edge represents the cellcell interaction. Each node, or cluster of nodes, has different traits and evolutionary trajectory, yet depend on each other, as shown in **Figure 1**. A network representation may allow us to understand the emergent properties of the cancer ecosystem. For example, one of intrinsic properties of biological network is "robust-yet-fragile" (Kitano, 2007). The removal of a single node may have little impact on the whole system. However, the weak perturbation of multiple nodes can lead to the system failure, even if these nodes are not deleted. A number of effective therapies in treating complex diseases may follow this principle (Xie et al., 2012). For instance, successful anti-psychotic drugs, such as clozapine, mediate their effects through binding entire families of serotonin and dopamine receptors. The clinical failures of many anti-psychotic drugs can be attributed to them being too selective as designed (Hopkins et al., 2006). In another example, the anti-cancer effect of HIV protease inhibitors is proposed to comes from their weak bindings to multiple kinases (Xie et al., 2011). Moreover, the perturbation of edges may be more effective than nodes to regulate the state transition of a non-linear dynamic system (Tong et al., 2012). An additional advantage of edge perturbation is that the cancer cell has little selection pressure to evolve into a drug resistant phenotype, as the cancer cell will not be killed directly by the drug. In a proof-of-concept study, blocking cell-cell communication inhibited the repopulation of cancer stem cells, thus enhancing the effectiveness of anti-cancer therapy (Kurtova et al., 2015). To capture the whole dynamic spectrum of the cancer ecosystem, mechanistic and quantitative dynamic simulation is needed. Coarse-grained dynamic modeling has successfully identified an optimized sequence of therapies to improve the survival of patients with metastatic castrate resistant prostate cancer (Gallaher et al., 2014). Agent-based modeling that has been successfully applied to study the dynamics of complex systems could be a powerful tool to integrate whole cell models into a dynamic model of the cancer ecosystem (An et al., 2009).

Multi-target therapy can be achieved by either polypharmacology or drug combination. Polypharmacology has several advantages over drug combination. Firstly, it is not a trivial task to optimize dosages and sequences of a drug combination. The highly heterogeneous nature of cancer cells makes optimization even more difficult. Secondly, potential drug-drug interactions may cause serious side-effects. A singe "dirty" drug may reduce the probability of this problem. Finally,

it is argued that polypharmacology is more likely to achieve desired selective profiles than the drug combination (Varshavsky, 1998). Compared with anti-virulence agents for bacteria, selectivity is particularly challenging in the development of anti-cancer therapy, as the normal cell is more similar to the cancer cell than the bacteria. The side-effect of anti-cancer therapy is mainly because the drug cannot distinguish the normal cell from the cancer cell. The cancer genes that harbor driver mutations can be either up-regulated or down-regulated. A polypharmacological agent can be designed in such a way that it is mutually exclusive—it binds both the up-regulated gene A and the down-regulated gene B (terms A+ and B-, respectively) (Varshavsky, 1998). Consequently, the agent will selectively bind to A+ in the cancer cell as B- is less competitive. In the normal cell, A and B will competitively bind the agent, thus the agent will have little impact on the function of either of them. In contrast, the combination of two drugs, which bind to A and B, respectively, will be less selective. Following the same principle above, a selective agent can be designed for genes that are all up-regulated through targeting the allosteric site of one gene (Varshavsky, 1998).

The identification of suitable therapeutic targets is only the starting point of anti-cancer drug development. It is often more challenging to discover molecules that are able to achieve desirable therapeutic effects with minimum side effects. In addition to conventional druggable targets, targeting the protein-protein interaction (PPI) interface could be an efficient strategy to modulate signal transduction, cell-cell communications, and cell-environmental interactions (Jubb et al., 2012; Mullard, 2012). For example, cadherin-p120 interactions may mediate the contact inhibition of locomotion (CIL) that controls cell growth. The loss of CIL leads to tumor growth and metastasis (Mayor and Carmona-Fontaine, 2010). Thus, enhancing cancer cell-cell interactions may inhibit tumor progression. Historically, it is difficult to design small molecule inhibitors to target the PPI interface, as it is flat, large, and featureless, and thus considered undruggable. Multiple weak binders instead of a single strong binder could be an alternative strategy to designing PPI modulators, as the PPI interface is characterized by "hot spots" that contribute the most to the binding free energy, and well-defined grooves or small pockets that are associated with a continuous epitope binding partner (Ma and Nussinov, 2014). As mentioned previously, multiple weak binders may be more effective in regulating biological systems than a single selective inhibitor with high binding affinity. Although drug combinations are a successful multi-target therapy strategy, possible drug-drug interactions may limit the number of drugs administrated together. Polypharmacology offers two alternatives to the drug combination. Polypharmacology aims to design "dirty" drugs that can bind to multiple receptors simultaneously. Its effectiveness in treating systematic diseases has been documented (Xie et al., 2012). For example, targeted polypharmacological agents have been successfully designed to modulate signaling transduction events (Apsel et al., 2008). It is possible for a small molecule to target multiple PPI interfaces, as they are promiscuous across fold space (Gao and Skolnick, 2010), and may share conserved hot spot residues (Keskin et al., 2005; Shulman-Peleg et al., 2007). The desired multi-target binding specificity can be achieved through the fine-tuned side chain geometry and chemistry of protein-ligand interactions (Gao and Skolnick, 2010).

In summary, with advances in whole genome sequencing, high-throughput techniques, systems biology, and cloud computing, information on the evolutionary dynamics of cancer, from a single cell to the whole population, is starting to emerge.

### References


Concurrently, progress in medicinal chemistry is expanding the druggable target space (e.g., through targeting the proteinprotein interaction interface and allosteric events). Putting these efforts together may open a new door to multi-target therapies for cancer treatment.

### Acknowledgments

This research was supported, in part, under National Institute of Health Grant R01LM011986. We appreciate the constructive comments and suggestions by the editor and reviewers.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2015 Xie and Bourne. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# What is synergy? The Saariselkä agreement revisited

### Jing Tang\*, Krister Wennerberg\* and Tero Aittokallio\*

Institute for Molecular Medicine Finland (FIMM), University of Helsinki, Helsinki, Finland

Many biological or chemical agents when combined interact with each other and produce a synergistic response that cannot be predicted based on the single agent responses alone. However, depending on the postulated null hypothesis of non-interaction, one may end up in different interpretations of synergy. Two popular reference models for null hypothesis include the Bliss independence model and the Loewe additivity model, each of which is formulated from different perspectives. During the last century, there has been an intensive debate on the suitability of these synergy models, both of which are theoretically justified and also in practice supported by different schools of scientists. More than 20 years ago, there was a community effort to make a consensus on the terminology one should use when claiming synergy. The agreement was formulated at a conference held in Saariselkä, Finland in 1992, stating that one should use the terms Bliss synergy or Loewe synergy to avoid ambiguity in the underlying models. We review the theoretical relationships between these models and argue that one should combine the advantages of both models to provide a more consistent definition of synergy and antagonism.

Keywords: definition of synergy, drug combinations, Bliss and Loewe models, interaction barometer, consensus agreement

### Introduction

Evaluation of interaction effects between biologically active agents has become an important topic in many disciplines, including pharmacology (Cokol et al., 2011; Miller et al., 2013; Chevereau and Bollenbach, 2015), biochemistry (Hu et al., 2011; Zhang and Viikari, 2014; Bunterngsook et al., 2015), and environmental sciences (Darling and Côté, 2008; Piggott et al., 2015). The interaction between multiple agents is often classified as either synergistic or antagonistic, depending on how much the observed combination response differs from the expected response under the null hypothesis that the two agents are non-interacting. Multiple reference models have been formulated based on a distinctive set of empirical or biological assumptions (see e.g., Lehár et al., 2009). These assumptions, albeit difficult to validate a priori due to the lack of precise knowledge of the mechanisms of action, are often justifiable as long as they provide biologically plausible reasoning about the nature of non-interaction. However, the inherent differences in the model assumptions have inevitably led to inconsistency in the quantification of the degree of interaction, contributing to a major source of confusion and controversy on the definitions of synergy and antagonism. The Saariselkä agreement, proposed more than 20 years ago, aimed at reaching a consensus on the terminology for characterizing the degree of interaction (Greco et al., 1992). Acknowledging the theoretical background of the major competing models, the Saariselkä agreement admitted that there is no single universally best reference model. Rather than continuing the debate on the

#### Edited by:

Thomas J. Anastasio, University of Illinois at Urbana-Champaign, USA

#### Reviewed by:

Jiang Liu, University of Southern California, USA Murat Cokol, Sabanci University, Turkey

#### \*Correspondence:

Jing Tang, Krister Wennerberg and Tero Aittokallio, Institute for Molecular Medicine Finland (FIMM), University of Helsinki, Tukholmankatu 8, FI-00290 Helsinki, Finland jing.tang@helsinki.fi; krister.wennerberg@helsinki.fi; tero.aittokallio@helsinki.fi

#### Specialty section:

This article was submitted to Experimental Pharmacology and Drug Discovery, a section of the journal Frontiers in Pharmacology

> Received: 29 June 2015 Accepted: 11 August 2015 Published: 01 September 2015

#### Citation:

Tang J, Wennerberg K and Aittokallio T (2015) What is synergy? The Saariselkä agreement revisited. Front. Pharmacol. 6:181. doi: 10.3389/fphar.2015.00181 appropriateness of the model assumptions, the agreement called for a compromise between the advocates of the models, and proposed a practical guideline for reporting synergy or antagonism, where the underlying reference model should be explicitly described to avoid any ambiguity.

### The Two Reference Models

For the rest of the review, we will use drug combination as an example of interaction data. A drug's effect y is often measured at a certain dose x as the percentage of biological response, i.e., x > 0, 0 < y < 1. Let us consider that, drug 1 at dose x<sup>1</sup> produces a response y<sup>1</sup> and drug 2 at dose x<sup>2</sup> produces a response y2. Next, we combine the two drugs at the dose pair (x1, x2) and observe a combination response y<sup>c</sup> . To quantify the degree of drug interaction, we need to formulate a reference model to answer the following question: if there is no interaction between the drugs, what would be the expected combination response yEXP at (x1, x2)? So far, two major reference model classes have been proposed, the Bliss independence model (Bliss, 1939) and the Loewe additivity model (Loewe, 1953). While each having its own logical basis, the underlying assumptions behind these two models are relatively distinct.

The Bliss independence model adopts a probabilistic perspective by treating a drug combination under noninteraction as a joint action of independent, yet competing perturbations by the individual drugs. Such a probabilistic independence allows the expected combination response to be computed as the product of the individual drug responses:

$$\mathbf{y}\_{\text{BESS}} = \mathbf{y}\_1 + (\mathbf{1} - \mathbf{y}\_1)\mathbf{y}\_2 = \mathbf{y}\_1 + \mathbf{y}\_2 - \mathbf{y}\_1\mathbf{y}\_2. \tag{1}$$

An observed combination response greater or smaller than yBLISS can be interpreted as a departure from the probabilistic independence, which thus implies an interaction between the two drugs. The Loewe additivity model, on the other hand, requires additional information about the dose-response relationships of the individual drugs. Namely, let y = f1(x) and y = f2(x) be the dose-response functions for drug 1 and drug 2, respectively. Then the doses at which each drug alone produces the expected response yLOEWE can be represented as, f −1 1 (yLOEWE) and f −1 2 (yLOEWE), where f −1 is an inverse function which maps the response y back to the dose x. Formally, the Loewe additivity model states that yLOEWE must satisfy:

$$\frac{\varkappa\_1}{f\_1^{-1}(\text{y\_{LOEWE}})} + \frac{\varkappa\_2}{f\_2^{-1}(\text{y\_{LOEWE}})} = 1. \tag{2}$$

The rationale behind Equation (2) is to fit non-interaction to the so-called sham experiment scenario, where a drug is combined with itself, that is, f1(x) = f2(x). According to Equation (2), one can derive yLOEWE = f(x<sup>1</sup> + x2) for the sham experiment, reflecting the intuition that combining two drugs of the same type should induce neither synergy nor antagonism.

### The Saariselkä Agreement

The assumptions and performance of the two reference models have been compared and discussed in many review articles (e.g., Berenbaum, 1989; Greco et al., 1995; Chou, 2006; Lee, 2010; Zhao et al., 2010). There have been attempts to distinguish the Bliss and Loewe models in terms of mechanistic implications (Shafer et al., 2008; Laskey and Siliciano, 2014; Chevereau and Bollenbach, 2015). The Bliss independence model is expected to hold true for non-interacting drugs that elicit their responses independently, e.g., by targeting separate pathways. Loewe additivity, in contrast, is more compatible with the case where the drugs have similar modes of action on the same targets or pathways. However, little is known about whether such mechanistic justifications for the Bliss and Loewe models reflect the reality. Further, with increasing understanding of drugs' modes of action, any "previously unexpected" interaction effect becomes more expected, which makes the reference models totally dependent on the temporal state of knowledge. As pointed out in the Saariselkä agreement (Greco et al., 1992), and also by many others, neither Loewe additivity nor Bliss independence is necessarily reflecting the expected modes of action of a drug combination (Fitzgerald et al., 2006; Yeh et al., 2006; Breitinger, 2012). Rather, Loewe and Bliss models should be used as data exploratory approaches, with a major purpose to identify potential synergistic drug combinations that warrant further mechanistic investigation, but not the other way around, i.e., using the mechanistic evidence to determine which reference model is more appropriate.

After reaching the common understanding on the model assumptions, the Saariselkä agreement allowed the researchers for the flexibility to choose a preferred reference model to evaluate interactions of multiple agents, with the only precondition that the names of the specific models need to be explicitly reported. Namely, depending on which model is used, a combination response greater or less than yEXP will be termed as Loewe synergy, Loewe antagonism, Bliss synergy or Bliss antagonism, respectively. Following these recommendations, the controversy over the definitions of synergy and antagonism seemed subsided. More recently, a variety of interaction assessment methods have been further developed and applied to a wide range of biological research fields. Notably, most of these methods can be traced back to the two basic model classes. For example, variants of the Loewe additive model include combination index (Lee et al., 2007; Chou, 2010), isobologram analysis (Tallarida, 2006) and response surface models (Greco et al., 1995; Kong and Lee, 2006); variants of the Bliss independence model include various synergy contour approaches (Fitzgerald et al., 2006; Zhao et al., 2014).

### What is Synergy?

Paradoxically, even with the clear distinction that has been made between the reference models, we feel that the fundamental question still remains unanswered, if not becoming even more serious: What is synergy after all? Since the expected combination responses yBLISS and yLOEWE most often are not identical (Berenbaum, 1989), choosing the model to use has become

```
Frontiers in Pharmacology | www.frontiersin.org September 2015 | Volume 6 | Article 181 |
```
a practical burden for a researcher who tries to draw solid biological conclusions out from the data. Due to the lack of practical guidelines, the model selection has become a personal preference or largely a convention that has been followed in a particular research field without clear reasons (Lee et al., 2007; Zhao et al., 2014). There has been a tendency to favor a model that yields a lower expected combination response, as it results in a higher likelihood of detecting stronger synergy. To make matters even worse, there has been often a dilemma when a drug combination is classified as synergistic according to one model but antagonistic according to the other (Cokol et al., 2011). The Saariselkä agreement, unfortunately, seem to have failed to provide any recommendations for solving these practical issues. What the Saariselkä agreement achieved was a compromise for accepting individualized claims, but the ultimate aim to advance the consensus knowledge on the degree of interaction has remained largely missing.

To ease the model selection burden, we propose here the use of new terminology that incorporates both of the two reference models, together with the single drug responses, to distinguish non-interaction, synergy and antagonism. With simple algebra, one can show that max(y1, y2) ≤ yBLISS. For the Loewe additivity model with a monotonically increasing dose-response relation, one can also show that max(y1, y2) ≤ yLOEWE. We note that, max(y1, y2) is also the expected response from a popular reference model, called highest single agent (HSA) model (Berenbaum, 1989). If the combination response y<sup>c</sup> is lower than max(y1, y2), then one would intuitively infer antagonism. Therefore, we may use max(y1, y2), to distinguish antagonism from non-interaction. Similarly, one can use the response of the less effective single drug, that is min(y1, y2), to further distinguish between weak and strong antagonisms. For distinguishing synergy from non-interaction the answer is less obvious, as it depends on the comparison between yBLISS and yLOEWE. There has been considerable interest in the mathematical relationships between the Bliss independence and the Loewe additivity models to understand how much difference in the characterization of drug interaction one can expect when choosing one model over another (see e.g., Goldoni and Johansson, 2007). In particular, two authors of the Saariselkä agreement have reported results from such comparisons (Drescher and Boedeker, 1995; Dressler et al., 1999). They showed that yLOEWE > yBLISS is generally observed for very steep dose-response curves, while yLOEWE < yBLISS when the curves become more flat. Since, yBLISS and yLOEWE differ in a complex way depending on the parameterization of the dose-response functions, we propose two cut-offs, min(yBLISS, yLOEWE) and max(yBLISS, yLOEWE), for characterizing synergistic combinations. We reason that the consistency between the Bliss independence and the Loewe additivity models should be indicative of the degree of synergy: If both the Bliss model and the Loewe model classify a drug combination as synergistic, that is, y<sup>c</sup> > max(yBLISS, yLOEWE), then we call it a strong synergy; If the combination is classified as synergistic according to one model only, that is, min(yBLISS, yLOEWE) < y<sup>c</sup> < max(yBLISS, yLOEWE), then it is called weak synergy. Finally, non-interacting drugs have max(y1, y2) < y<sup>c</sup> < min(yBLISS, yLOEWE), reflecting our view that non-interaction should also be defined similarly as a range, rather than a single point as in the individual reference models. Given such a classification, one may continue to develop statistical testing methods for evaluation of its significance for replicate data. To facilitate better understanding of these definitions, we designed an interaction barometer that enables a systematic comparison of these proposed interaction terms along an axis of drug combination response y<sup>c</sup> (**Figure 1**).

The benefits of adopting the proposed terminology for the degree of interaction are two-fold. First, the definitions of synergy and antagonism are based on a simultaneous evaluation of the two reference models, as well as the individual drug responses. Such a data-driven approach avoids any pre-defined preference either for the Bliss independence or Loewe additivity when characterizing drug interactions, and thus it minimizes the biases toward either of the models. This is consistent with the idea that any synergy model should be treated as an exploratory ranking statistic for prioritization of the most potent combinations for further evaluation, rather than a "true model" for explaining synergy or antagonism mechanisms. Further, this terminology enables a more intuitive definition of non-interaction, under which the combination response may be higher than the single drug response, but not as high as the expected responses from the Bliss independence and the Loewe additivity models. Note that a drug combination falling into such an interval would be classified as antagonistic according to both of the two models, but since it produces a higher response than the single drugs, one would rather characterize it as an additive effect or non-interaction. For the sake of clarity, we would call it non-interaction, and in fact discourage the terms additive or additivity since these may be confused with the additivity implicated by the Loewe additivity model. The interval of non-interaction shown in **Figure 1** is positioned at the center of the barometer as a gray zone for those drug combinations with no clear evidence in support of either synergy or antagonism.

Secondly, the different interaction terms are positioned along the common response axis (e.g., measured as the percentage inhibitions of cell growth), which makes it easier to relate the degree of drug interaction with its outcome in the drug response. With the proposed interaction barometer (**Figure 1**), one can immediately tell the differences between the drug combination response y<sup>c</sup> and the responses of individual drugs (y1, y2), as well as the expected combination responses (yBLISS, yLOEWE) based on the two reference models. The clear correspondence between the degree of synergy and the combination response is in many ways superior to the use of an interaction index, such as combination index or other similar approaches (see e.g., Lee et al., 2007; Lee, 2010), which tend to be less obvious to interpret in terms of response boosting. For example, a combination index of 0.1 has been considered as a very strong synergy by Chou (2006), but how much extra response the synergy can produce for the drug combination is difficult to tell. In contrast, with the interaction barometer one can easily visualize the levels of boosted response of the combination compared to the single drugs or reference models. From the model development perspective,

the graphical representation of the competing reference models in the interaction barometer may facilitate a more systematic comparison among different approaches. For example, when new reference models are introduced, one can always position its expected combination response onto the barometer to enable better understanding of its relationships with the existing models.

### Synergy vs. Efficacy

So far, we have merely discussed about the assessment of drug interactions using the difference between the observed combination response and its theoretical expectation, i.e., y<sup>c</sup> − yEXP, to classify a drug combination as synergistic or antagonistic. However, a drug combination can be also classified as either effective or ineffective based solely on its actual combination response y<sup>c</sup> . It is important not to confuse these two concepts, synergy and efficacy, as the nomenclatures are related but should not be treated the same. Synergy is a measure of the degree of drug interaction, while efficacy is a measure of phenotypic response of a drug combination. It is possible that a drug combination is highly synergistic, while its actual response may be insufficient to reach therapeutic efficacy. On the other hand, a drug combination that exhibits strong response does not necessarily imply a synergistic interaction. For instance, only one of its component drugs may produce the response, while the other one is simply lowering the adverse effect of the first drug without affecting its on-target activity. In preclinical testing, a drug combination with strong synergy and efficacy should be prioritized for further mechanistic investigation, with an additional requirement of tolerable toxicity profile (Fitzgerald et al., 2006). Accordingly, the dosages of a drug combination are also important factor for clinical feasibility and for maintaining acceptable side effects. For instance, the concept of therapeutic synergy compares the therapeutic windows of the single agents to that of their combinations, instead of using compound efficacies alone (Kashif et al., 2014). However, the main focus of this review was the definition of synergistic interaction, and we refer those readers interested in the therapeutic significance of synergy in drug discovery to previous reviews (Fitzgerald et al., 2006; Sucher, 2014).

### Conclusion

The definition of synergistic interaction is still under debate. After a careful investigation of the Bliss independence model and the Loewe additivity model, we argue that, without jeopardizing the validity of both models, a more consistent terminology for classifying synergy and antagonism can be made. By comparing the observed combination response with the expected combination responses from the two models, as well as the single drug responses, one can classify the drug combination into five categories including strong antagonism, weak antagonism, non-interaction, weak synergy, and strong synergy. We propose the use of the interaction barometer to visualize the degree of interaction on the common axis of drug response, which has been shown to facilitate the interpretation and comparison between different combinations. We view our efforts as a continuation to what the Saariselkä agreement started more than 20 years ago but has not yet concluded: a consensus on concepts and terminology for interaction assessment. We acknowledge that our proposal is not yet solving the practical issues for analyzing real data which typically contain combination responses tested at different dose ranges. How to maximize the benefits of the interaction barometer to summarize the interaction patterns of a drug combination would be a source of future research initiatives. We hope that such a classification scheme will raise more discussions about the standardization of the interaction assessment, toward finally reaching a consensus not only on the definition itself, but also on the other important issues, such as the experimental design of combination experiments, their quality control and the statistical evaluation of synergy and antagonism.

### Funding

This work was supported by the Academy of Finland (grants 272437, 269862, 279163, and 292611 for TA, 277293 for KW); Cancer Society of Finland (JT, TA, and KW).

### References

Berenbaum, M. C. (1989). What is synergy? Pharmacol. Rev. 41, 93–141.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2015 Tang, Wennerberg and Aittokallio. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Identifying problematic drugs based on the characteristics of their targets

Tiago J. S. Lopes 1, 2, <sup>3</sup> \*, Jason E. Shoemaker 1, 3, Yukiko Matsuoka1, 4 , Yoshihiro Kawaoka1, 2, 3 and Hiroaki Kitano1, 4, 5, 6, 7 \*

<sup>1</sup> Japan Science and Technology Agency ERATO Kawaoka Infection-Induced Host Responses Project, Minato-ku, Japan, <sup>2</sup> Department of Pathobiological Sciences, School of Veterinary Medicine, Influenza Research Institute, University of Wisconsin-Madison, Madison, WI, USA, <sup>3</sup> Division of Virology, Department of Microbiology and Immunology, Institute of Medical Science, University of Tokyo, Tokyo, Japan, <sup>4</sup> The Systems Biology Institute, Tokyo, Japan, <sup>5</sup> Sony Computer Science Laboratories, Inc., Tokyo, Japan, <sup>6</sup> Integrated Open Systems Unit, Okinawa Institute of Science and Technology, Okinawa, Japan, <sup>7</sup> Laboratory for Disease Systems Modeling, RIKEN Center for Integrative Medical Sciences, Yokohama, Japan

#### Edited by:

Thomas J. Anastasio, University of Illinois at Urbana-Champaign, USA

#### Reviewed by:

Giovanna Cenini, Universität Bonn, Germany Ghanshyam Upadhyay, City College of New York, USA

#### \*Correspondence:

Tiago J. S. Lopes, Influenza Research Institute, University of Wisconsin, 575 Science Dr., Madison, WI 53711, USA dasilvalopes@wisc.edu; Hiroaki Kitano, The Systems Biology Institute, Falcon Building 5F, 5-6-9 Shirokanedai, Tokyo 108-0071, Japan kitano@sbi.jp

#### Specialty section:

This article was submitted to Experimental Pharmacology and Drug Discovery, a section of the journal Frontiers in Pharmacology

> Received: 18 June 2015 Accepted: 17 August 2015 Published: 01 September 2015

#### Citation:

Lopes TJS, Shoemaker JE, Matsuoka Y, Kawaoka Y and Kitano H (2015) Identifying problematic drugs based on the characteristics of their targets. Front. Pharmacol. 6:186. doi: 10.3389/fphar.2015.00186 Identifying promising compounds during the early stages of drug development is a major challenge for both academia and the pharmaceutical industry. The difficulties are even more pronounced when we consider multi-target pharmacology, where the compounds often target more than one protein, or multiple compounds are used together. Here, we address this problem by using machine learning and network analysis to process sequence and interaction data from human proteins to identify promising compounds. We used this strategy to identify properties that make certain proteins more likely to cause harmful effects when targeted; such proteins usually have domains commonly found throughout the human proteome. Additionally, since currently marketed drugs hit multiple targets simultaneously, we combined the information from individual proteins to devise a score that quantifies the likelihood of a compound being harmful to humans. This approach enabled us to distinguish between approved and problematic drugs with an accuracy of 60–70%. Moreover, our approach can be applied as soon as candidate drugs are available, as demonstrated with predictions for more than 5000 experimental drugs. These resources are available at http://sourceforge.net/projects/psin/.

Keywords: multi-target drugs, drug safety, target validation, machine learning, protein networks, supervised learning

### Introduction

New compounds are traditionally discovered by using large biological screening techniques to identify substances that cause the desired effects. While this approach has been effective for years and produced the drugs used today, technological advances are shifting the drug discovery process toward a more rational approach, with computational drug-design and pathway analysis playing major roles. With the costs of compound design dramatically increasing and most of these funds being spent on drugs that never make it to market (Munos, 2009; Scannell et al., 2012), there is a clear need for new technologies to develop more specific, less toxic compounds.

Recently, in silico analyses have been successfully applied throughout the drug discovery pipeline. Examples include methods to help understand the changes caused by candidate compounds in protein interaction networks (Csermely et al., 2005; Yildirim et al., 2007), and algorithms to develop specific ligands that inhibit the activity of pathogen proteins (Fleishman et al., 2011; Whitehead et al., 2012). In addition, computational analyses in key studies have revealed the off-targets of drug candidates and predicted important side effects (Keiser et al., 2007, 2009; Campillos et al., 2008; Yamanishi et al., 2008; Liu et al., 2011; Lounkine et al., 2012). These studies have shown that computational analyses are an essential part of drug discovery. Yet the early identification of problematic drugs remains a major challenge.

Here, we propose a method to distinguish between compounds that are safe and those likely to be harmful. For this purpose, we considered the targets of more than 1800 approved and problematic drugs (i.e., withdrawn from market, or halted in development due to safety concerns). To study the properties of these targets, we created a protein similarity network (PSIN), in which the proteins are connected only if their sequences are similar. We found that the centrality measures of the PSIN network clearly indicated which human proteins are likely to cause harmful effects if their activities are modulated by drugs; our analysis suggested that ∼5000 human proteins had characteristics that resembled those of targets of problematic drugs. Next, by using machine learning techniques, we developed an index (called the Rejection Score) to quantify the likelihood of a candidate drug being problematic. Although some substances were difficult to classify and obtained intermediate scores, most were consistent with their status of being approved or problematic.

Finally, based on the targets of more than 5000 experimental substances from major databases (∼700 of which are currently undergoing pre-clinical or clinical evaluation), we predicted which drugs have a high likelihood of approval by regulatory agencies. This process of validation of individual proteins and assignment of rejection scores to candidate drugs should improve gene target selection and candidate prioritization in drug development.

### Materials and Methods

### Protein Sequence Comparison

To perform the protein sequence comparisons required to assemble the PSIN, we used a stand-alone version of BLAST toolkit v2.2 (Camacho et al., 2009), and the human protein sequence database obtained from Uniprot (released in August 2012). We removed all splicing variants from this database, leaving only the first protein isoforms. We used the PSI-BLAST algorithm and the BLOSUM distance matrix, with a gapopen cost value of 11, gap-extension cost of 1, and minimum expectation value of 1e-03. We set the E-value threshold for inclusion in the multipass model at 1e-05, and six PSI-BLAST iterations or less (if the results converged before; i.e., no further sequences could be discovered in the database using the profile as the input query).

### Databases

The drugs, their targets, and status (approved, withdrawn, illicit, experimental, etc), were obtained from Drugbank (Wishart et al., 2006), from the Therapeutic Target Database (Zhu et al., 2012), and from ChEMBL (Gaulton et al., 2012). We merged all three databases into one dataset, with drugs containing targets from all three databases. While the first two databases have information about the legal status of the compounds, ChEMBL has only an indication of "therapeutic" or "non-therapeutic," and whether a therapeutic drug contains a black-box warning. Therefore, to study the characteristics of approved and problematic drugs. we used only the drugs for which we had legal status information from DrugBank or TTD, and the therapeutic drugs from ChEMBL (**Supplementary Table S1**).

From the ChEMBL database, we considered only proteins targeted by a compound if they had an IC50-value < 30 nM. Additionally, several drugs did not have targets present in the PSIN or in the protein-protein interaction network PPI. Because these drugs targeted proteins that were isolated from the rest of the network, only drugs with at least one target present in the PSIN or PPI were considered, and for drugs with multiple targets, those targets not present in either the PSIN or the PPI were removed before any analysis was done.

For the protein-protein interaction analysis, we used HIPPIE (Human Protein-Protein Interaction Reference Schaefer et al., 2012).

### Assigning Drugs Approved or Problematic Labels

The legal status of the drugs from the two databases was highly heterogeneous, containing approved drugs (those available by prescription or over the counter), illicit, and withdrawn drugs (those removed from the market or that had their development halted due safety or efficacy concerns). For our study, we were interested in understanding what distinguishes drugs successfully used to treat patients from those causing drug attrition or those that were withdrawn from the market due harmful effects. Thus, we first identified drugs that were withdrawn from the market in several countries and classified them as problematic. Second, drugs that were discontinued during clinical trials due to safety or efficacy issues were also considered problematic. In contrast, we considered a drug "successful" if it was available for purchase. Additionally, after detailed inspection of drugs classified as "illicit" in the DrugBank database, we verified that they were mainly used to treat psychological disorders (e.g., anxiety, schizophrenia, and insomnia), and had the potential for abuse and addiction. Since these substances can be obtained in most countries when prescribed by clinicians, we concluded that most of them were not in fact illicit, but rather "controlled substances," in which case, we also considered them approved drugs (**Supplementary Table S2**).

### Centralities, Averages, and Relevance Measures

For each protein network, we calculated the betweenness of a node v as:

$$B(\nu) = \sum \frac{sij(\nu)}{sij}, \quad \text{with i } \ne j, \,\nu \ne i \text{ and } \nu \ne j$$

where sij is the number of shortest paths between the nodes i and j and sij (v) is the fraction of those shortest paths passing through node v.

Burt's Constraint was calculated as:

$$C(i) = \sum\_{j} (p\_{ij} + \sum\_{q} p\_{iq} p\_{qj})^2, \quad \text{with } q \neq i, j, \text{ and } j \neq i$$

where piqpqj is the product between the proportional strength of the node i's relationship with node q, and the proportional strength of the node q's relationship with node j. The details of these calculations in their original sociological context were reported by Burt (1992, 2004).

In addition, when considering multiple targets of the same drug, we transformed all individual PSIN network measures to the log10 scale and then combined their centrality measures by calculating their arithmetic means.

### Implementation, Data Analysis, and Pre-processing

The computations involving pre-processing and machine learning classifiers were performed by using the Weka suite for data mining (Frank et al., 2004); our code was written in Java and all algorithms were used with their default parameters. The Support Vector Machine was implemented using LibSVM, with the code available at https://weka.wikispaces.com/LibSVM (visited in August 2015).

The statistical and network analyses were performed in R. Additionally, we used the iGraph package (Csardi and Nepusz, 2006) for the network analysis, the poweRlaw package for powerlaw fits, and the ROCR (Sing et al., 2005) package to create the ROC curves.

Pre-processing involved four steps: (1) not all proteins had centrality values in the PPI or in the PSIN, hence, we filled those missing values with the mean of the training and test sets separately by using the Weka function ReplaceMissingValues; (2) we had to over-sample the smaller class because our dataset contained more instances from the approved class than from the problematic class. Hence, we used the SMOTE (Chawla et al., 2002) algorithm for this task, with an oversample proportion of 500% and 8 nearest neighbors. We used the Tomek links (Tomek, 1976) method to remove instances whose nearest neighbors belonged to the opposite class. This strategy proved very effective relative to other pre-processing alternatives to deal with unbalanced, overlapping datasets (Batista et al., 2004). (3) we removed the instances that were on the "border" of different classes, i.e., instances that were the nearest neighbors of several instances from different classes (see Figures 1, 2 of Batista et al., 2004); and (4) we ran a preliminary cleaning step by using Multilayer Perceptron exclusively on the training set to remove the misclassified instances (we used the RemoveMisclassified routine from the Weka package).

For the training and testing procedure, we removed the drugs that targeted the same set of proteins. For example, if a hypothetical drug in the test set targeted proteins A, B, and C, then all other compounds in the training set that targeted A, B, and C were discarded. Additionally, we removed ∼100 drugs that targeted the same proteins but had conflicting classifications (i.e., some were approved and others problematic). This ensured that we had no redundant instances in the dataset, and that the same targets were not simultaneously in the training and test sets. After calculating the mean centrality measures of all drug targets, they were scaled to the interval [0,1] by using the R package Reshape http://had.co.nz/reshape - visited in June 2015.

### Results

### Network Characteristics

A protein similarity network is distinct from a protein-protein interaction network (PPI) because in the former, neighbor proteins do not necessarily interact or regulate each other's activities; instead, two proteins are connected only if their amino acid sequences are similar. Although other protein networks exist (Weston et al., 2004; Camoglu et al., 2006; Zhang and Grigorov, 2006; Atkinson et al., 2009; Rattei et al., 2010; Valavanis et al., 2010), they suffer from shortcomings such as the use of small protein datasets, employing information other than amino acid sequences, not being specific for human proteins, or not using signature-based methods to assess protein similarity. These shortcomings led us to create a network with the characteristics required to study the properties of drug targets. We used PSI-BLAST (Altschul et al., 1997) to query and compare the ∼20,000 human protein sequences in the Uniprot database. BLAST searches were not reciprocal (i.e., searches with "protein A" identified "protein B" as similar, but searches with "protein B" did not necessarily identify "protein A" as similar). Therefore, to establish a link between two nodes (proteins) in the PSIN we considered only bidirectional hits.

The PSIN has ∼17,000 proteins connected by ∼1,700,000 edges. The network does not have a single large component; rather, it has more than 800 smaller connected components. We used the degree (the number of edges—i.e., neighbors each node has) to quantify the connectivity of the PSIN. Its nodes have an average of 200 connections, with the most connected having ∼2600 neighbors. We verified that, similar to PPI networks, the degree distribution of the PSIN also fits in the power-law distribution, where many nodes have a few connections and a few nodes have many connections (**Supplementary Figure 1**). Nodes with up to 500 neighbors were connected to other nodes of similar degree, and above this point, the nodes were usually connected to those with 400–500 connections (**Figure 1A**).

In the PSIN, the proteins are connected to each other by similarity of one or more shared domains. We verified that most of the low-degree nodes and their neighbors belonged to the same families, because these proteins had domains found only in a few other proteins. For instance, the family of alphadefensin proteins, responsible for Gram-negative antibacterial activity, formed a cluster of only five proteins connected to each other through their exclusive and characteristic defensin domain. In contrast, highly connected proteins had a mixture of common and rare domains. For example, the notch1 protein has several repetitions of its 7 domains, which are connected to 1445 proteins from more than 50 families (**Figure 1B**).

In summary, the degree distribution in the PSIN resembled the power-law, like other protein networks. Low-degree proteins comprised rare domains and were mainly connected to members of the same family; high-degree proteins were composed of both rare and common domains, and were connected to members of several different families.

### Network Characteristics of Drug Targets

By using the PSIN and a PPI database, we searched for characteristics that discriminated between the targets of approved and problematic drugs.

We obtained all drugs and their reported targets from three major drug-target databases (see Materials and Methods), and merged all three databases taking into account the different drug names and synonyms; the merged dataset contained 1802 drugs for which their legal status was available (approved, withdraw, illicit, etc., **Supplementary Table S1**), and more than 5000 experimental drugs (**Supplementary Table S2**). Next, we assigned a simplified class label to each of the 1802 drugs, indicating whether the drug was considered safe and was marketed (approved), or if it had had its development halted or was withdrawn from the market (problematic).

We observed that the targets of approved and problematic drugs largely overlapped (**Figure 2A**), and there were more reported targets in the combined databases for the approved drugs than for the problematic drugs (**Figure 2B**). This is due to the strict requirements for drug approval by regulatory agencies, since before going to market, companies must provide detailed reports about modes of action, and after a compound is released, researchers from academia often report additional targets.

For each target of the approved and problematic drugs, in addition to the degree (the number of neighbors a protein has in the network), we calculated their betweenness, closeness centrality, and Burt's constraint in the PSIN and PPI networks. The betweenness describes how central a node is, counting the number of shortest paths that must pass through that node to connect the other nodes in the network. The closeness centrality from a node measures how many steps are necessary to reach every other node. The Burt's constraint, was first employed in a socio-psychological context, where the author studied the location of individuals in a large social network and quantified which individuals are in a position of advantage, located between groups, and have access to information and resources from different environments (Burt, 1992, 2004) (**Figure 2C**).

Compared with proteins targeted by the approved drugs, those targeted by problematic compounds had a significantly higher degree in both networks, and much lower closeness centrality and Burt's constraint values (**Figure 3**; for each centrality measure, One-Way ANOVA, p < 0.0001, followed by Tukey's HSD test; **Supplementary Figure 2**). In contrast, we observed no significant differences in the betweenness values of proteins targeted by problematic and those targeted by approved drugs in the PSIN or in the PPI network. These findings indicate that while targets of approved drugs have protein domains that are not shared among many other proteins and are involved in fewer interactions, targets of problematic drugs have domains that are more common throughout the proteome and have more protein interactions reported.

**Figure 3** suggests that proteins targeted by approved compounds and problematic compounds have characteristics similar to problematic targets. This finding could shed light on why some drugs are approved, while others are rejected, even though they target the same protein. For example, VEGF receptors, which are involved in blood vessel growth, are popular anti-cancer targets. Usually, the drugs targeting these receptors

FIGURE 2 | (A) Although most targets of approved drugs are exclusive, the problematic targets are almost entirely covered by the approved category. Between parentheses are the number of singleton proteins in the PSIN. (B) Approved and problematic drugs have different numbers of reported targets. While most problematic drugs have only one target reported, approved drugs have several—identified either by the community after the drug is marketed or by companies as part of the drug-approval process. (C) The Burt's constraint was proposed in a sociological context to study positions of advantage for individuals in a group. In this simple example, if the nodes are individuals, on the left no node can negotiate or bargain with the others, since they all have alternative connections. However, on the right, if a structural hole exists, Node 1 is in a better position, since the other two nodes may not be aware of each other's existence;hence, Node 1 is less "constrained" than the other two. In a protein similarity context, proteins with low constraint values are generally those with several common domains, located between different protein families. In contrast, proteins with large constraint values are the peripheral nodes, with a few domains shared among only a few other proteins.

FIGURE 3 | (A–D) In general, targets of problematic drugs have high degrees and closeness centralities in the PSIN and PPI networks. However, their betweenness values are not significantly different from the targets of approved drugs in either protein network (One-Way ANOVA, \*\*\*p << 0.0001 and \*p > 0.05, sample sizes for each group are the same as depicted in Figure 2A). The closeness from the targets of both networks was close to two main values, differing by only decimal digits; therefore, we rounded the values to their closest integer, namely 17 or 19 in the PSIN and 14 or 18 in the PPI. While three PSIN centrality measures were found to be strong indicators of the differences between targets of problematic and approved drugs, the centrality measures of the PPI network could also detect these differences, albeit in a moderate fashion (Tukey's Honest Significance Difference—Supplementary Figure 2). Overall, this likely stems from the fact that the current PPIs still have only ∼10,000 proteins and numerous false-positive interactions; with new proteins and high-quality interactions being constantly added, we expect this to change in the future.

cause major side effects (Roodhart et al., 2008) including hypertension, coagulation disorders, and neurotoxicity. We found that the three VEGF receptors had characteristics similar to other targets of problematic drugs, suggesting that compounds that target proteins with these characteristics will either be unapproved or, if approved, will likely cause harmful effects.

Next, we asked how many proteins share the characteristics of the approved and problematic targets. We observed that ∼65% of problematic drug-targets had a degree value > 110 and a Burt's constraint <0.025, whereas only ∼22% of approved drugtargets had these characteristics. In the PPI, ∼62% of problematic targets had a degree value >11 and a Burt's constraint <0.1, compared with ∼33% of approved drug-targets with these

characteristics. In general, these values are ones that show the highest separation between the targets of problematic and approved drugs (**Supplementary Figure 3**). When we consider the centrality measures separately, in the PSIN, ∼7600 proteins had degree and Burt's constraint values similar to those of problematic targets, and in the PPI, ∼4600 had degree and constraint values similar to other proteins of the problematic group. When we consider the characteristics of both networks together, ∼1200 proteins had measures that closely resembled those of problematic drug targets. Notably, this group has several popular targets of anticancer drugs, including almost all members of the cyclin-dependent kinase family, aurora kinases, and Pim/PLK serine threonine kinases. We then attempted to verify whether compounds targeted neighbor proteins in the PSIN. We built a contingency table counting the number of compounds that targeted neighbor proteins and the number of compounds that did not target neighbor proteins in the PSIN and tested its statistical significance by using Fisher's exact test. We obtained a p > 0.05, confirming that compounds often target proteins with no detectable sequence similarity (Keiser et al., 2007; Apsel et al., 2008; Yamanishi et al., 2008).

Taken together, our results indicate that the centrality measures calculated from the PSIN (and to a lesser extent from the PPI) can be used to distinguish between individual targets of approved and problematic drugs. These characteristics define the "danger zone" for therapeutic modulation, that is, they serve as indictors that modulating the activity of these proteins may be harmful to human health.

### Classifying Multi-target Drugs

Since drug targets of approved and problematic drugs could be distinguished individually, we asked whether evaluating the characteristics of all proteins targeted by a compound would help to predict the compound's safety and consequently, its approval or rejection. For classification and prediction tasks, supervised learning algorithms use a training set with examples assigned to different classes, and after a training phase, these algorithms attempt to predict the classes of instances they have not seen before. In our case, our training set comprised drugs, their targets, and their status (approved or problematic).

We built a dataset by calculating the centrality measures used above for each target of each drug; however, since most drugs have multiple targets, we combined the centrality measures of these targets using the means of their individual measures (see Materials and Methods). Overall, our dataset comprised 1802 drugs: 1445 approved and 357 problematic (**Supplementary Table S1**). As in most real-life scenarios, this dataset is characterized by the imbalance between the number of approved vs. problematic drugs–a characteristic that is notably difficult for machine learning algorithms (Batista et al., 2004). Therefore, we pre-processed the dataset to increase the sensitivity of the classifiers to the characteristics of the problematic class (see Materials and Methods), and for the classification routine, we compared the performance of 14 machine learning classifiers, namely KStar (Cleary and Trigg, 1995), Naive Bayes (John and Langley, 1995), J48 (Quinlan, 1993), Thresold Selector (Witten et al., 2011), Multilayer Perceptron (MLP)(Bishop, 1995), JRip (Cohen, 1995), IB1 (Aha et al., 1991), PART (Frank and Witten, 1998), END (Dong et al., 2005), Random Tree (Breiman, 2001), Rotation Forest (Rodríguez and Kuncheva, 2006), Random Forest (Breiman, 2001), Decorate (Melville and Mooney, 2003), and Support Vector Machines (SVM)(Cortes and Vapnik, 1995).

We asked how the classifiers perform if we use only the centrality measures from the PSIN, from the PPI, or a combination of both. We divided the input dataset into 70% of instances for training and 30% for prediction, with no overlapping drugs between them. Drugs that bind the same set of protein targets were removed to prevent obvious redundancies during classifier evaluation (see Materials and Methods). This procedure was repeated 100 times to quantify the prediction accuracy for each set of centrality measures. We then verified that although the topological characteristics of the PPI network could moderately distinguish between individual approved and problematic drugs, the predictive power of the machine learning algorithms was highest when using only the centrality measures of the PSIN network (**Supplementary Figure 4**), therefore, for subsequent analyses we used only the PSIN.

After close inspection of the drug-targets dataset, we realized that confounding factors of the drug-binding protein data might affect classifier performance (e.g., differing numbers of binding partners per drug, missing drug-targets in the protein networks). Therefore, we designed three tests to determine whether the PSIN data could enhance classifier performance over the performance obtained when only these confounding factors were considered. In the first test, we shuffled the class-labels (i.e., approved and problematic labels and the targets of all drugs—always keeping the same proportions as the original datasets—and compared them to the performances obtained when using the standard dataset, by using the 70%–30% division for training and testing sets, respectively. We verified that while most classifiers had an area under the ROC curve (AUC) close to 0.7 when trained to the complete dataset, the classifiers generally performed close to random guessing (AUC∼0.5) when trained using random datasets (p < 0.01, comparing the AUCs of each classifier, Wilcoxon two-sided signed-rank test).

In the second test, we developed a more stringent procedure wherein we created 40 randomized datasets also by shuffling the labels of the proteins in the PSIN, but here, for each single cross-validation, we first randomly split the standard dataset (i.e., that derived from the PSIN) into the 70%–30% training–test sets and determined the AUC. Next, we took each of the 40 randomized datasets and carefully divided them into training and testing data, making sure that we selected the same drugs that were used for training and testing with the standard data, and calculated the AUC. The 40 AUCs from the randomized data represented our null distribution, that is, the expected AUC achieved when drugs can bind any proteins. From this null distribution, we determined the likelihood that the AUC achieved by the standard data happened by chance, by comparing all runs of the randomized data to the standard dataset, and we verified that the overall classification procedure had better performance than random datasets in more than 88% of all comparisons (**Supplementary Figure 5**).

In the final test, we shuffled only the PSIN protein labels, removing the class distinctions discussed in **Figure 3**, while keeping the same network topology distributions (i.e., the power-law degree distribution). Again, we divided the dataset into 70%–30% for training and testing, respectively, repeating this procedure 100 times, and compared the results to the standard dataset. We observed that the classifiers performed considerably better than all of the shuffled networks (**Supplementary Figure 6**).

Together, these results demonstrate that the PSIN and the machine learning classifiers can overcome the effects of confounding factors, and distinguish multi-target problematic and approved drugs based solely on the network characteristics of their targets. Some of the algorithms outputted only binary classifications (END, Random Tree, SVM), or had the same underlying base classifier (J48); therefore, for further classifications, we used three algorithms (KStar, MLP, and Rotation forest) that were built using different underlying principles and outputted a probability that a drug belonged to the problematic class. This approach should compensate for any inevitable biases that all algorithms have.

### Predicting Drug Safety

After analyzing the capabilities of the classifiers, we used them as a prediction tool for new multi-target drugs. For fairness, all previous tests to study the characteristics of the classifiers had been performed without parameter or dataset optimization. However, to use as a prediction tool, it is desirable to fine-tuned the input dataset and include only the most meaningful examples in the training set.

While the approved drugs are generally compounds approved by regulatory agencies and successfully commercialized, the problematic drugs were deemed problematic for one or more of 10 different reasons (**Supplementary Table S3**). We wanted to test whether individually removing each of these 10 reasons from the input set would improve the classification.

For this purpose, we created 10 different datasets containing all of the approved and problematic drugs except those in each of the 10 groups that led to the drug failure. Each of these 10 datasets was then randomly divided into training and test sets by using the 70%–30% proportion, and the AUC of the classification was calculated. This procedure was repeated 100 times and our results showed that removal of one group of problematic drugs (those deemed "withdrawn") considerably improved the classification (p < 0.001, One-Way ANOVA, followed by Tukey's HSD test). Most likely, these compounds targeted proteins that were also targeted by approved drugs, but—in contrast to approved drugs—had harmful effects. Thus, by removing these drugs from the dataset, we removed a confounding factor and consequently reduced the false positives and increased the overall drug classification.

Next, from the complete dataset (with 1802 drugs), we selected one drug at a time and used the drugs from the optimized dataset as a training set. Importantly, we also removed from the training set any drug that targeted the same proteins as the drug being evaluated (the predictions are available in **Supplementary Table S4**).

We found that the classification of existing and experimental drugs into two classes (i.e., approved and problematic) could be over-simplistic; for drug development, it is more informative to quantify how likely a compound is to cause harm (Evans et al., 2009). Therefore, we created an index, which we named the "Rejection Score" (RS), by using the average of the probabilities calculated by use of the three chosen classifiers. We used this index to indicate whether a compound was predicted to be safe (RS close to 0.0) or more likely to be toxic (RS close to 1.0).

We found that 55% of approved drugs had a RS < 0.02 (**Figure 4A**, **Table 1**); yet, only 23% of problematic drugs had RSs close to 0.02. Conversely, 23% of approved drugs and 61% of problematic compounds had RSs greater than 0.9. Beyond this point, we observe a sharp increase in the number of problematic drugs and a slow increase in the number of approved drugs; hence, drugs with RSs of 0.9–0.95 could be considered "highrisk" compounds that are likely to cause strong side effects (nonetheless, their use may be warranted to treat life-threatening diseases). It is important to note that these cut-offs are arbitrary and may vary depending on the risk that is deemed acceptable. Moreover, given that drugs have distinct numbers of targets, we tested the rejection score and observed that it moderately correlated with its number of targets (**Supplementary Figure 7**).

Next, we asked whether the Rejection Scores reflected the known adverse reactions of marketed drugs. To answer this question, we analyzed the labels and package inserts of 245 drugs obtained from the SIDER database (Kuhn et al., 2010). We chose drugs with RSs covering the full range of predicted scores (i.e., from 0.0 to 1.0); for all selected drugs, we listed the indicated precautions, contraindications, adverse reactions, warnings, and where available, boxed warnings (**Supplementary Table S5**).

Several drugs with an RS of 0.9–1.0 had associated warnings and cautions regarding the risk of severe reactions, for example, beta-blocker drugs (Pirbuterol, Atenolol, Alvimopan), which can cause life-threatening reactions including heart failure, bradycardia, and angina. The side effects of beta-blockers have been known for decades (Frishman, 1988), and the high RSs of these drugs likely stem from the fact that their target proteins (Beta-1 and Beta-2 adrenergic receptors) have characteristics similar to those of other problematic drug targets, namely a high-degree (500+), and a low Burt's constraint (<0.006).

At the other end of the scale, we found that drugs with RSs of 0.0–0.2 and of 0.2–0.4 have been commercially available for years; an example is Bumetanide (RS 0.26), a diuretic that may cause profound water loss and electrolyte depletion only if used in excess. Bumetanide's reported targets include proteins involved in the transport of potassium, chloride, and sodium, and these transporters have low PSIN degree values (∼24) and high Burt's constraints (∼0.14), which are characteristics of proteins that are relatively isolated in the network and that are similar to only a few others. Another example is Diazepam (RS 0.0005), a benzodiazepine broadly prescribed since 1963 to treat anxiety and insomnia that can cause unpleasant, but manageable side effects, such as nausea, skin rashes, and headache (Riss et al., 2008). The targets of Diazepam, gamma-aminobutyric acid (GABA) receptors, have low PSIN degree values (46), and a high Burt's constraint (0.0826), and form a cluster in which these

proteins connect to each other and to a few other GABA receptor subtypes.

compounds. Each chart contains the number of drugs of the respective group.

In addition to drugs classified as approved or problematic, we also calculated the RS of ∼5000 experimental drugs. Experimental drugs listed in TTD and Drugbank are those currently being evaluated in clinical and pre-clinical trials, that have already undergone such evaluation, did not make it to


TABLE 1 | Percentage of drugs classified according to their Rejection

a In this dataset, 357 drugs were deemed "Problematic"; 1445 were deemed "Approved." See Supplementary Table S4 for details.

market, or are being tested for alternative indications. We also predicted the RSs of non-therapeutic substances from ChEMBL. The respective compounds selected from all three databases are listed in **Supplementary Table S3**, and their RS predictions are shown in **Supplementary Table S6**.

For experimental drugs, while more than 50% of compounds were predicted to have high rejection scores (>0.95), only 15– 20% of compounds in clinical trials had scores similar to those of approved drugs (**Figure 4B**), suggesting that only a small number of these candidate compounds may be approved and not cause severe side effects. These observations are consistent with the known high attrition rates observed in the pharmaceutical industry (Kola and Landis, 2004).

To summarize, most of our predictions seem to match the status of marketed drugs and their reported adverse reactions, and in general, our results show that drugs with high RSs were more likely to be discontinued in their development or be withdrawn from the market after commercialization.

### Discussion

Here, we investigated the characteristics of proteins targeted by approved and problematic drugs. We found that they have distinguishing characteristics that can be readily identified by using protein similarity and protein-protein interaction networks. In addition, we used machine learning methods to devise a score to quantify the risks of a drug being harmful to human health. We used this approach to predict the safety of several drugs and found that, for the most part, the prediction is consistent with their status of approved or problematic.

Given the prediction accuracy of ∼70%, it is unlikely that pharmaceutical companies would use such prediction results for Go/No-Go decision on compounds that may be already in the late pre-clinical or clinical stages. However, such information may be highly valuable in creating an overall portfolio of the candidate compounds from their very large chemical libraries. If a pharmaceutical company uses our method to evaluate millions of compounds and decides to moderately bias their choices according to the ratings presented here, the long-term outcome may be very different. This will be increasingly important when assessing the suitability of multiple-target drugs—i.e., it is essential to efficiently identify combinations of targets and compounds that are both safe and effective.

Interestingly, while attempting to find general properties of individual proteins that may be associated with side effects, we found characteristics from protein networks that clearly distinguish between targets of approved and problematic drugs (**Figure 3**). The degree, betweenness, and closeness centrality of proteins are well known measures and are broadly used in network studies. The Burt's constraint is already used in sociological studies and interestingly, revealed itself as a strong indicator of protein druggability (**Figure 3B**). Together, these network characteristics define a new axis along which drug targets can be assessed for their viability. In addition, together with other considerations (e.g., existence of accessible binding pockets, and the location and time-point of expression), this methodology can help validate safe targets whose modulation is therapeutically relevant (Bunnage, 2011). Further, the correlation between the RSs and the severity of the reported adverse reactions was not perfect; some drugs with a high RS had only mild adverse reactions, while others with low RSs had clear warnings about their potential harm. A potential explanation is that some drugs targeting highly connected proteins were approved despite their known adverse reactions (e.g., anti-cancer drugs), while relatively safe drugs may not have been developed due to business decisions. Thus, the Rejection Score should not be considered in isolation, but in conjunction with the network centralities of the drug's targets. For instance, a drug with low RS targeting high-degree proteins suggests that this compound is likely to cause moderate or severe side effects but still has the potential to be approved, depending on its indication. Conversely, if a compound has a high RS and targets high-degree proteins, this is a strong indication that this compound will be problematic because only a few or no other drugs targeting proteins with similar characteristics have been successfully commercialized.

When developing therapeutic compounds it is difficult to assess beforehand which proteins can be targeted without causing major side effects but still overcome cell tolerance to perturbations (Kitano, 2007; Hopkins, 2008). Moreover, understanding and anticipating side effects can be complex. For instance, the weight gain, diabetes, and cardiovascular problems experienced by patients being treated with antipsychotic drugs have causes that remain unclear (De Hert et al., 2012), although mounting evidence suggests that a compound's side effects are caused by modulation of its primary and secondary targets (Xie et al., 2009; Correll et al., 2011; Lounkine et al., 2012).

Some cases were difficult to classify. A notable example was Thalidomide (RS = 0.0004), a drug used to treat morning sickness and removed from the market after it was associated with birth defects (Stephens et al., 2000). Due to its inhibition of blood vessel growth (D'Amato et al., 1994), this compound was investigated to treat cancer (Verheul et al., 1999) and appeared to improve the survival of multiple myeloma patients (Singhal et al., 1999). With a low RS, this drug targeted proteins with characteristics of both approved and problematic targets, indicating that in addition to network characteristics, it is essential to verify that the stage of development and the location for target inhibition is appropriate (i.e., inhibiting angiogenesis in tumors is desirable, but during limb development is catastrophic).

Naturally, our approach has limitations. First, there are falsepositives in both the PPI and PSIN networks. In the former, false-positives exist mainly due different experimental setups (Venkatesan et al., 2009), and in the latter, the algorithm used to build the network has a statistical cut-off, which, although strict, does not completely prevent proteins from being set as neighbors by chance. Second, the list of targets reported for each drug is incomplete and it varies from reporting several targets for approved drugs, to only a few for problematic compounds. As recently demonstrated, known drugs have multiple protein off-targets, bound with enough specificity to interfere with their functions (Campillos et al., 2008; Yamanishi et al., 2008; Keiser et al., 2009; Lounkine et al., 2012). Additionally, the targets reported in the database might show not only direct interaction, but also indirect activation or repression. Therefore, future studies should take into account the nature of the interaction between the compound and the members of pathways, as well as which isoforms of the proteins are the actual targets of the drugs (here, as in most pathway databases and protein networks, we considered only the first splicing variant). Another limitation of the approach presented here is the mismatch between the RSs of drugs and their reported adverse effects. The reason for this mismatch is known, but its resolution is markedly complicated. Drugs approved despite their known severe adverse effects, as well as drugs deemed problematic for business-related reasons, are confounding factors even for the most sophisticated classifiers. Ideally, one would separate drugs by indication or class, train classifiers using only the drugs in each group, and finally check the predictions against drug labels and the complete documentation for the drugs in the training set. However, this is not practical at present because no database has enough positive and negative examples from each drug class, and only a fraction of all results and adverse reactions observed during clinical trials are available at present (Wieseler et al., 2013). Therefore, to benefit from the methods presented here, it is important to always consider the RS of drug candidates together with the centrality of the protein targets. Nevertheless, the quality of protein networks and of the drug-target databases is constantly improving, and despite these limitations, the methods presented here predicted the safety or danger of 60–70% of known drugs (**Figure 4A**). Moreover, we believe that as more comprehensive post-translational and structural information becomes available, its integration to the PSIN will enhance its predictive capabilities, furthering our understanding of the mechanisms of drug action and its effects on protein structures.

Finally, our prediction that several experimental drugs may not be approved provides evidence that the approaches presented here can be used long before a drug reaches the end of the development pipeline (Scannell et al., 2012); in fact, as soon as the targets of a compound are determined, we would recommend that that compound be subjected to the procedures described here.

### Author Contributions

TL designed the study, implemented the algorithms, analyzed the data and wrote the manuscript with feedback from all of the authors; JS helped elaborate the tests, implemented part of the algorithms and analyzed the data; YM helped conceive the study and analyzed the data; HK and YK provided the funding and coordinated the research efforts.

### Funding

This work was supported by the Japanese Science and Technology Agency, through the ERATO Kawaoka influenza host-response project.

### Acknowledgments

We thank Susan Watson for editing the manuscript, and Takeshi Hase, Ronaldo Prati, and Miguel Andrade-Navarro for constructive comments and fruitful discussions.

### Supplementary Material

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fphar. 2015.00186

Supplementary Figure 1 | We verified that the PSIN degree distribution could be fitted in the power-law. The figure depicts the degree distribution (black circles) and the curve-fit (in red), with the modifying exponent parameter indicated in the figure. The fit is satisfactory for most nodes, but those with a higher degree have some deviation from the fitted curve. Additionally, we can observe the degree distribution for the HIPPIE database, and for a sequence of 15,000 random numbers chosen with replacement from the interval [1, 2000].

Supplementary Figure 2 | Tukey's Honest Significant Difference test indicating the differences between the groups of drug targets. APP stands for Approved drug-targets, PRO Problematic drug-targets, and BOTH are the proteins targeted by drugs from both groups. The sample sizes of each group are the same as those shown in Figure 2A.

Supplementary Figure 3 | The figure shows the distribution of the degree and the Burt's constraint of the targets of approved and problematic drugs. In the PSIN, the highest separation between targets of approved and problematic drugs occurs at a degree value of 110; 75–80% of approved drug targets have values lower than this value, whereas only 30–35% of problematic drug targets have degree values smaller than 110. Additionally, we can observe that the separation of targets of problematic and approved drugs is much

narrower in the PPI than in the PSIN, making it more difficult to distinguish the potential therapeutic targets based on the centrality measures of the PPI network.

Supplementary Figure 4 | Boxplots with popular measures used to evaluate classifiers and different datasets. We can verify that for most measures and classifiers, the PSIN alone yields the best results.

Supplementary Figure 5 | We created randomized datasets that maintained the dependencies between the samples by shuffling the labels of the proteins in the PSIN and recalculating the average network topologies for each drug to create a single randomized dataset. In total, we created 40 randomized datasets. (A) Depicts an example, for a single cross-validation where we randomly split the standard dataset (i.e., that derived from the PSIN) into a training set and a test set and determined the AUC. Next, we took each of the 40 randomized datasets, carefully split them into training and testing data making sure we selected the same drugs that were used for training and testing with the standard data, and calculated the AUC. The 40 AUCs from the randomized data represent our null distribution, that is, the expected AUC achieved when drugs can bind any proteins. From this null distribution, we determined the likelihood that the AUC achieved by the standard data happened by chance. Thus, the p-value for a single iteration is the fraction of AUCs from the randomized data that was greater than or equal to the AUC achieved with the standard data; for example, a p-value of 0.1 would indicate that in 10% of cases the randomized data yielded an AUC equivalent to the standard run (B), Dataset #0 is the standard run, the remaining are the randomized datasets. We then repeated this for 120 cross-validation iterations (totaling 4920 calculations) and calculated the average and the median p-values observed across all iterations of the three classifiers used to formulate the Rejection Score. (C) Depicts histograms of the p-values achieved from 120 cross-validation iterations.

Supplementary Figure 6 | Depicted are the Area Under the ROC Curves (AUCs) where the Complete Dataset (i.e., original PSIN), is compared with a dataset in which we kept the networks intact, but shuffled the protein labels. We observe that most classifiers obtained higher AUCs than the shuffled datasets, indicating that even though confounding factors (e.g., number of targets per drug) play a role in the classification, the classifiers were still able to provide significant added value to the overall procedure. In each panel, the groups were compared by using the Wilcoxon two-sided signed-rank test, and in all cases, they were distinct with p << 0.01.

Supplementary Figure 7 | We observed a negative correlation (−0.35) between the number of targets a drug has and its rejection score (if the number of targets was in the logarithmic scale, the correlation was −0.51), yet, most drugs with only a few or with many targets still had similar rejection scores. The figure also suggests that we were not capturing the characteristics of only one or two of the databases.

Supplementary Table 1 | Drugs, targets, and their regulatory status.

Supplementary Table 2 | Experimental drugs and their targets.

Supplementary Table 3 | Criteria used for the drug categorization.

Supplementary Table 4 | Computational classification of marketed drugs.

Supplementary Table 5 | Drugs, Rejection Scores, and reported adverse effects.

Supplementary Table 6 | Computational classification of experimental drugs.

### References


Apsel, B., Blair, J. A., Gonzalez, B., Nazif, T. M., Feldman, M. E., Aizenstein, B., et al. (2008). Targeted polypharmacology: discovery of dual inhibitors of tyrosine and phosphoinositide kinases. Nat. Chem. Biol. 4, 691–699. doi: 10.1038/nchembio.117


explain the side effects of cetp inhibitors. PLoS Comput. Biol. 5:e1000387. doi: 10.1371/journal.pcbi.1000387


**Conflict of Interest Statement:** Tiago J. S. Lopes, Yoshihiro Kawaoka, and Hiroaki Kitano have a patent application for the drug classification approach presented in this manuscript. Jason E. Shoemaker and Yukiko Matsuoka declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2015 Lopes, Shoemaker, Matsuoka, Kawaoka and Kitano. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.