# THE EMERGING DISCIPLINE OF QUANTITATIVE SYSTEMS PHARMACOLOGY

EDITED BY: Tarek A. Leil and Sergey Ermakov PUBLISHED IN: Frontiers in Pharmacology

#### *Frontiers Copyright Statement*

*© Copyright 2007-2015 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-88919-642-5 DOI 10.3389/978-2-88919-642-5

## 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

## **THE EMERGING DISCIPLINE OF QUANTITATIVE SYSTEMS PHARMACOLOGY**

## Topic Editors:

**Tarek A. Leil,** Bristol-Myers Squibb, USA **Sergey Ermakov,** Bristol-Myers Squibb, USA

The promise of quantitative systems pharmacology: converting human physiology and pharmacology into computer models that can be used to generate insight into disease progression and therapeutic intervention. Copyright, Lightspring/Shutterstock.com

In 2011, the National Institutes of Health (NIH), in collaboration with leaders from the pharmaceutical industry and the academic community, published a white paper describing the emerging discipline of Quantitative Systems Pharmacology (QSP), and recommended the establishment of NIH-supported interdisciplinary research and training programs for QSP. QSP is still in its infancy, but has tremendous potential to change the way we approach biomedical research. QSP is really the integration of two disciplines that have been increasingly useful in biomedical research; "Systems Biology" and "Quantitative Pharmacology". Systems Biology is the field of biomedical research that seeks to understand the relationships between genes and biologically active molecules to develop qualitative models of these systems; and Quantitative Pharmacology is the field of biomedical research that seeks to use computer aided modeling and simulation to increase our

understanding of the pharmacokinetics (PK) and pharmacodynamics (PD) of drugs, and to aid in the design of pre-clinical and clinical experiments. The purpose of QSP modeling is to develop quantitative computer models of biological systems and disease processes, and the effects of drug PK and PD on those systems. QSP models allow testing of numerous potential experiments "in-silico" to eliminate those associated with a low probability of success, avoiding the potential costs of evaluating all of those failed experiments in the real world. At the same time, QSP models allow us to develop our understanding of the interaction between drugs and biological systems in a more systematic and rigorous manner. As the need to be more cost-efficient in the use of research funding increases, biomedical researchers will be required to gain the maximum insight from each experiment that is conducted. This need is even more acute in the pharmaceutical industry, where there is tremendous competition to develop innovative therapies in a highly regulated environment, combined with very high research and development (R&D) costs for bringing new drugs to market (~\$1.3 billion/drug). Analogous modeling & simulation approaches have been successfully integrated into other disciplines to improve the fundamental understanding of the science and to improve the efficiency of R&D (e.g., physics, engineering, economics, etc.). The biomedical research community has been slow to integrate computer aided modeling & simulation for many reasons: including the perception that biology and pharmacology are "too complex" and "too variable" to be modeled with mathematical equations; a lack of adequate graduate training programs; and the lack of support from government agencies that fund biomedical research. However, there is an active community of researchers in the pharmaceutical industry, the academic community, and government agencies that develop QSP and quantitative systems biology models and apply them both to better characterize and predict drug pharmacology and disease processes; as well as to improve efficiency and productivity in pharmaceutical R&D.

**Citation:** Leil, T. A., Ermakov, S., eds. (2015). The Emerging Discipline of Quantitative Systems Pharmacology. Lausanne: Frontiers Media. doi: 10.3389/978-2-88919-642-5

# Table of Contents


Athan Spiros, Patrick Roberts and Hugo Geerts


Sergey Ermakov, Peter Forster, Jyotsna Pagidala, Marko Miladinov, Albert Wang, Rebecca Baillie, Derek Bartlett, Mike Reed and Tarek A. Leil

## Editorial: The emerging discipline of quantitative systems pharmacology

Tarek A. Leil\* and Sergey Ermakov

*Bristol-Myers Squibb, Clinical Pharmacology and Pharmacometrics/Exploratory Clinical and Translational Research, Princeton, NJ, USA*

Keywords: systems pharmacology, systems biology, pharmacometrics, clinical pharmacology, modeling and simulation, in silico modeling, biomedical research

Quantitative Systems Pharmacology (QSP) has emerged recently as an approach that integrates knowledge coming from multiple disciplines including drug pharmacology, systems biology, physiology, mathematics and biochemistry. QSP was formally defined as a discipline and endorsed in the NIH White Paper (Sorger et al., 2011) in 2011. It has emerged at a time when the pharmaceutical industry is facing growing challenges in efficiency and productivity in R&D. QSP has the potential to help overcome some of these challenges. QSP models allow researchers to evaluate multiple hypotheses in-silico that would otherwise need to be evaluated experimentally. There is an expectation that the use of QSP will reduce the cost of R&D and the risks associated with uncertainties and gaps in our knowledge while bringing new therapies to patients.

QSP models are typically perceived as a research tool for hypothesis generation in drug discovery and exploratory clinical development; however, recently the US FDA used a QSP model to evaluate a proposed drug regimen for a new biologic therapy (Peterson and Riggs, 2015). In their communication with NPS Pharma, the FDA's Clinical Pharmacology division used a published QSP model of the calcium homeostasis system (Peterson and Riggs, 2010) to recommend an alternate dosing regimen for NATPARA, an injectable parathyroid hormone replacement drug used to control low blood calcium in patients with hypoparathyroidism. The use of a QSP model by the FDA to recommend an alternate dosing regimen to a sponsor highlights one of the important future applications of QSP models in regulatory interactions, and also represents an important milestone for the field. It is the first public instance of a QSP model being used by a regulatory agency to make a clinical recommendation to a sponsor. In the future, it is anticipated that it will be sponsors that leverage the utility of QSP models to support their own clinical decision making with regulatory agencies.

In the present research topic entitled, "The Emerging Discipline of Quantitative Systems Pharmacology," we provide an introduction to the developing field of QSP with a series of articles that describe models in different disease areas; showing how these models can be used to evaluate important research questions in pharmaceutical R&D. The research topic starts with a perspective article by Leil and Bertz (2014) that describes the history of how modeling tools where used in pharmacology and in drug development. The authors point out two major reasons for the growing importance of the QSP approach, (i) difficulty in finding new targets for therapies, and (ii) increasing cost and time required for developing a successful drug. Integrated into the drug development process, QSP modeling could become an effective tool to facilitate R&D; for example to translate knowledge between experimental systems (e.g., animal to human), and to predict the effects of multiple therapeutic interventions in combination; a task that would be inefficient using only clinical experimentation. The other articles in the research topic go on to demonstrate applications of QSP models to influence decision making in biomedical research.

One of the important applications of QSP modeling in pharmaceutical R&D is optimization of clinical dose and schedule. Oncology is one of the disease areas where the narrow therapeutic window of most therapeutics demands fine tuning of dose and schedule. Utilization of high doses

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

> \*Correspondence: *Tarek A. Leil, tarek.leil@bms.com*

#### Specialty section:

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

> Received: *29 May 2015* Accepted: *12 June 2015* Published: *30 June 2015*

#### Citation:

*Leil TA and Ermakov S (2015) Editorial: The emerging discipline of quantitative systems pharmacology. Front. Pharmacol. 6:129. doi: 10.3389/fphar.2015.00129* of anti-angiogenesis therapy can result in rapid suppression of angiogenesis and hypoxia leading to tumor shrinkage. Paradoxically, this can subsequently lead to reduced drug delivery to the tumor and resumption in tumor angiogenesis, followed by progression of tumor growth. The paper by Sharan and Woo (2015) discusses how to delay or prevent this from happening by optimizing the dose and regimen of the anti-angiogenesis therapies with a QSP model of angiogenesis.

Another important application of QSP models is to provide mechanistic explanations for clinical data that are often counterintuitive to the perceived mechanism of action (MoA) of a drug. Despite the fact that two SGLT2 inhibitors are already approved for use in patients, questions remain regarding their MoA and why efficacy is lower than expected based on their high potency and selectivity for SGLT2. The papers by Demin et al. (2014) and Lu et al. (2014) explored this issue using mechanistic models of renal tubular filtration and transport, incorporating the PK and MoA of SGLT2 inhibitors. Lu et al. (2014) proposed two possible explanations for the lower than expected efficacy; the residual activity of SGLT2 following inhibition in the renal tubules, and the compensatory effect of SGLT1. Demin et al. (2014) supported this hypothesis, but also offered an alternative in which the sites of action of SGLT2 inhibitors are located not in the lumen of the kidney's proximal tubules where the concentration of SGLT2 inhibitor is high, but perhaps in the proximal tubule where the concentration of inhibitor is lower. Complex dose-response dependencies are often encountered in many disease areas, for instance, in treatments of schizophrenia as investigated by Spiros et al (Spiros et al., 2014). With the use of a sophisticated QSP model of cognitive impairment in schizophrenia, the authors predicted an inverse U-shape doseresponse with glycine that is a consequence of the shifting balance between excitation and inhibition in the cortical network.

The application of mechanistic models for prediction of target dependent or independent toxicity in secondary tissues has been the focus of QSP for many years, as this is the most common reason for termination of the development of otherwise efficacious therapies. In order to predict toxicity using a mechanistic model, it is useful to incorporate a physiologically based PK (PBPK) model to predict drug concentrations in the target organ. Woodhead et al. (2014) described the use of a PBPK/toxicity model of drug induced liver injury (DILI) to evaluate the impact of bile salt export pump (BSEP) inhibition on hepatotoxicity in rats and humans. The DILI model was used to predict the responses to BSEP inhibitors with and without clinical hepatotoxicity. In accordance with the observed

## References


clinical data, the model predicted that bosentan, but not telmisartan, will cause mild hepatocellular ATP decline and serum ALT elevation. Similar to the research by Woodhead et al. (2014), Chetty et al. also relied on PBPK to predict drug concentrations in tissue, linking these concentrations to targetmediated pharmacodynamic (PD) effects (Chetty et al., 2014). They did so using the Simcyp PBPK simulator, a tool that has a built-in PBPK model and allows users to input drug specific parameters that have been measured in-vitro to predict in-vivo plasma and tissue PK. Simcyp has become an important tool for pharmaceutical R&D and for communicating with regulatory agencies regarding the PK of investigational drugs and their potential PK-related drug-interactions. Chetty et al. described the use of Simcyp to predict the tissue concentrations of four different drugs, metoprolol, nifedipine, triazolam, and zolpidem (Chetty et al., 2014). They demonstrated how polymorphisms in drug metabolizing enzymes would have an effect on the concentration of drugs in the target tissue and the subsequent impact on pharmacodynamics.

One of the technical issues that potentially limit more widespread use of QSP in biomedical research is the lack of an accepted standard modeling tool to facilitate sharing of models between researchers. The tool should permit evaluation of experimental scenarios of interest in a flexible computational environment for conducting efficient high throughput simulation. A potential solution was implemented in the web-based virtual systems pharmacology (ViSP) platform described in the article by Ermakov et al. (2014). The salient feature of ViSP is the use of a model in the form of an executable file. Matched with a full set of model parameters this executable becomes independent of the model structure and the software that were used to develop the model while preserving flexibility in the input parameters. These characteristics could be useful in the future when the utilization and sharing of QSP models becomes more widespread.

In conclusion, we would like to emphasize the emerging potential that QSP holds for biomedical research and in particular to improve decision making in pharmaceutical R&D. In order for this to occur, QSP will need to present more examples of the value that it can bring to the process of hypothesis generation and testing; examples like those included in this research topic.

## Acknowledgments

We would like to thank the contributing authors to this research topic for providing interesting and relevant research articles.


SGLT2 in renal glucose reabsorption in humans. Front. Pharmacol. 5:274. doi: 10.3389/fphar.2014.00274


schizophrenia: exploring glycine modulation of excitation-inhibition balance. Front. Pharmacol. 5:229. doi: 10.3389/fphar.2014. 00229

Woodhead, J. L., Yang, K., Siler, S. Q., Watkins, P. B., Brouwer, K. L. R., Barton, H. A., et al. (2014). Exploring BSEP inhibition-mediated toxicity with a mechanistic model of drug-induced liver injury. Front. Pharmacol. 5:240. doi: 10.3389/fphar.2014.00240

**Conflict of Interest Statement:** Both of the authors are employees of Bristol-Myers Squibb. 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 Leil and Ermakov. 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.

## Quantitative Systems Pharmacology can reduce attrition and improve productivity in pharmaceutical research and development

## **Tarek A. Leil \* and Richard Bertz**

Clinical Pharmacology and Pharmacometrics, Exploratory Clinical and Translational Research, Bristol-Myers Squibb, Princeton, NJ, USA

#### **Edited by:**

Chiranjib Chakraborty, Galgotias University, India

#### **Reviewed by:**

Filippo Caraci, University of Catania, Italy Robert A. McArthur, McArthur and Associates GmbH, Switzerland

#### **\*Correspondence:**

Tarek A. Leil, Clinical Pharmacology and Pharmacometrics, Exploratory Clinical and Translational Research, Bristol-Myers Squibb, P.O. Box 4000, Princeton, NJ 08543-4000, USA e-mail: tarek.leil@bms.com

The empirical hypothesis generation and testing approach to pharmaceutical research and development (R&D), and biomedical research has proven very effective over the last half-century; resulting in tremendous increases productivity and the rates of approval for new drug applications at the Food and Drug Administration (FDA). However, as discovery of new therapeutic approaches for diseases with unmet medical need becomes more challenging, the productivity and efficiency of the traditional approach to drug discovery and development is diminishing. Innovative approaches are needed, such as those offered by Quantitative Systems Pharmacology (QSP) modeling and simulation. This "systems" approach to modeling and simulation can be used to guide the hypothesis generation and testing process in pharmaceutical R&D, in a manner similar to its adoption in other industries in the past. Embedding QSP into the existing processes within pharmaceutical discovery and development will be required in order to realize the full beneficial impact of this innovative approach.

**Keywords: Quantitative Systems Pharmacology, pharmaceutical R&D, modeling and simulation, Systems Biology, pharmacometrics, drug discovery and development**

## **INTRODUCTION**

One of the oldest documents recording the process of drug discovery is the Ancient Egyptian "materia medica" dating to the sixteenth century B.C. This approach to drug discovery, based on empirical evidence from the natural world, is known as Pharmacognosy (de Pasquale, 1984), and was the primary means of drug discovery until the middle of the twentieth century. Advances in biochemistry, molecular and cellular biology, and medicinal chemistry during middle of the twentieth century resulted in a shift to a more hypothesis driven, mechanism-based approach to drug discovery (Drews, 2000). This approach includes mining of data from human epidemiology studies, combined with nonclinical *in vitro* and *in vivo* experiments to demonstrate the validity of a therapeutic target. The addition of high throughput chemical synthesis and screening permits identification of target selective, high affinity compounds. This hypothesis driven approach to pharmaceutical research and development (R&D) has been a tremendous advance relative to the ancient methods of Pharmacognosy, and resulted in a dramatic increase in the percentage of new drug applications (NDAs) approved by the Food and Drug Administration (FDA) since the early 1960s (**Figure 1A**; U.S. Food and Drug Administration, 2013b). However, there appears to be a decrease in the productivity (**Figure 1B**) of this approach that appears to be due to at least two major factors. One being the increasing difficulty in finding novel therapeutic targets, either for diseases with well established standards of care or those with unmet medical need (Pammolli et al., 2011; Scannell et al., 2012), and the other being the increasing cost associated with discovery and drug development of new drugs. The current average cost to bring a drug to market is \$1.5 billion, over 10 times higher than the cost in the 1970s (DiMasi and Grabowski, 2007; Scannell et al., 2012). The need for improved productivity in the pharmaceutical industry has been recognized by the FDA, with the establishment of its "Critical Path Initiative" in 2004. This initiative was intended to improve the drug and medical device development processes, the quality of evidence generated during development, as well as the outcomes of clinical use of these products (Woodcock and Woosley, 2008).

## **COMPUTER MODELING AND SIMULATION AS TOOLS TO IMPROVE PRODUCTIVITY**

A decline in productivity is to be expected for any industry as it matures, and the competition from established products increases. This decline in productivity is typically associated with increasing development costs, partly due to the difficulty in differentiating one product from another in the marketplace. Industries must find innovative ways to increase the probability of commercial success while at the same time decreasing development costs. Most industries eventually realize the value of computer aided modeling and simulation as one of the means for achieving both of these objectives. Computer aided modeling and simulation allows testing of numerous potential scenarios "*in silico*" to eliminate those associated with a low probability of success, avoiding the tremendous costs of evaluating all of those failed scenarios in the real world. Today, aerospace, automotive,

electronics, and other industries routinely incorporate modeling and simulation into their R&D processes (Woltosz, 2012). The pharmaceutical industry has been slow to integrate computer aided modeling and simulation for many reasons: including the perception that biology and pharmacology are "too complex" to be modeled with mathematical equations; a lack of adequate graduate training programs for pharmaceutical modeling and simulation scientists; and the lack of support from government funding agencies for academic research in computer aided modeling and simulation approaches for biomedical research. However, in the last decade, both the FDA and National Institutes of Health (NIH) have recognized the value of modeling and simulation in increasing productivity in biomedical research and pharmaceutical R&D. The FDA established its Pharmacometrics Division in 2009 to promote and evaluate the use of modeling and simulation approaches in regulatory submissions to the agency (U.S. Food and Drug Administration, 2013a). In 2011, the NIH published a white paper describing the emerging discipline of Quantitative Systems Pharmacology (QSP) modeling, and recommended the establishment of NIH-supported interdisciplinary research and training programs for QSP (Sorger et al., 2011). QSP modeling and simulation is a new term to describe the integration of two disciplines that have been increasingly useful in biomedical research and pharmaceutical R&D; "Systems Biology" and "Quantitative Pharmacology." Systems Biology is the field of biomedical research that seeks to characterize biological networks of interactions, including those between genes and biologically active molecules to develop models of these systems that are usually qualitative in nature. Quantitative Pharmacology (a.k.a. Pharmacometrics) is the field of biomedical research that seeks to use computer aided modeling and simulation to increase our understanding of the pharmacokinetics (PK) and pharmacodynamics (PD) of drugs, and to aid in the design of pre-clinical and clinical experiments. The purpose of QSP modeling is to develop quantitative computer models of biological systems and disease processes; as well as the effects of drug PK and PD on those systems.

## **PHARMACEUTICAL R&D AND QUANTITATIVE SYSTEMS PHARMACOLOGY**

Pharmaceutical R&D is a stepwise process where investment in further characterizing the pharmacology of a candidate molecule is incrementally increased as confidence in the molecule's probability of regulatory and commercial success increases. Investment initially starts with *in vitro* biochemical and pharmacology studies; then moves to animal pharmacology and toxicology studies, then to human healthy volunteer pharmacology and toxicology studies; and finally to large and expensive patient efficacy and safety studies. QSP models are based on the fundamental understanding of biological pathways, disease processes, and drug mechanisms of action. Therefore, they are very effective tools for integration of prior collected biological/pharmacological knowledge, formulation of pharmacological hypotheses, and for efficient translation between the various experimental models within pharmaceutical R&D. Key milestones in R&D where QSP models will be critical to increasing the probability of success will be in the target identification stage, the transition from pre-clinical to first in man studies, the transition from healthy volunteer to patient studies, and the transition from adult to pediatric (**Figure 2A**). These milestones represent the point where knowledge from one set of experimental models (e.g., animal) must be effectively translated to another set (e.g., human). QSP can facilitate this translation by formal integration of the knowledge from the original experimental model and generation of hypotheses for potential outcomes in the next experimental model. Computer aided simulations to guide the design of experiments intended to test those hypotheses.

In addition to their utility in translation between experimental models, QSP allows prediction of the effects of multiple therapeutic interventions in combination. As it becomes clear in many therapeutic areas that modulation of single drug targets is less effective (e.g., oncology, virology), the cost of testing all of the potential combinations in the clinic is prohibitive. QSP can provide a framework in which to evaluate these potential combinations prior to testing in the clinic, by providing a fundamental

#### **FIGURE 2 | Continued**

**(A)** Incorporation of Quantitative Systems Pharmacology modeling and simulation in pharmaceutical R&D. Drug discovery and development is a long and complex process with numerous transition periods where effective translation from one experimental model to the next is a challenge. Major transitions occur when moving to first in man studies and first in pediatric studies. The cycles of application of QSP modeling and simulation are defined by (1) integration of experimental data and biological knowledge to develop QSP models; (2) generation of hypotheses for potential outcomes in future experiments; (3) testing of

systems and quantitative understanding of how these different mechanisms will interact.

## **PHARMACEUTICAL R&D AND QUANTITATIVE SYSTEMS PHARMACOLOGY**

Published examples describing the use of QSP modeling and simulation to facilitate biomedical research and pharmaceutical R&D have been increasing in recent years. Most of these publications have been focused on PK, since the processes that govern drug absorption, distribution, metabolism, and excretion are better established compared to those that govern disease biology and PD (Edginton et al., 2008). Strougo et al. (2012) demonstrated the use of these physiologically based PK (PBPK) models for prediction of PK in children prior to the conduct of the first pediatric clinical studies. There are software packages that can be licensed with PBPK models incorporated that allow prediction of *in vivo* drug PK based on the *in vitro* properties of the molecule (Kuentz et al., 2006; Jamei et al., 2009). QSP models that predict both PK and PD are much more complex, and tend to be disease area specific. Vega-Villa et al. (2013) published a QSP model of the nitric oxide metabolic pathways and demonstrated the models ability to predict toxic methemoglobin levels in humans treated with nitric oxide. Geerts et al. (2013a) published a QSP model of cognitive deficit in schizophrenia and were able to simulate the enhancement of cognition with clozapine and risperidone, as well as the worsening of cognition with γ -aminobutyric acid (GABA) modulators lorazepam and flumazenil. There are software packages that can be licensed that allow prediction of both PK and PD for a variety of drugs and mechanisms of action, but at much greater expense compared with those used for PBPK alone (Shoda et al., 2010; Eissing et al., 2011). Agoram and Demin (2011) reported on the development and application of a QSP model of the PD of the 5-lipoxygenase (5-LO) pathway. This QSP model has been used to explain the complex PK–PD relationship of zileuton, a marketed 5-LO inhibitor (Karelina et al., 2010). The model was able to demonstrate the mechanism behind the longer duration, but similar magnitude, of action with the 600 vs. 400 mg dose of zileuton. More important for its utility in drug discovery, this QSP model could be used to predict the PK and PD of a new molecule or combination of molecules intended to modulate another component of the 5-LO pathway based solely on the biophysical properties of the molecule and its potency at the target. One could then design a series of *in vivo* experiments to validate the hypotheses generated from these predictions.

those hypotheses with experiments that have been designed via simulation from QSP models. **(B)** Integration of QSP models in pharmaceutical R&D process. The model development team should be a sub-team of existing drug discovery and development teams. The goals of the model development team are to develop QSP models for their particular disease area and/or to apply existing QSP models to facilitate key milestones in discovery and development. The model development process should be rigorous and stepwise, so that models that are developed can used broadly in the disease area and can be used to communicate with regulatory agencies.

Quantitative Systems Pharmacology holds great promise in being able to uncover innovative therapeutic paradigms for complex multi-factorial diseases such as Alzheimer's and diabetes. Because these diseases involve multiple physiological processes and can affect multiple organs, QSP can provide an integrated understanding of the pathology as well as the possible complex counter-intuitive results of therapeutic intervention.

## **INTEGRATION OF QUANTITATIVE SYSTEMS PHARMACOLOGY INTO PHARMACEUTICAL R&D**

In order to leverage QSP to accomplish this, it must be properly integrated into the decision making process in pharmaceutical R&D. As mentioned above, it is possible to license PBPK and QSP models to facilitate decision making in pharmaceutical discovery and development. However, since licensing is often limited to a few specialized functions, this approach decreases the flexibility and utility of QSP across the different functions within R&D. QSP models should be integrated into the processes of discovery and development within pharmaceutical companies in order to maximize their potential benefit on R&D efficiency and productivity. To better integrate these models into the existing processes, they can be developed internally within the drug discovery and development teams. In addition, the teams must be organized and educated to support this integration. **Figure 2B** shows the integration of the QSP model development team with the drug discovery and development team. Development of a QSP model for prediction of clinical PK and PD should be a rigorous and stepwise process, with three main steps:


In addition to the modeling engineer, who will incorporate the equations into the model based on the agreed upon scope, critical individuals on the model development team include the following:


Once a QSP model is developed, it should be regularly updated with relevant internal and external research. In addition, it should be made available to all scientists in R&D that work in that particular disease area. The QSP model can be made flexible such that it can be readily adapted to other species that may be of interest in drug discovery (e.g., rat, monkey, rabbit, etc.). This would facilitate translation of experiments between these species and human.

## **OUTLOOK FOR QUANTITATIVE SYSTEMS PHARMACOLOGY IN PHARMACEUTICAL R&D**

Quantitative Systems Pharmacology holds great promise in being able to uncover innovative therapeutic paradigms for complex multi-factorial diseases such as Alzheimer's and multiple sclerosis. Because these diseases involve multiple physiological processes and can affect multiple organs, QSP can provide an integrated understanding of the pathology as well as the possible complex results of therapeutic intervention. QSP thus offers pharmaceutical R&D an innovative way to conduct at drug discovery and development, particularly in diseases that are poorly translated from animal disease models. Geerts et al. (2013b) recently published an article on how QSP, when combined with phenotypic screening and preclinical animal models, could be used to address the bottleneck in both cognitive and neuropsychiatric drug discovery and development for Alzheimer's disease. For such complex diseases that are poorly translated from animal disease models, target-focused drug discovery holds little promise for finding successful therapies. However, pharmaceutical companies are large and bureaucratic, and dramatic changes to the direction in which R&D is conducted may be adopted rather slowly.

There are smaller pharmaceutical and biotech companies that are fully integrating QSP into their biomedical R&D processes. One example is Merrimack Pharmaceuticals in Cambridge, MA, USA; founded by MIT professor of Biology and Biological Engineering Michael Yaffe. Merrimack states on their website, "We are a Systems Biology company. We believe that improving cancer care requires a systems-based understanding of the dynamic interactions within a cancer cell and its environment" (Merrimack Pharmaceuticals Inc., 2014). The scientific and financial communities will be watching these small companies, and if their model for a systems-based approach to R&D is successful, the pressure will be increased for "large pharma" to more fully adopt such innovative approaches into their discovery and development processes.

## **REFERENCES**


ofMedicalProductsandTobacco/CDER/ucm167032.htm [accessed October 13, 2014].


**Conflict of Interest Statement:** The authors are employees of Bristol-Myers Squibb.

*Received: 15 August 2014; accepted: 24 October 2014; published online: 10 November 2014.*

*Citation: Leil TA and Bertz R (2014) Quantitative Systems Pharmacology can reduce attrition and improve productivity in pharmaceutical research and development. Front. Pharmacol. 5:247. doi: 10.3389/fphar.2014.00247*

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

*Copyright © 2014 Leil and Bertz. 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.*

## Systems pharmacology approaches for optimization of antiangiogenic therapies: challenges and opportunities

## **Satish Sharan and Sukyung Woo\***

Department of Pharmaceutical Sciences, College of Pharmacy, The University of Oklahoma Health Sciences Center, Oklahoma City, OK, USA

#### **Edited by:**

Tarek A. Leil, Bristol-Myers Squibb, USA

#### **Reviewed by:**

Xiaogang Wu, Institute for Systems Biology, USA Jeffrey S. Barrett, Sanofi Pharmaceuticals, USA

#### **\*Correspondence:**

Sukyung Woo, Department of Pharmaceutical Sciences, College of Pharmacy, The University of Oklahoma Health Sciences Center, CPB 331, 1110 N. Stonewall Avenue, Oklahoma City, OK 73117, USA e-mail: sukyung-woo@ouhsc.edu

Targeted therapies have become an important therapeutic paradigm for multiple malignancies. The rapid development of resistance to these therapies impedes the successful management of advanced cancer. Due to the redundancy in angiogenic signaling, alternative proangiogenic factors are activated upon treatment with anti-VEGF agents. Higher doses of the agents lead to greater stimulation of compensatory proangiogenic pathways that limit the therapeutic efficacy of VEGF-targeted drugs and produce escape mechanisms for tumor. Evidence suggests that dose intensity and schedules affect the dynamics of the development of this resistance. Thus, an optimal dosing regimen is crucial to maximizing the therapeutic benefit of antiangiogenic agents and limiting treatment resistance. A systems pharmacology approach using multiscale computational modeling can facilitate a mechanistic understanding of these dynamics of angiogenic biomarkers and their impacts on tumor reduction and resistance. Herein, we discuss a systems pharmacology approach integrating the biology of VEGF-targeted therapy resistance, including circulating biomarkers, and pharmacodynamics to enable the optimization of antiangiogenic therapy for therapeutic gains.

**Keywords: systems pharmacology, biomarkers, dose selection, antiangiogenic therapy, resistance, targeted therapies, biologically effective dose, bed**

## **INTRODUCTION**

Therapeutic intervention in diseases takes place within a milieu of factors, including drug pharmacokinetics, signaling pathways, mechanisms of drug action, and compensatory processes. Studying any single pathway, mechanism of action, or interactive process in isolation has limited value in improving our understanding of the complexity of disease physiology. This is evident in the redundancy of signaling networks, feedback, and cross talk between multiple regulatory processes (Moriya et al., 1996; Pawson and Warner, 2007; Logue and Morrison, 2012). A systems approach is required to quantitatively integrate underlying disease and information contributing to treatment response and resistance. Systems pharmacology combines large scale experimental studies, pharmacokinetics, mechanisms of action, signaling pathways, adaptation mechanisms, biomarker, and pharmacodynamic data in a quantitative framework utilizing computational methods. This approach can facilitate understanding of disease systems, their mechanisms of action and pathways, and hypothesis development (Agoram, 2014). Systems pharmacology further offers a tool for translational considerations from non-clinical models to patients, realizing the bench to bedside paradigm (Allerheiligen, 2010; Kreeger and Lauffenburger, 2010; van der Graaf and Benson, 2011; Demin et al., 2013; Rogers et al., 2013; Vicini and van der Graaf, 2013; Visser et al., 2014). Herein we discuss systems pharmacology approaches to achieve a mechanistic understanding of the dynamics of circulatory biomarkers for antiangiogenic agents, thereby guiding selection of doses that can maximize the therapeutic benefits.

## **CHALLENGES IN ANTIANGIOGENIC THERAPIES**

Angiogenesis is critical for tumor growth and metastasis. VEGF signaling is an extensively studied pathway for blocking tumor angiogenesis. Several antiangiogenic agents targeting the VEGFpathways have been approved and are important modalities in the management of advanced cancers. Bevacizumab, a therapeutic antibody targeting VEGF and various VEGF receptor tyrosine kinase inhibitors (TKIs), have shown clinical benefit in solid tumors. However, the benefits of VEGF-targeted agents are shortlived and resistance to anti-VEGF agents rapidly emerges after an initial response phase, leading to restored tumor growth and progression. This rapid development of resistance to therapy constitutes a major clinical obstacle to providing extended therapeutic benefits with this class of drugs. Thus, effective strategies are needed to delay or prevent resistance to VEGF antiangiogenics.

Resistance to antiangiogenic agents arises through multiple mechanisms, including the activation of compensatory responses that are mediated by malignant cells and stroma cells within the microenvironment. Angiogenesis is a highly adaptive biological process. Tumors can resume angiogenesis and progress using diverse angiogenic signaling, including VEGF, FGF, HGF, PDGF, PlGF, and several proangiogenic cytokines. Numerous compensatory angiogenic factors are upregulated upon anti-VEGF therapy in a dose-dependent manner (Ebos et al., 2007), suggesting that dose intensity and frequency influence the development of therapy resistance. Higher doses of anti-VEGF therapy can create favorable conditions for metastasis by upregulating these growth factors (Ebos et al., 2009a). This emphasizes the importance of finding optimal dosing schedules for anti-VEGF therapy. The current dosing approach does not consider the best way to delay or prevent resistance to VEGF-targeted therapy, and thereby improve patient survival beyond a few months (Jubb et al., 2006; Azad et al., 2008; Cannistra, 2008).

## **BIOLOGICALLY EFFECTIVE DOSES OF ANTI-VEGF THERAPY**

Oncology drug development often involves the maximumtolerated dose (MTD)-based paradigm, even when data suggest that a drug maximally inhibits its target at lower doses. The recent analysis by the U.S. Food and Drug Administration (FDA) showed that inappropriate dose selection was the major cause of postmarketing requirements for oncology drugs approved between 2011 and 2013 (Prowell, 2014). Clinically recommended doses are often derived based on their safety profiles. Toxicity has been the primary end point for conventional dose-finding strategies (Parulekar and Eisenhauer, 2004; Le Tourneau et al., 2009). Since antiangiogenic therapies are mostly cytostatic in nature, they do not always conform to the concept that MTD produces maximum benefits (Sleijfer and Wiemer, 2008). Studies have revealed better therapeutic benefits when lower doses of antiangiogenic therapies were used in combination with other treatments (Kabbinavar et al., 2003; Huang et al., 2012). Similar results were observed with other targeted therapies, such as mammalian target of rapamycin kinase inhibitor, in which a lower dose of 25 mg was selected as the recommended dose for treatment after testing 25, 75, and 250 mg doses (Atkins et al., 2004). Likewise, in a study of 24 consecutive Phase I clinical trials, in which 97.7% of participants received targeted agents, patients receiving lower (≤25% MTD) doses responded as well as those patients receiving medium (25– 75% MTD) or high (≥75% MTD) doses (Jain et al., 2010). These findings support the concept that higher doses are not necessarily the most effective.

Higher doses of anti-VEGF therapies can lead to pronounced anti-vascular effects and, subsequently, hypoxia in the tumor, e.g., treatment-induced hypoxia. Treatment-induced hypoxia stimulates several compensatory biological processes to circumvent continued VEGF inhibition, leading to resistance to therapy (Harris, 2002; Casanovas et al., 2005; Drevs et al., 2005; Kerbel, 2005; Mizukami et al., 2005; Hendriksen et al., 2009; Casanovas, 2011). This excessive pruning also leads to reduced delivery of therapies into the tumor (Jain, 2005; Van der Veldt et al., 2012; Van der Veldt and Lammertsma, 2014). Tumors have abnormal vasculature, which leads to an abnormal blood supply that produces hypoxic regions in the tumor. Hypoxia has been also implicated in tumor progression by increasing genomic instability (Nelson et al., 2004) and selection of more malignant cancer stem cells with increased metastatic potential (Bottaro and Liotta, 2003; Conley et al., 2012). Therefore, antiangiogenic therapy can produce more challenges than benefits, if it is inappropriately administered (Huang et al., 2013b; Jain, 2013, 2014). This is consistent with RK Jain's vascular normalization concept, in which the judicious use of antiangiogenic drugs can lead to more efficient delivery of drugs and oxygen to the tumor cells (Jain, 2005). Utilization of the vascular normalization strategy has been shown to improve cancer immunotherapy (Huang et al., 2012, 2013a) and survival in glioblastoma patients (Sorensen et al., 2012; Emblem et al., 2013). Therefore, there is a critical need to find the biologically effective dose (BED) that balances between normalization and excessive anti-vascular effects from antiangiogenic agents, since suboptimal and higher doses can fail to alleviate hypoxia. Further, the BED can minimize stimulation of alternative, compensatory proangiogenic signals in response to treatment-induced hypoxia, and thus limit the rapid development of treatment resistance, extending the therapeutic benefits of antiangiogenic agents (Jubb et al., 2006; Azad et al., 2008; Cannistra, 2008).

## **DYNAMICS OF CIRCULATING ANGIOGENIC BIOMARKERS**

The transient effects of antiangiogenic therapy predominantly result from a redundancy in the angiogenesis signaling that mediates tumor escape from anti-VEGF therapy. Many of the signaling molecules (circulating angiogenic factors, or CAF) within these compensatory pathways can be detected systemically in patients treated with VEGF-targeting agents. For example, increases in VEGF and PlGF, and decreases in VEGFR2 can be observed. These changes are considered a "class" effect of VEGFtargeted therapies (Jain et al., 2009). Many of these observed CAF changes are recapitulated in tumor-bearing mice in a dosedependent manner, and are correlated with antitumor activity (Ebos et al., 2007). Thus, CAF are increasingly recognized as important pharmacodynamic biomarkers for better understanding the treatment response and aiding in the identification of the optimal dosing schedules for VEGF-targeted therapy (Huang et al., 2013b). Understanding the molecular interactions between therapy-induced CAFs and resistance to VEGF-targeted agents can inform the development of strategies to delay or overcome resistance to antiangiogenic therapy (Jain et al., 2009; Clarke and Hurwitz, 2013).

These circulating biomarkers are dynamic, altered over the course of treatment by variables including *in vivo* drug concentrations (PK), changes in the tumors (e.g., antitumor effect and disease progression), the biological turnover of signaling molecules, compensatory mechanisms, tumor-independent CAF induction by normal cells in the host body, and the development of resistance. These diverse contributing factors create uncertainty when attempting to use dynamic biomarkers. Mathematical modeling can play an important role in understanding and utilizing the biomarkers to find the optimum biological dose and schedule which can delay the onset of therapy resistance (Duda et al., 2013). A recent study showcased the utility of computational models in identifying dosing schedules to manipulate the dynamics of the development of resistance to EGFR-targeted therapy (Foo et al., 2012; Dolgin, 2014). A systems pharmacology approach using multiscale computational modeling offers a tool to integrate the biology of response and resistance to VEGF-targeted therapy, including circulatory biomarkers and the pharmacokinetics/pharmacodynamics of antiangiogenic drugs (**Figure 1**), to optimize therapeutic gains.

## **SYSTEMS PHARMACOLOGY APPROACH TO ANTIANGIOGENIC THERAPY**

The major challenge in developing a systems pharmacology model is how to integrate the dynamics outside the cell (phar macokinetics) with their downstream effects in terms of protein

formation or pharmacodynamic effects. PK/PD modeling has been used to explain the relationship between pharmacokinetics and the end downstream effects. What is missing is the mechanistic information in between. Limiting our investigation to antiangiogenic therapy, we anticipate three major challenges to filling this gap: (1) determining the interaction of ligands to their receptors and perturbation by drug molecules, (2) integrating the ensuing signal from these receptors with the downstream protein production machinery, and (3) accounting for interaction between various cell types, that produces pharmacodynamics responses and resistance.

## **DRUG-TARGET INTERACTION**

Receptor occupancy theory is well-developed and can be readily utilized to integrate this process (Black and Leff, 1983; Black et al., 1985; Mager and Jusko, 2008; Chen et al., 2009; Goodman and Redberg, 2014). We must be mindful that biology is complex and there are many subtypes of ligands, receptors, and coreceptors which have varying degrees of affinity and modulatory functions. Ligand-receptor interaction for angiogenesis involves the VEGF family of ligands (VEGF-A, B, C, D, and PlGF), three main receptors (VEGFR-1, -2, and -3), co-receptors NRP-1 and NRP-2, and heparan sulfate proteoglycans. NRP- and, -2 and proteoglycans play modulatory roles in ligand-receptor interaction; even VEGF-A is alternatively spliced to form VEGFA121, VEGFA145, VEGFA165, and VEGFA189 (Hoeben et al., 2004; Koch et al., 2011; Tugues et al., 2011). Popel and colleagues have contributed extensively to our understanding of the kinetics and interaction of VEGF ligands and receptors (Stefanini et al., 2010; Finley et al., 2011, 2013; Finley and Popel, 2012, 2013; Tan et al., 2013). Although a potential contribution of other ligand and receptor isoforms and families may be recognized, several studies have simplified these interactions by accounting for the most important VEGF ligand, VEGF-A (165), and receptor VEGFR-2 interaction as the major players in angiogenesis (Sharan and Woo, 2014; Zhang et al., 2014). The combination of competitive ligand receptor binding and an inhibitory Hill function model can be used to explain the VEGF-induced VEGFR activation and inhibitor-induced VEGFR inactivation (Sharan and Woo, 2014).

## **SIGNAL TRANSDUCTION**

Signaling pathways are an important component of a systems pharmacology model, which links receptor-ligand interaction to pharmacodynamic outputs (Iyengar et al., 2012). VEGF binding to its receptors led to the phosphorylation of the tyrosine kinase domain, which in turn initiated the canonical downstream signaling cascades involved in proliferation, migration, survival, and permeability (Tugues et al., 2011). Ordinary differential equation (ODE)-based models, also termed mechanistic or physicochemical models (Birtwistle et al., 2013; Zhang et al., 2014), are often used to describe the canonical signaling cascades. The advantage of this approach being more mechanistic can help a personalized medicine paradigm by incorporating information related to genomic variation and mutation (Iyengar et al., 2012). The limitation of this approach is the currently incomplete mechanistic knowledge of several mediators, signaling processes, and parameter identifiability. Given the incomplete mechanistic knowledge of several mediators and signaling processes, or in the absence of measurement of mediator signaling molecules, a more empirical quantitative logic (QL; Kirouac et al., 2013; Kirouac and Onsum, 2013) or transduction model (Mager and Jusko, 2001) can be utilized to characterize signal transduction. The QL approach has been elegantly explained by Kirouac and Onsum (2013) in building multiscale models which capture the features of oncogenic signaling networks. The transduction model has the flexibility of handling multiscale events with different transit time parameters to account for time needed for signal transduction from receptor-ligand interaction to nucleus, time for cell machinery to form proteins, and to show the pharmacodynamic effects on tumor growth. It is vital to take a balanced approach between mechanistic representation of the signaling pathway and the model's predictive power (Sharan and Woo, 2014).

## **THERAPEUTIC AND COMPENSATORY RESPONSES TO ANTI-VEGF THERAPY**

Ideally, signaling events are linked to tumor growth kinetics. Tumor inhibitory effects of anti-VEGF agents can be described by adapting well-established models (Simeoni et al., 2004, 2013; Ribba et al., 2014). In addition to tumor growth inhibition, systemically quantifiable biomarkers such as CAFs can serve as an important measurement to identify disease progression, make dose selections, or stratify patients. Indirect response models can be effectively used to capture inhibition, stimulation, and turnover rates of biomarkers modulation (Mager et al., 2003). We can use non-linear feedback regulation to account for compensatory increases in circulatory biomarkers in response to treatment-induced hypoxia by anti-VEGF agents.

The contributions of host cells and stroma cells within the tumor microenvironment have been increasingly recognized to play an important role in cancer progression and treatment (Ebos et al., 2007, 2009b; Kerbel and Ebos, 2010; Jain, 2013; Stroh et al., 2014). This should be explored using a systems pharmacology model. Antiangiogenic therapies have been shown to upregulate various growth factors in healthy cells and are dosedependent in non-tumor-bearing mice (Ebos et al., 2007). This dose-dependency is also observed in healthy human volunteers (Lindauer et al., 2010). Thus, it is important to characterize tumor and host cell contributions to CAF modulation and to provide mechanistic information for interpreting biomarker data in respect to antiangiogenic treatments.

## **APPLICATION OF A SYSTEMS PHARMACOLOGY MODEL OF CAFs FOR DOSE OPTIMIZATION**

We have recently developed a systems pharmacology model that uses sunitinib as the test drug to quantify the link between *in vivo* drug concentrations (PK), target–drug interactions, the biological target pathway, antitumor activity, and compensatory signals leading to treatment resistance (**Figure 1**). We used the most frequently studied CAFs, including VEGF, PlGF, and sVEGFR2. Our model predictions were consistent with the time- and dosedependent changes in these hypoxia-derived CAFs following sunitinib given to mice at various dosages (Ebos et al., 2007). We then tested our model within a clinical setting to explain VEGF changes in patients with cancer who experience different treatment outcomes; we found that our predictions were consistent with the observed VEGF changes in patients receiving sunitinib for the treatment of metastatic renal cancer (Kontovinis et al., 2009). The stimulation/inhibition capacity and the hill coefficients of VEGF, PlGF, and sVEGFR2 in mice were similar to those reported in humans, indicating that system-specific parameters for conserved physiological processes such as angiogenesis are comparable across species (Sharan and Woo, 2014).

**FIGURE 2 | (A)** Relationships of therapeutic efficacy and modulation of VEGF and PlGF biomarkers to sunitinib doses. Percentage reduction in tumor volume (•) and fold change in PlGF (N) and VEGF () are shown at various doses of sunitinib at the end of study. At the dose of 40 mg/kg/day, ∼75% of tumor volume was reduced, with minimal upregulation of hypoxia-dependent CAF. Further dose escalation resulted in marginal therapeutic gain (<5%), but significant upregulation of CAF, which may indicate excessive anti-vascular effects. **(B)** Utilization of CAF biomarkers in the selection of biological dose of antiangiogenic drugs. The fold change in VEGF and PlGF may serve as a surrogate marker for excessive anti-vascular effects and, in turn, potential for emerging resistance. This illustrates how the biologically effective dose may be selected in a manner which does not invoke significant hypoxia and involves little stimulation of hypoxia-dependent CAF. Monitoring multiple CAFs will be advantageous, as each factor has a different dynamic range. PlGF has a wider dynamic range than VEGF, and results in higher fold change at the same dose. This provides an advantage over VEGF, because PlGF changes are more likely detectable even at lower doses.

Our model allows us to delineate CAF changes in the tumor microenvironment and host body during VEGF-targeted therapy and to assess their impacts on tumor response and resistance to therapy. This provided insight into the possible ways to utilize the CAF for better dose guidance of these therapies, either alone or in combination. We found a relationship of tumor reduction and compensatory increase in VEGF and PlGF with increasing sunitinib doses in xenograft mice (**Figure 2A**). The increase in proangiogenic factors was directly related to the dose, suggesting that these CAF can be used as biomarkers to determine the optimal dose for antiangiogenic drugs. For example, in an A431 xenograft mouse model, the maximum benefit from sunitinib treatment may be achieved at a dose of 20–40 mg/kg/day of sunitinib, as such doses produce no significant changes in VEGF or PlGF levels. Further dose escalation resulted in marginal therapeutic gain (<5%), but significant upregulation of hypoxia-dependent CAF, which may indicate excessive anti-vascular effects. As such, these CAF could be used to construct a therapeutic index for antiangiogenic agents. **Figure 2B** illustrates this CAF biomarker-based paradigm for dose selections in particular, balancing between antitumor effects and CAF changes. The CAF modulation may serve as a surrogate marker reflecting the anti-vascular effects of antiangiogenic treatment. Ligands of tyrosine kinase receptors have been found to confer resistance by engaging survival signals redundant to those of targeted kinase (Wilson et al., 2012). If we assume that higher changes in compensatory signals are associated with higher likelihood of early onset of resistance, antiangiogenic doses may be increased up to the level at which the CAF increase from their baseline is minimal (e.g., <2-fold for VEGF).

Increasing sunitinib doses also led to VEGF and PlGF stimulation with different magnitudes (**Figure 2B**). This differential stimulation of pro-angiogenic factors can be exploited to facilitate dose finding. Given the inter-individual variability and heterogeneity in tumor response, monitoring multiple biomarkers, rather than relying on a single marker, would be advantageous. We found that PlGF changes were ∼2-fold higher than VEGF changes at the same dose (Sharan and Woo, 2014). This finding suggests that PlGF has a wider dynamic range than VEGF, and can ensure better detection of its change even at lower doses. Thus, monitoring PlGF and VEGF can aid in ensuring that the therapy does not fall below the minimum effective dose nor go above the excessive anti-vasculature dose. This finding is consistent with the recent study in which increased PlGF, but not VEGF, was associated with patients responding to cediranib (Batchelor et al., 2013). While we illustrated the CAF-based dose-finding strategy using VEGF and PlGF, other CAF could be used, as different tumor types and drug targets can invoke different CAF dynamics. Many concepts and the mathematical framework are broadly applicable among several tumor types and different antiangiogenic agents, and could serve as a paradigm for determining the optimal dose of targeted therapies.

Antiangiogenic therapies are often administered in combination with chemotherapy. There is increasing interest in combining antiangiogenics with other targeted therapies in order to improve therapeutic outcomes. However, since the clinical doses of many targeted therapies are determined based on MTD rather than BED, we cannot easily deduce the dosage and schedules of combinations from single agent studies. When antiangiogenic therapies are combined with drugs of same class, excessive overlapping toxicities have resulted (Azad et al., 2008, 2009). In addition, antiangiogenic therapies at higher doses could reduce the efficacy of concomitant cytotoxic agents, most likely due to reduced drug delivery by excessive vessel pruning (Van der Veldt et al., 2012). The CAF biomarker-based approach could also be useful for determining the optimal dose of combination therapy. In combination therapy, the role of antiangiogenic drugs may be focused on vascular normalization. Other therapeutics can be targeted toward killing tumor cells. In such cases, it is desirable for antiangiogenic drugs to be administered at lower doses at which the stimulation of compensatory signaling is minimal.

## **CONCLUSION**

Therapy-induced CAF can be effectively utilized as pharmacodynamic biomarkers to find the optimal biological dose for antiangiogenic drugs. This will ensure that the therapy maintains minimum effective dose levels, without invoking much compensatory response from the system. Routine incorporation of biomarkers into future clinical trials will be critical for the optimization of anti-VEGF agents and development of next generation of antiangiogenic regimens. Biomarker studies can be augmented by imaging studies, such as dynamic contrast-enhanced MRI (DCE-MRI), or other imaging techniques to monitor vessel integrity, permeability of blood vessels, and tumor perfusion (Murukesh et al., 2010). Thus, future strategies will require circulating biomarkers and imaging with an integrated multi-scale computational tool to guide optimal dose selection for antiangiogenic agents. As more data become available in the future, with the advance of high throughput methods like genomic data, the main challenge will be vertically integrating those data. Systems pharmacology will offer a tool to vertically integrate knowledge from pharmacokinetics, mechanisms of action, genomics, biomarkers, toxicokinetics, and pharmacodynamics. This approach yields an informed perspective from which we can streamline drug discovery and development. The knowledge gained from this approach can provide an in-depth understanding and, hence, a better approach for achieving enduring therapeutic benefits from antiangiogenic therapy.

## **ACKNOWLEDGMENT**

Research reported in this publication was supported by the National Institutes of General Medical Sciences of the National Institutes of Health under Award Number P20GM103639. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

## **REFERENCES**


growth kinetics in xenograft models after administration of anticancer agents. *Cancer Res.* 64, 1094–1101. doi: 10.1158/0008-5472.CAN-03-2524


**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.

*Received: 05 September 2014; accepted: 06 February 2015; published online: 20 February 2015.*

*Citation: Sharan S and Woo S (2015) Systems pharmacology approaches for optimization of antiangiogenic therapies: challenges and opportunities. Front. Pharmacol. 6:33. doi: 10.3389/fphar.2015.00033*

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

*Copyright* © *2015 Sharan and Woo. 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.*

## Analysis of the efficacy of SGLT2 inhibitors using semi-mechanistic model

#### *Oleg Demin Jr. <sup>1</sup> \*, Tatiana Yakovleva1, Dmitry Kolobkov2 and Oleg Demin2*

*<sup>1</sup> Laboratory Alpha, Institute for Systems Biology Moscow, Moscow, Russia*

*<sup>2</sup> Institute for Systems Biology Moscow, Moscow, Russia*

### *Edited by:*

*Tarek A. Leil, Bristol-Myers Squibb, USA*

#### *Reviewed by:*

*Sukyung Woo, University of Oklahoma Health Science Center, USA Satyaprakash Nayak, Pfizer Inc., USA*

#### *\*Correspondence:*

*Oleg Demin Jr., Laboratory Alpha, Institute for Systems Biology Moscow, Nauchny proezd, 20, Bldg. 2, Technopark Slava, 117246 Moscow, Russia e-mail: demin\_ jr@insysbio.ru*

The Renal sodium-dependent glucose co-transporter 2 (SGLT2) is one of the most promising targets for the treatment of type 2 diabetes. Two SGLT2 inhibitors, dapagliflozin, and canagliflozin, have already been approved for use in USA and Europe; several additional compounds are also being developed for this purpose. Based on the *in vitro* IC50 values and plasma concentration of dapagliflozin measured in clinical trials, the marketed dosage of the drug was expected to almost completely inhibit SGLT2 function and reduce glucose reabsorption by 90%. However, the administration of dapagliflozin resulted in only 30–50% inhibition of reabsorption. This study was aimed at investigating the mechanism underlying the discrepancy between the expected and observed levels of glucose reabsorption. To this end, systems pharmacology models were developed to analyze the time profile of dapagliflozin, canagliflozin, ipragliflozin, empagliflozin, and tofogliflozin in the plasma and urine; their filtration and active secretion from the blood to the renal proximal tubules; reverse reabsorption; urinary excretion; and their inhibitory effect on SGLT2. The model shows that concentration levels of tofogliflozin, ipragliflozin, and empagliflozin are higher than levels of other inhibitors following administration of marketed SGLT2 inhibitors at labeled doses and non-marketed SGLT2 inhibitors at maximal doses (approved for phase 2/3 studies). All the compounds exhibited almost 100% inhibition of SGLT2. Based on the results of our model, two explanations for the observed low efficacy of SGLT2 inhibitors were supported: (1) the site of action of SGLT2 inhibitors is not in the lumen of the kidney's proximal tubules, but elsewhere (e.g., the kidneys proximal tubule cells); and (2) there are other transporters that could facilitate glucose reabsorption under the conditions of SGLT2 inhibition (e.g., other transporters of SGLT family).

**Keywords: SGLT-2, systems pharmacology modeling, Type 2 diabetes mellitus (T2DM), dapagliflozin**

## **INTRODUCTION**

Type 2 diabetes mellitus (T2DM) is a metabolic disorder that is characterized by hyperglycemia resulting from insulin resistance and relative lack of insulin. Current therapies for T2DM primarily address endocrine pathogenesis of insulin resistance and b-cell dysfunction. Consequently, many patients receive multiple glucose-lowering therapies and eventually require exogenous insulin administration. However, in clinical practice patients often fail to meet the targets for glycemic control (Ali et al., 2013), most frequently because of adverse effects–including hypoglycemia and weight gain–caused by therapeutic agents (Neumiller, 2014). Therefore, in order to avoid direct influence on insulin resistance and b-cell dysfunction, T2DM can be treated with inhibition of glucose reabsorption.

Glucose is reabsorbed to blood in the proximal tubules of the kidneys during the formation of primary urine. Where the reabsorption process is inhibited, glucose is excreted in urine and blood glucose concentration is reduced. Patients treated with inhibition of glucose reabsorption have a low risk of hypoglycemia because the mechanism of the treatment is independent of insulin release or endogenous glucose production (Misra, 2013). In the kidneys, families of glucose transporters and sodium-dependent glucose co-transporters are involved in glucose reabsorption. For example, sodium-dependent glucose cotransporter 2 (SGLT2) is thought to contribute to 90% of glucose reabsorption in kidneys (Liu et al., 2012), which makes it a promising target in T2DM therapy.

Currently, two SGLT2 inhibitors are marketed (dapagliflozin, canagliflozin) and several more are in development (e.g., empagliflozin, tofogliflozin). Based on inhibitory concentration 50% (IC50) values measured *in vitro* (Ohtake et al., 2012) and plasma concentrations measured in clinical trials (Yang et al., 2013), the marketed dosage of dapagliflozin was predicted to inhibit SGLT2 almost completely and thereby reduce glucose reabsorption by 90%. However, clinical trials demonstrated that dapagliflozin induced 50–80 g of urinary glucose excretion per day, which corresponded to only 30–50% inhibition of reabsorption (Komoroski et al., 2009a; Kasichayanula et al., 2011a). To explain these findings, several hypotheses were published (Liu et al., 2012). The following hypotheses explain the lower than predicted efficacy of SGLT2 inhibitors:


To explore these hypotheses further, it is possible to adopt a systems pharmacology modeling (SPM) approach.

Many mathematical models that describe development and treatment of T2DM appear in the literature, and Ajmera et al. (2013) present a detailed review of these models. Several mathematical models describing SGLT2 inhibitors are also presented. For example, three models describe the pharmacokinetics (PK) and pharmacodynamics (PD) of SGLT2 inhibitors in animals (Yamaguchi et al., 2011, 2012, 2013), two population PK models exist for empagliflozin (Riggs et al., 2013) and dapagliflozin (van der Walt et al., 2013), and Maurer et al. (2011) describes a PK/PD model for dapagliflozin in rats and humans. However, a model that describes the concentration of SGLT2 inhibitors at the potential site of action (i.e., the lumen of proximal tubule in the kidneys) is yet to be published. The level of a compound in plasma may differ significantly from that in the kidneys; therefore, prediction of the concentration of SGLT2 inhibitors in the lumen of the kidney's proximal tubules is important for understanding the PD effect of the drug.

The aim of this study was to construct a model that describes the active secretion of SGLT2 inhibitors from plasma into the lumen of the proximal tubules, reverse reabsorption, and urinary excretion. Using an SPM approach, our objective was to test the hypotheses used to explain the discrepancy between expected and observed levels of glucose reabsorption following administrations of SGLT2 inhibitors (see above). We also aimed to compare the efficacies of different SGLT2 inhibitors by simulating their concentration level in the lumen of the kidney's proximal tubules and estimating the level of inhibition produced during treatment in humans.

### **METHODS**

A family of semi-mechanistic PK/PD models was developed to describe administration, degradation, transport, glomerular filtration, active secretion, reverse reabsorption, and urinary excretion of 5 SGLT2 inhibitors (dapagliflozin, canagliflozin, ipragliflozin, empagliflozin, tofogliflozin). The assumptions used for model development are presented in **Table 1**. The models describe the PK of the drugs, inhibition of glucose reabsorption mediated by SGLT2, and levels of inhibition in transporters of the SGLT family. The models for dapagliflozin, ipragliflozin, and tofogliflozin, which have identical structures, include 4 compartments: plasma, peripheral compartment (tissues, organs), lumen of the kidney's proximal tubules, and urine. The models for canagliflozin and empagliflozin have the same structure and include 3 compartments: plasma, lumen of the kidney's proximal tubules, and urine. The rate equations for each model are similar, but many of the parameter values are specific to a particular drug. We chose to add the peripheral compartment, and a rate equation describing transport between plasma and peripheral compartment, to models describing activity of dapagliflozin, ipragliflozin, and tofogliflozin in order to achieve a better description of plasma PK data (see "Models verification strategy" Section in Supplementary Materials). The effects of compounds on SGLT2 and other transporters are described as functions of inhibitor concentration in the lumen of the kidney's proximal tubules. Schemes of the models are presented in **Figure 1**. Using the family of models, we were able to quantify, analyze, and compare PK and PD characteristics of the drugs. A system of ordinary differential equations, rate equations, and explicit functions are presented in **Table 2**. Model development and fitting procedures were performed using the DBSolve Optimum package (Gizzatkulov et al., 2010).

## **IDENTIFICATION OF MODEL PARAMETERS**

The strategy of model verification is presented in the Supplementary Materials. The dapagliflozin, ipragliflozin, and tofogliflozin models include 19 parameters: 4 physiological, 9 drug specific PK, and 6 drug specific PD. In the right hand side of canagliflozin or empagliflozin models, 17 parameters are included: 4 physiological, 7 drug specific PK, and 6 drug specific PD. The parameter values are specific to each SGLT2 inhibitor model, with the exception of 4 physiological parameters (*Vplasma*, *Vlumen*, *GFR*, and Q*urine*) that have equal values in each model. The physiological parameter values were taken (or calculated) from the literature, as were those for 2 PK parameters (*F* and *fup*) and 6 PD parameters (IC50 for SGLT1–6). The remaining 7 or 5 parameters were fitted against PK data on the dynamics of compounds in plasma and urine. The 95% confidence intervals

**Table 1 | Assumptions for models development.**



were calculated for fitted parameters. Parameter values, 95% confidence intervals, and the source of data for their identification are presented in Supplementary Table 1.

## **COMPARTMENT VOLUMES**

Volume of blood plasma was taken from the literature. The volume of the lumen of the kidney's proximal tubules was calculated using experimental data (see Supplementary Materials). In the models that described dapagliflozin, ipragliflozin, and tofogliflozin, distribution volume in peripheral compartments was fitted against PK data. Compartment volumes and the source of identification are presented in Supplementary Table 1.

## **SIMULATIONS**

Model simulations are presented as curves (simulated with optimal parameters values) with shadows (95% confidence bands) or as bars (simulated with optimal parameters values) with error bars (95% confidence bands). The description of simulating the 95% confidence bands is presented in the Supplementary Materials.

## **RESULTS**

All data available in the literature was used for model verification and validation (**Table 3**). Models were calibrated against clinical data obtained in trials where a single dose of SGLT2 inhibitor was administered. By applying the verification strategy described in the Supplementary Materials, we found that: (i) the active secretion of dapagliflozin, canagliflozin, and ipragliflozin from plasma to kidney lumen was equal to zero; (ii) the reabsorption of empagliflozin was equal to zero; and (iii) both active secretion and reabsorption of tofogliflozin were equal to zero. The appropriate choice of model parameters was validated against data obtained from single and multiple dose administration. **Figures 2**, **3** represent examples of dapagliflozin model calibration using plasma PK and urine recovery data respectively, where data was obtained from single dose trials. **Figures 4**, **5** also represent examples of dapagliflozin model validation against plasma PK and urine recovery data respectively, but data was obtained from multiple dose trials. The 95% confidence bands captured the variability in clinical data (**Figures 4**, **5**); therefore, we concluded that model validation was of satisfactory quality. Several additional figures demonstrating the quality of verification and validation in dapagliflozin, canagliflozin, ipragliflozin, empagliflozin, and tofogliflozin models are presented in the Supplementary Materials (Supplementary Figures 1–16). These figures demonstrate that the models enable a satisfactory description of plasma PK and recovery of SGLT2 inhibitors in urine.

## **CONCENTRATION OF SGLT2 INHIBITORS IN THE LUMEN OF THE KIDNEY'S PROXIMAL TUBULES**

To analyze and compare the concentration of drugs in the lumen of the kidney's proximal tubules–the predicted site of inhibition–we simulated time profiles of the drugs resulting



from oral administration at various dosages. **Figure 6** demonstrates that concentration levels of tofogliflozin, ipragliflozin, and empagliflozin are higher than levels of other inhibitors following administration of marketed SGLT2 inhibitors at labeled doses and non-marketed SGLT2 inhibitors at maximal doses (approved for phase 2/3 studies).

In order to investigate which of five compounds is most likely to accumulate at the potential site of action, we also simulated concentration levels of SGLT2 inhibitors in the lumen of the kidney's proximal tubules following a standardized dosage of the compounds equal to 20 mg. In this scenario, the model predicted that concentrations of empagliflozin and tofogliflozin were higher



than the levels of dapagliflozin, canagliflozin, and ipragliflozin (Supplementary Figure 17).

We also simulated and compared time profiles of SGLT2 inhibitor concentrations in plasma and lumen of the kidney. For all SGLT2 inhibitors, our model predicted that the concentration in the lumen of the kidneys was higher than in plasma (**Figures 7**, **8**, Supplementary Figures 18–20).

## **LEVELS OF SGLT2 INHIBITION DURING TREATMENT WITH SGLT2 INHIBITORS**

We simulated the average level of inhibition of glucose reabsorption mediated by SGLT2 and compared our findings with experimentally measured levels of glucose reabsorption inhibition following treatment with compounds (**Figures 9**, **10**, Supplementary Figures 21–23). In model simulations, the average level of inhibition of glucose reabsorption mediated by SGLT2 was higher than the

**FIGURE 3 | Example of verification of the dapagliflozin model using urine data.** Cumulative amount of dapagliflozin recovered in urine following simulation of a single administration of 50 mg. Curve represents model simulation and dots represent experimental data. Colors of dots correspond to different data sources: black—Obermeier et al. (2010); blue— Kasichayanula et al. (2011a); red—Kasichayanula et al. (2013).

**FIGURE 4 | Example of validation of the dapagliflozin model using plasma data.** Total level of dapagliflozin in plasma following simulation of a single administration, but at different doses (Yang et al., 2013). Curve represents model simulations and dots represent experimental data. Colors of curves and dots correspond to the different doses: black—5 mg, blue—10 mg. Model simulations are presented with 95% confidence bands.

level of total glucose reabsorption inhibition, and almost equal to 100%, following treatment with all experimental doses of SGLT2 inhibitors (**Figures 9**, **10**, Supplementary Figures 21–23).

To compare the efficacies of each compound, the level of SGLT2 inhibition was simulated after administration of equal doses (20 mg) (Supplementary Figure 24), labeled doses of marketed inhibitors, and maximal doses of other inhibitors approved for phase 2/3 studies (**Figure 11**). For all compounds, administration of equal or labeled/maximal approved doses resulted in almost 100% inhibition of SGLT2. Therefore, the difference in efficacy of the tested SGLT2 inhibitors was not significant.

## **EFFECT OF SGLT2 INHIBITORS ON OTHER TRANSPORTERS FROM THE SGLT FAMILY**

To predict the effect of dapagliflozin, canagliflozin, ipragliflozin, empagliflozin, and tofogliflozin on other members of the SGLT family, the inhibition level of these transporters in the lumen of the kidney's proximal tubules was simulated following administration of labeled doses of marketed inhibitors and maximal doses of other inhibitors approved for phase 2/3 studies. We found that the influence of the inhibitors on other transporters was substantial–approximately 50% or above for particular compounds (Supplementary Figure 25–29). The most effective inhibitor of SGLT1 was canagliflozin (Supplementary Figure 25), while ipragliflozin and tofogliflozin produced the strongest inhibition of SGLT3 (Supplementary Figure 26). While all other SGLT2 inhibitors have a similar effect on SGLT4, dapagliflozin is a weak inhibitor of this co-transporter (Supplementary Figure 27). Ipragliflozin and empagliflozin are the strongest SGLT5 inhibitors (Supplementary Figure 28), while canagliflozin has the strongest effect on SGLT6 (Supplementary Figure 29).

**FIGURE 6 | Concentration of SGLT2 inhibitors in the lumen of the kidney's proximal tubules.** Level of SGLT2 inhibitors in the lumen of the kidney's proximal tubules following simulation of multiple administrations of labeled doses of marketed SGLT2 inhibitors and maximal doses of other SGLT2 inhibitors approved for phase 2/3 studies. Colors of curves correspond to different compounds: black—10 mg QD dapagliflozin; blue—300 mg QD canagliflozin; red—25 mg QD empagliflozin; green—300 mg QD ipragliflozin; pink—40 mg QD tofogliflozin. Model simulations are presented with 95% confidence bands.

## **DISCUSSION**

SGLT2 is a promising target for the treatment of T2DM because its inhibition could lower glucose levels without directly influencing insulin resistance or b-cell function. However, in a

compartments following simulation of a single administration of 40 mg. Colors of curves correspond to various compartments. Black—total plasma concentration; blue—unbound plasma concentration; red—concentration in lumen. Model simulations are presented with 95% confidence bands.

previous study, administration of SGLT2 inhibitors resulted in only 30–50% inhibition of glucose reabsorption in human subjects even at high doses, which was in contrast to the expected inhibition of 90% based on *in vitro* potency and plasma levels (Liu et al., 2012). In our study, we applied a systems pharmacology modeling approach to test several hypotheses that attempt to explain the discrepancy between the expected and observed levels of glucose reabsorption following administration of SGLT2 inhibitors. We developed models, simulated concentrations of SGLT2 inhibitors in the lumen of the kidney's proximal tubules, and compared the efficacy of these compounds in terms of levels of *in vivo* SGLT2 inhibition.

Models for SGLT2 inhibitors were verified and validated against all available clinical data describing the PK of the compounds in plasma and urine. Similar to the model of Maurer et al. (2011), we introduced a peripheral compartment into our model of dapagliflozin. We also included this compartment for ipragliflozin and tofogliflozin models. Despite the inclusion of a peripheral compartment in a previously published model of empagliflozin (Riggs et al., 2013), we chose a model without peripheral compartment for empagliflozin based on Akaike Information Criterion (AIC) (Akaike, 1974). The AIC we calculated for the empagliflozin model without a peripheral compartment was lower than that of the model with peripheral compartment (343 vs. 345, respectively). Our canagliflozin model was also developed on the basis of one compartmental PK model.

Applying a verification strategy described in the Supplementary Materials, we were able to estimate the models' parameter values for active secretion of the drugs from plasma to urine and reabsorption from urine to plasma. Indeed, we found that the active secretion from plasma to the lumen of the

of different doses of dapagliflozin (Komoroski et al., 2009a). Model simulations are presented with 95% confidence bands.

**reabsorption inhibition levels during treatment with dapagliflozin.** Comparison of the average inhibition of glucose reabsorption mediated by

kidneys was equal to zero for dapagliflozin, canagliflozin, and ipragliflozin, but non-zero reabsorption should be taken into account to describe available clinical data. The reabsorption of empagliflozin was equal to zero, but non-zero secretion contributes substantially to the balance of the drug between plasma and urine. Both active secretion and reabsorption of tofogliflozin were equal to zero. Plasma and urine PK data are not sufficient to evaluate the unique values of the rate constants that are responsible for reabsorption and secretion simultaneously. These two parameters are correlated with each other. For tofogliflozin, there is lack of clinical data available to evaluate at least one of these parameters. To understand whether the uncertainty in these parameter values affects the description of clearance in the model, the total systemic and renal clearances measured in clinical trials were compared with those calculated in the model. The total systemic clearance of dapagliflozin was 265 mL/min in the model and 207 mL/min in clinical trials (Boulton et al., 2013). The renal clearance of dapagliflozin was 3.95 mL/min in the model and 3–5 mL/min in clinical trials (Komoroski et al., 2009a,b; Kasichayanula et al., 2011a). The renal clearances of empagliflozin and ipragliflozin in the model were 35 and 1.88 mL/min respectively, while they were 30–37 mL/min for empagliflozin (Heise et al., 2013a,b) and 1–3 mL/min for ipragliflozin (Veltkamp et al., 2011; Zhang et al., 2013) in clinical data. The total systemic clearance of tofogliflozin was equal to 9.9 L/h in the model and 10 L/h in clinical trials (Schwab et al., 2013). The renal clearance of tofogliflozin was equal to 20 mL/min in the model and 25.7 mL/min in clinical trials (Schwab et al., 2013). Therefore, the values produced by our models are similar to experimentally measured values.

Comparison of the average inhibition of glucose reabsorption mediated by

presented with 95% confidence bands.

**FIGURE 11 | Levels of SGLT2 inhibition after drug administration.** Levels of SGLT2 inhibition following simulation of multiple administrations of labeled doses of marketed SGLT2 inhibitors and maximal doses of other SGLT2 inhibitors approved for phase 2/3 studies. Colors of curves correspond to different compounds: black—10 mg QD dapagliflozin; blue—300 mg QD canagliflozin; red—25 mg QD empagliflozin; green—300 mg QD ipragliflozin; pink—40 mg QD tofogliflozin. Model simulations are presented with 95% confidence bands.

One hypothesis proposed by Liu et al. (2012) states that the low efficacy of SGLT2 inhibitors could be explained by the low concentration of compounds at the potential site of action–the lumen of kidney's proximal tubules. Our model suggests that concentrations of SGLT2 inhibitors in the lumen are higher than the corresponding unbound and total concentrations in plasma. This can be explained in terms of reabsorption of water in kidney's proximal tubules (Panchapakesan et al., 2009). A decrease in water volume in the lumen leads to an increase in the concentration of the compounds in lumen than in the plasma. In the model, this is reflected by the ratio of glomerular filtration and urine excretion (urine formation) rates. Thus, we conclude that the concentration of compounds in the lumen of the kidney's proximal tubules is actually relatively high, and we reject the hypothesis that the low efficacy of SGLT2 inhibitors is caused by their low concentration at the site of action. Our model predicts that concentrations of empagliflozin and tofogliflozin in kidney lumen are higher than concentrations of dapagliflozin, canagliflozin, and ipragliflozin following administration of equal doses. This prediction is supported by data obtained from clinical trials, which indicates that approximately 20% of the empagliflozin and tofogliflozin dose is recovered in urine (Brand et al., 2012; Zell et al., 2013) in comparison to about 1% for the other SGLT2 inhibitors (Kasichayanula et al., 2011a; Devineni et al., 2013; Zhang et al., 2013).

The model shows that administration of labeled doses of marketed inhibitors (dapagliflozin and canagliflozin) and of maximal doses approved for phase 2/3 studies of tofogliflozin, ipragliflozin, and empagliflozin leads to almost complete inhibition of SGLT2 *in vivo*. The average inhibition level of glucose reabsorption mediated by SGLT2 (the average SGLT2 inhibition level) predicted by the model is higher than glucose reabsorption inhibition level measured in clinical trials, and almost equal to 100%. Thus, if the lumen of the kidney's proximal tubules is indeed the site of action of SGLT2 inhibitors, and SGLT2 is the main transporter that facilitates reabsorption of glucose, then there is a contradiction between the simulations produced by the model and the clinical data. This contradiction leads us to support two hypotheses proposed by Liu et al. (2012). Firstly, that the potential site of action of SGLT2 inhibitors is not in the lumen of the kidney's proximal tubules, but in the proximal tubule cells. Secondly, that there are other transporters that could facilitate glucose reabsorption under the conditions of SGLT2 inhibition. Both hypotheses appear reasonable within the framework of our current model. To determine which of these two hypotheses is true, further modeling supported by additional experimental data is required.

The model was applied to compare the efficacy of different SGLT2 inhibitors in respect of inhibiting other transporters in the SGLT family that are expressed in kidneys. Simulating the administration of equal doses of compounds (20 mg), labeled doses of marketed inhibitors, and maximal doses of other inhibitors approved for phase 2/3 studies showed that levels of SGLT2 inhibition are similar for all compounds and almost equal to 100%. This is caused by high concentrations of SGLT2 inhibitors in the lumen of the kidneys and low *in vitro* IC50 values. The effect of SGLT2 inhibitors on other transporters is also rather strong because of its high concentration in kidney's lumen. Canagliflozin is the strongest SGLT1 and SGLT6 inhibitor, and this is caused by canagliflozin's low IC50 (the lowest of all the SGLT2 inhibitors considered here). The most effective SGLT3 inhibitors are ipragliflozin and tofogliflozin because of their high *in vitro* potencies. All compounds have similar effects on SGLT4, except dapagliflozin. Dapagliflozin has the weakest influence on SGLT4 because of its low concentration in the kidney's lumen. Ipragliflozin and empagliflozin are the strongest SGLT5 inhibitors because of their low IC50 of for SGLT5, and the high concentration of ipragliflozin in the lumen.

In conclusions, we have developed a systems pharmacology model for SGLT2 inhibitors that enables the estimation of its concentration in the lumen of the kidney's proximal tubules (the potential site of SGLT2 action) and the prediction of SGLT2 inhibition levels during treatment in humans. We have shown that the concentration of SGLT2 inhibitors in the lumen of the proximal tubules is high, and that the level of SGLT2 inhibition during treatment in humans is almost 100%. Based on the results of our model, two explanations for the observed low efficacy of SGLT2 inhibitors were supported: (1) the site of action of SGLT2 inhibitors is not in the lumen of the kidney's proximal tubules, but elsewhere (e.g., the kidneys proximal tubule cells); and (2) there are other transporters that could facilitate glucose reabsorption under the conditions of SGLT2 inhibition (e.g., other transporters of SGLT family). It was found that all SGLT2 inhibitors have similar efficacy.

## **SUPPLEMENTARY MATERIAL**

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

## **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.

*Received: 30 April 2014; accepted: 09 September 2014; published online: 13 October 2014.*

*Citation: Demin O Jr., Yakovleva T, Kolobkov D and Demin O (2014) Analysis of the efficacy of SGLT2 inhibitors using semi-mechanistic model. Front. Pharmacol. 5:218. doi: 10.3389/fphar.2014.00218*

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

*Copyright © 2014 Demin, Yakovleva, Kolobkov and Demin. 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.*

## Use of systems pharmacology modeling to elucidate the operating characteristics of SGLT1 and SGLT2 in renal glucose reabsorption in humans

#### *Yasong Lu1 \*, Steven C. Griffen2†, David W. Boulton3 and Tarek A. Leil <sup>1</sup>*

*<sup>1</sup> Quantitative Clinical Pharmacology, Clinical Pharmacology and Pharmacometrics, Exploratory Clinical and Translational Research, Bristol-Myers Squibb, Princeton, NJ, USA*

*<sup>2</sup> Diabetes Development Center, Global Clinical Research, Bristol-Myers Squibb, Princeton, NJ, USA*

*<sup>3</sup> Clinical Pharmacology and Pharmacometrics, Exploratory Clinical and Translational Research, Bristol-Myers Squibb, Princeton, NJ, USA*

#### *Edited by:*

*Maria Angela Sortino, University of Catania, Italy*

#### *Reviewed by:*

*Volker Vallon, University of California San Diegio and VA San Diego Healthcare System, USA Robinson Sabino-Silva, Federal University of Uberlandia, Brazil*

#### *\*Correspondence:*

*Yasong Lu, Clinical Pharmacology and Pharmacometrics, Exploratory Clinical and Translational Research, Bristol-Myers Squibb, MS# E14-08, Route 206 and Province Line Road, Princeton, NJ 08543, USA e-mail: yasong.lu@bms.com*

### *†Present address:*

*Steven C. Griffen, JDRF, New York, USA*

In the kidney, glucose in glomerular filtrate is reabsorbed primarily by sodium-glucose cotransporters 1 (SGLT1) and 2 (SGLT2) along the proximal tubules. SGLT2 has been characterized as a high capacity, low affinity pathway responsible for reabsorption of the majority of filtered glucose in the early part of proximal tubules, and SGLT1 reabsorbs the residual glucose in the distal part. Inhibition of SGLT2 is a viable mechanism for removing glucose from the body and improving glycemic control in patients with diabetes. Despite demonstrating high levels (in excess of 80%) of inhibition of glucose transport by SGLT2 *in vitro*, potent SGLT2 inhibitors, e.g., dapagliflozin and canagliflozin, inhibit renal glucose reabsorption by only 30–50% in clinical studies. Hypotheses for this apparent paradox are mostly focused on the compensatory effect of SGLT1. The paradox has been explained and the role of SGLT1 demonstrated in the mouse, but direct data in humans are lacking. To further explore the roles of SGLT1/2 in renal glucose reabsorption in humans, we developed a systems pharmacology model with emphasis on SGLT1/2 mediated glucose reabsorption and the effects of SGLT2 inhibition. The model was calibrated using robust clinical data in the absence or presence of dapagliflozin (DeFronzo et al., 2013), and evaluated against clinical data from the literature (Mogensen, 1971; Wolf et al., 2009; Polidori et al., 2013). The model adequately described all four data sets. Simulations using the model clarified the operating characteristics of SGLT1/2 in humans in the healthy and diabetic state with or without SGLT2 inhibition. The modeling and simulations support our proposition that the apparent moderate, 30–50% inhibition of renal glucose reabsorption observed with potent SGLT2 inhibitors is a combined result of two physiological determinants: SGLT1 compensation and residual SGLT2 activity. This model will enable *in silico* inferences and predictions related to SGLT1/2 modulation.

**Keywords: systems pharmacology model, SGLT, dapagliflozin, renal glucose reabsorption, glucosuria, diabetes mellitus**

## **INTRODUCTION**

In the kidney, plasma glucose is continuously filtered by glomeruli and reabsorbed along the proximal tubules. Under normal physiological conditions, the reabsorption is almost complete. The reabsorption is mediated primarily by two sodium-glucose cotransporters (SGLTs), SGLT1 and SGLT2. In the kidney, SGLT2 is located in the early part (S1/S2 segments) of the proximal tubules, and is recognized as a low affinity, high capacity pathway for renal glucose reabsorption. SGLT1, on the other hand, is located in the distal part (S3 segment) of the proximal tubules, and is characterized as a high affinity, low capacity pathway (Wright, 2001; Mather and Pollock, 2011). SGLT2 is believed responsible for 80–90% of renal glucose reabsorption, and SGLT1 for the rest (10–20%) in healthy humans under normal physiological conditions (DeFronzo et al., 2012). SGLT2 has been identified as a viable target for improving glycemic control in diabetes. Two potent and selective SGLT2 inhibitors, dapagliflozin and canagliflozin, have been approved by many regulatory agencies for treating type 2 diabetes mellitus (T2DM).

Given the overwhelming contribution (*>*80%) of SGLT2 to renal glucose reabsorption, it has been expected that SGLT2 inhibitors, at sufficient exposures, would reduce renal glucose reabsorption by over 80%. This expectation, however, appeared to be contradicted by the clinical observations that only 30–50% of inhibition in glucose reabsorption was achieved with dapagliflozin and canagliflozin (Komoroski et al., 2009a; Devineni et al., 2013; Washburn and Poucher, 2013). To explain this apparent contradiction, several hypotheses, from peculiar pharmacokinetics of an inhibitor in the kidney (Liu et al., 2012) to SGLT1 compensation (Haddish-Berhane et al., 2010; Maurer et al., 2011; Abdul-Ghani et al., 2013), have been proposed. These hypotheses are yet to be tested.

Recently, Hummel et al. (2011) used a quantitative *in vitro* electrophysiological study to generate hypotheses about the relative contributions of human SGLTs to renal glucose reabsorption. Hummel et al. found that the two human transporters have similar apparent affinity for D-glucose (5 mM for hSGLT2 vs. 2 mM for hSGLT1), and inferred that the capacity of hSGLT1 for renal glucose reabsorption may be over 50% of hSGLT2 under normal conditions in humans. As such, the difference in the contribution to renal glucose reabsorption between the two cotransporters may be less profound than previously perceived.

Despite a large body of research in SGLTs and renal glucose reabsorption, the quantitative understanding of the characteristics of these cotransporters in humans remains limited (Vallon, 2011). Assessments in this regard have largely relied on fragments of data, insufficient to account for all key variables (e.g., SGLTs activities, plasma glucose levels, pharmacokinetic profiles of SGLTs inhibitors), and empirical, static mathematical models that do not account for the dynamic processes of renal glucose filtration, reabsorption, and transfer along tubular lumen over time. Consequently, a quantitative, holistic characterization has not yet been formulated.

Systems pharmacology modeling is a powerful tool for data and knowledge integration and hypothesis testing, and for providing quantitative understanding of a pharmacological target or pathway and insights into "what-if" scenarios that may not be feasibly obtained experimentally. For SGLTs-mediated renal glucose reabsorption, Yamaguchi et al. reported simplified systems pharmacology models in mice (Yamaguchi et al., 2012) and rats (Yamaguchi et al., 2011), and Haddish-Berhane et al. (2010) presented a conference poster on a minimal systems pharmacology model in humans with limited evaluation against clinical data on dapagliflozin (Komoroski et al., 2009a).

This report presents a systems pharmacology model that was developed based on renal physiology and a robust clinical data set, with emphasis on SGLTs-mediated glucose reabsorption in the proximal tubules. The model was evaluated against several external clinical data sets. It is anticipated that the model will be valuable in:


## **MATERIALS AND METHODS**

### **STUDIES AND DATA SETS**

The studies and data sets used for model calibration and evaluation are listed in **Table 1**. For more details, the reader is referred to the original reports.

The DeFronzo et al. (2013), Polidori et al. (2013) and Wolf et al. (2009) studies employed stepped hyperglycemic clamp (SHC) procedures, and the Mogensen study (1971) was conducted at fixed, elevated plasma glucose levels. The clinical approach of artificially maintaining a constant plasma glucose concentration allowed us to ignore the potential impacts of renal glucose reabsorption on plasma glucose concentration, hence simplifying the process of model development. Simulations using the systems pharmacology model with fixed glucose levels will provide "clean" illustrations of SGLTs operating characteristics. A more comprehensive model integrating renal glucose reabsorption and glucose-insulin homeostasis will be reported elsewhere (Lu et al., 2014).

The mean data from each study were used for model calibration or evaluation. The data in DeFronzo et al. (2013) were available from an internal database owned by Bristol-Myers Squibb/AstraZeneca. We excluded from analysis those data points where the actual plasma glucose level deviated 25% or more from the corresponding group means. These data points appeared at the steps with target glucose level ≥450 mg/dL, and represented only 17% of total data points at those steps. This exclusion should abolish potential undue influences of excessive variability in the data on parameter estimation.

## **MODEL STRUCTURE**

The model structure, shown in **Figure 1**, was developed based on the renal physiology and pharmacological understanding of SGLTs inhibition. The model describes the disposition of glucose as well as SGLTs inhibitors, if applicable, with emphasis on glomerular filtration and tubular reabsorption. The proximal convoluted tubules (PCT) were divided equally into six sequential sub-segments (PCT1-6), and the proximal straight tubules (PST) were divided equally into three sub-segments (PST1-3). The division allowed a more accurate description of the luminal glucose concentration as the filtrate progresses through tubular segments and the amount of UGE over time. The number of sub-segments was chosen to achieve an approximate agreement between predicted and observed UGE in a healthy subject under normal conditions. The distal tubules were not included due to their irrelevance to glucose reabsorption. The glomerular filtrate flowed from PCT1 through PST3 and drained into the urinary bladder. A urine compartment was added for collecting urine and urinary glucose.

Along the proximal tubules, filtered glucose was continuously reabsorbed. It was assumed that the absorption was mediated by SGLT2 in the PCT (PCT1-6) and by SGLT1 in the PST (PST1-3). The maximum reabsorption rate of SGLT2 (Vmax2) was uniformly distributed among the PCT1-6 sub-segments, and likewise for the Vmax of SGLT1 (Vmax1) among the PST1-3 subsegments. In each sub-segment, glucose was reabsorbed via a Michaelis-Menten process as Equation (1):

$$R\_{\dot{\jmath}} = \frac{V\_{\text{max},\dot{\jmath}} \times C\_{\text{glu},\dot{\jmath}}}{K\_m + C\_{\text{glu},\dot{\jmath}}} \tag{1}$$

where the subscript j is an index for a tubular sub-segment, and for a given sub-segment, *R* is the glucose reabsorption rate



(mass/time), *K*<sup>m</sup> denotes glucose affinity for SGLT1 or SGLT2, and *Cglu* represents luminal glucose concentration.

The mass of reabsorbed glucose was directed to another compartment (glucose reabsorbed) instead of the plasma glucose compartment. This approach is appropriate for scenarios where renal glucose recovery does not affect plasma glucose level, such as: (1) experimental procedures that fix plasma glucose levels (Mogensen, 1971; Wolf et al., 2009; DeFronzo et al., 2013; Polidori et al., 2013), and (2) subjects with normal glucose tolerance who can efficiently dispose the absorbed mass to maintain plasma glucose constant at the fasting state.

In the case where an SGLTs inhibitor, e.g., dapagliflozin or canagliflozin, was administered, the unbound portion of the inhibitor in plasma was freely filtered via glomeruli. The inhibitor then traveled through the tubular sub-segments and the urinary bladder, and was excreted to the urine compartment, similar to glucose but without tubular reabsorption. Within each subsegment, the inhibitor competed with glucose for SGLT1/2, and hence competitively inhibited glucose reabsorption. The reabsorption rate (*R*∗ *<sup>j</sup>* ), with competitive inhibition of an inhibitor, in a given sub-segment became:

$$R\_j^\* = \frac{V\_{\text{max},j} \times C\_{\text{glu},j}}{K\_m \times \left(1 + \frac{C\_{d \text{reg},j}}{K\_i}\right) + C\_{\text{glu},j}} \tag{2}$$

where *Cdrug* denotes luminal SGLTs inhibitor concentration, and *K*<sup>i</sup> is the affinity of the inhibitor to SGLTs. [See Supplementary Materials for further expansion of Equations (1) and (2)].

To obtain a similar time course of plasma concentration of canagliflozin with the dosing regimen (100 mg QD for 8 days) described in Polidori et al. (2013), a two-compartment pharmacokinetic (PK) model was developed based on the mean PK data reported (Devineni et al., 2013). For dapagliflozin, observed mean PK after 10 mg, QD, for 7 days was reported in DeFronzo et al. (2013). Interpolation of the observed dapagliflozin PK provided an input into the plasma inhibitor compartment (**Figure 1**) to allow description of dapagliflozin inhibition of tubular glucose reabsorption.

## **MODEL PARAMETERS AND CALIBRATION**

The physiological parameters, such as volumes, flow rates, glucose affinity for SGLTs (Km), and glucose reabsorption capacities (Vmax) are listed in **Table 2**, and SGLT2 inhibitor physicochemical parameters and binding affinity for SGLTs (Ki) are in **Table 3**. Most of the parameters were from the literature, measured in each of the respective studies, or based on reasonable assumptions, except for Vmax, Km, and Ki, whose values were calibrated. For parameter calibration, literature values were taken as the starting points (see **Tables 2, 3**), and then fine-tuned to allow the model predictions to be consistent with the mean UGE data in DeFronzo et al. (2013). Because a satisfactory agreement between the predictions and the observations of UGE could not be achieved over the entire plasma glucose range of 100–550 mg/dL, the calibration was focused on the data in the clinically relevant range, 100–400 mg/dL. In the end, only the set of calibrated values were considered physiologically plausible and accepted if it adequately described the DeFronzo et al. (2013) as well as the other three data sets (Mogensen, 1971; Wolf et al., 2009; Polidori et al., 2013).

The potential influences of diabetes and SGLT2 inhibition on the parameters to be calibrated were considered during parameterization. Renal SGLTs expression and activity may change in response to SGLT2 inhibition and/or diabetes. In the wildtype mouse, SGLT2 protein expression was enhanced with the treatment of empagliflozin, a selective SGLT2 inhibitor, without upregulation of mRNA (Vallon et al., 2014). In the diabetic state, the expressions of SGLT2 mRNA and protein have been found upregulated significantly relative to the respective controls in genetically modified mice (Vallon et al., 2014), diabetic rats, (Freitas et al., 2008; Tabatabai et al., 2009), and humans (Rahmoune et al., 2005). For renal SGLT1, however, the response is more diverse, with increased, unchanged, or reduced expression and/or activity observed in animals (Vallon and Thomson, 2012; Vallon et al., 2014). It is challenging to incorporate these potential changes in SGLT activity in the model for two reasons: (i) limited quantitative understanding in humans regarding these changes, and (ii) adequate calibration of parameters for these changes is not supported by available data. For simplification, therefore, the Vmax1, Km, and Ki values were assumed consistent between the healthy and diabetics, and Vmax2 was allowed to adjust between the healthy and disease state. The Vmax2 in healthy subjects was estimated as a proportion of that in diabetics, and the value of the proportion was calibrated using the DeFronzo et al. (2013). The potential impact of SGLT2 inhibition on Vmax, Km, and Ki values was ignored.

Although Vmax in humans has generally been reported as the sum of Vmax1 and Vmax2, with difficulty in separating the two components, it is worth pointing out that in our study, quantitative separation of Vmax1 and Vmax2 was feasible without an assumption of the value of Vmax1/Vmax2 ratio, because the calibration data set (DeFronzo et al., 2013) encompassed scenarios with and without perturbation of SGLT2 activity. Such a separation was achieved previously in rats with the aid of mathematical modeling (Yamaguchi et al., 2011).

## **MODEL EVALUATION**

Once it was calibrated using the DeFronzo et al. (2013), the model was evaluated for its predictivity against three data sets from different sources (Mogensen, 1971; Wolf et al., 2009; Polidori et al., 2013). The parameters were held constant for the evaluation unless they were study specific, in which case they were adjusted per the study conditions as listed in **Tables 1, 2**. The Ki values of canagliflozin, necessary for simulating the Polidori et al., conditions (Polidori et al., 2013), are listed in **Table 3**.

#### **SIMULATIONS AND EXPLORATIONS**

## *Renal glucose reabsorption and UGE vs. loss-of-function mutation of SGLTs*

Numerous mutations in SGLT1 (Martin et al., 1996; Lam et al., 1999) and SGLT2 have been identified in humans (Santer et al., 2003; Kleta et al., 2004; Calado et al., 2008; Yu et al., 2011). The mutations in SGLT1 disrupt the trafficking of SGLT1 from the endoplasmic reticulum to the plasma membrane (Lam et al., 1999), and the mutations in SGLT2 reduce SGLT2 expression in the apical side of PCT (Yu et al., 2011). These mutations are likely to reduce the Vmax of these cotransporters. It is yet to be clarified to what extent the function of SGLTs in the kidney is affected by a given mutation. Simulations using our systems pharmacology model can provide theoretical, quantitative relationships between a reduction in Vmax and glucose reabsorption or UGE in an otherwise healthy person. To enable these simulations, the mean daily plasma glucose profile in the healthy subjects from Freckmann et al. (2007) was used as an input to the plasma glucose compartment of our model. The quantitative SGLTs function-UGE relationships will be instrumental to mapping renal SGLTs genotypes to their apparent functions.

## **Table 2 | Physiological parameters.**


## *Sensitivity of renal glucose reabsorption and UGE to SGLT1 kinetics (Vmax 1 and K<sup>i</sup> 1)*

The analysis of sensitivity of renal glucose reabsorption and UGE to SGLT1 kinetics will help clarify these questions: (1) How strong is the influence of an alteration of SGLT1 kinetics on renal glucose reabsorption and UGE in the healthy state? (2) How strong is the influence in the diabetic state? (3) From drug discovery perspective, without consideration of its effect on intestinal SGLT1, will an SGLT1/2 dual inhibitor induce stronger glucosuria than a highly selective SGLT2 inhibitor, e.g., dapagliflozin? The analysis was conducted with simulations in a naive healthy subject and a T2DM subject with or without SGLT2 inhibition under the SHC procedure used by DeFronzo et al. (2013) with a target plasma glucose range of 100–350 mg/dL. With all other parameters held constant, we first evaluated how a decrease in Vmax1 would affect renal glucose reabsorption and UGE; and likewise, we then evaluated how changes in Km1, or Ki1 with the presence of SGLTs inhibition, would affect renal glucose reabsorption and UGE.

## **SOFTWARE**

Processing of the raw data from the DeFronzo et al. study (2013) was conducted using S-PLUS 8.1 version 3.4 (TIBCO, Palo Alto,

**Table 3 | SGLT2 inhibitor-specific parameters.**


*\*Reference: Bristol-Myers Squibb/AstraZeneca report (2010): Summary of clinical pharmacology studies: Dapagliflozin (BMS-512148). BMS Document Control Number 930047848.*

CA) on a UNIX platform. Model development and simulations were performed using Berkeley Madonna version 8.3.18 (Berkeley Madonna Inc., Berkeley, CA).

## **RESULTS**

## **MODEL CALIBRATION USING DeFronzo et al. (2013)**

The UGE data from DeFronzo et al. (2013) allowed estimation of Vmax, glucose Km, and dapagliflozin Ki values for SGLTs in the healthy subjects and T2DM patients. In the model, the healthy and T2DM subjects were differentiated by their Vmax2 for describing the DeFronzo et al. conditions. The parameter estimates are presented in **Tables 2, 3**. The model performance is demonstrated in **Figure 2**. The model adequately described the cumulative (**Figures 2A,B**) and step-wise UGE data (**Figures 2C,D**) at baseline and in the first 4 h (where the target plasma glucose escalated from 100 to 350 mg/dL) after dapagliflozin treatment. From 4.67 h onward (where the target plasma glucose increased from 400 to 550 mg/dL), the model prediction of UGE in the dapagliflozin treated groups was slightly lower than the observed. The glucose concentrations in the tubular sub-segments PCT1-6 and PST1-3 in T2DM patients at baseline and treated with dapagliflozin are illustrated in Figure S1. At baseline, the tubular glucose concentration tapers along the proximal tubules with plasma glucose level up to 23 mM. With further increase in plasma glucose, the tubular glucose level becomes more uniform as the reabsorption approaches saturation. After dapagliflozin treatment, however, glucose is increasingly concentrated along the proximal tubules.

### **MODEL EVALUATION**

The calibrated model was evaluated for its predictive performance relative to three separate clinical data sets (Mogensen, 1971; Wolf et al., 2009; Polidori et al., 2013). The predictions are overlaid with corresponding observations in **Figure 3**. The predictions agreed well with the observed data, indicating that the model is plausible, has reasonable accuracy, and can be used for inference and prediction.

## **SGLTs OPERATING CHARACTERISTICS FOR THE DeFronzo et al. (2013) CONDITIONS**

### *SGLTs relative contributions to renal glucose reabsorption*

The model derived step-wise amount of glucose reabsorbed by renal SGLT1 and SGLT2 in the healthy subjects at baseline and after dapagliflozin treatment is shown in **Figures 4A,B**, and the relative contributions of the two pathways at each step are in **Figures 4C,D**. At near normal glycemic levels (∼100 mg/dL) at baseline (without SGLT2 inhibition), SGLT2 contributed to 90% of total reabsorption and SGLT1 10%. The 90%/10% split became 80%/20% with plasma glucose escalated to over 200 mg/dL. With the presence of dapagliflozin, the contribution of SGLT2 declined and SGLT1 became the more predominant reabsorption pathway; the relative contributions varied with plasma dapagliflozin concentration over time. Similar results were obtained in the T2DM patients, for whom the relative contributions of SGLT1 and SGLT2 before and after dapagliflozin treatment are illustrated in **Figures 4E,F**.

## *SGLTs operation efficiency*

The calculated operation efficiency (defined as reabsorption rate/Vmax × 100% for either SGLT1 or SGLT2) for both SGLTs in the healthy subjects is plotted in **Figure 5A**. At the plasma glucose level of ∼100 mg/dL, SGLT2 and SGLT1 were operating at ∼40 and 20% of their respective Vmax. The operation efficiency increased with plasma glucose (and thereby filtered glucose load) for both SGLTs, with the slope for SGLT1 being much steeper than for SGLT2. The operation efficiency at plasma glucose ≥400 mg/dL reached 97% for SGLT1 and 81% for SGLT2. Even with plasma glucose as high as 550 mg/dL, SGLT2 operated at just 89% of its capacity. Dapagliflozin treatment lowered SGLT2 operation efficacy to as low as 10%, and drove SGLT1 operation to over 90% of its capacity. Similar results were found in the T2DM patients (**Figure 5B**).

## **SIMULATIONS AND EXPLORATIONS**

## *Renal glucose reabsorption and UGE vs. loss-of-function mutation of SGLTs*

Simulations were conducted to establish quantitative relationships between loss of function (i.e., reduction in Vmax) of SGLT2 or SGLT1 and renal glucose reabsorption as well as UGE in an otherwise healthy subject with normoglycemia (plasma glucose ranging from 80 to 125 mg/dL with a time-weighted average of 90 mg/dL). The simulation results for SGLT2 are in **Figure 6A** and SGLT1 in **Figure 6B**. A 50% loss of function for SGLT2 caused UGE of 4.5 g per day, and 100% of loss of function resulted in 79 g UGE per day. The total glucose reabsorption was lowered by 17, 32, and 49% for a 75, 87.5, and 100% loss of SGLT2 function, respectively. Loss of SGLT1 function caused much less severe glucosuria, 1.2 g at 50% and 16 g at 100% of loss of function. The total glucose reabsorption was reduced by only 10% with complete loss of SGLT1 activity.

### *Sensitivity of renal glucose reabsorption and UGE to SGLT1 kinetics*

In a naive healthy subject with a plasma glucose level of ∼100 mg/dL, the elimination of SGLT1 activity, either through driving Vmax1 to zero or Km1 to infinity, inhibited renal glucose

reabsorption by only 10%. This was consistent with the result above for SGLT1 loss-of-function mutation (**Figure 6B**). In the diabetic state with a mean plasma glucose level up to 250 mg/dL, simulations suggested a slightly stronger influence (up to 20% lowering) on renal glucose reabsorption.

For simulations in the diabetic state with the presence of SGLT2 inhibition, the inhibitor was assumed identical to dapagliflozin, except that its Ki1 was allowed to change. The sensitivity of UGE to Vmax1 at 10, 14, 17, and 20 mmol/h (corresponding to a 50, 30, 15, and 0% reduction in SGLT1 capacity) is illustrated in **Figure 7A**. A mild to moderate, depending on the glucose level, increase in UGE was expected with reduction in Vmax1. At the plasma glucose level of 150–250 mg/dL, roughly equivalent to the range of average levels in real-life T2DM patients, a 30% reduction in Vmax1, presumed to be clinically well tolerated (Abdul-Ghani et al., 2013), augmented glucosuria by up to 30%. The sensitivity of UGE to Ki1 is depicted in **Figure 7B**. The tested Ki1 values ranged from 6 to 10,000 nM. The 6 nM represented a 20× selectivity (similar to the SGLT1/2 dual inhibitor LX4211 Zambrowicz et al., 2013) for SGLT2 (0.3 nM) vs. SGLT1 (6 nM) for an SGLTs inhibitor which is otherwise identical to dapagliflozin. UGE was found to be insensitive to Ki1.

## **DISCUSSION**

Highly selective and potent SGLT2 inhibitors, such as dapagliflozin and canagliflozin, have demonstrated significant and clinically meaningful effects on glycemic control in T2DM patients. It has been puzzling that SGLT2 inhibitors inhibit renal glucose reabsorption by only 30–50% clinically, despite the overwhelming contribution of SGLT2 to renal glucose reabsorption (80–90%) under normal conditions. Several hypotheses have been proposed for this apparently discrepant observation (Haddish-Berhane et al., 2010; Maurer et al., 2011; Pfister et al., 2011; Liu et al., 2012; Abdul-Ghani et al., 2013), and most of them are focused on the compensatory effect of SGLT1. The hypothesis of SGLT1 compensation has recently been confirmed in mice (Rieg et al., 2014). However, due to the differences in the experimental conditions in mice (Rieg et al., 2014) and in clinical trials (Komoroski et al., 2009a; Devineni et al., 2013; Heise et al., 2013; Washburn and Poucher, 2013) (see Table S1), extrapolation of the Rieg et al. (2014) finding to the clinic is not straightforward. As a whole, this situation indicates that, despite tremendous advances in the basic biology of SGLTs and pharmaceutical development targeting SGLT2, the roles of these transporters in renal glucose reabsorption, especially in humans, have yet to be clarified in a quantitative, mechanistic manner. To this end, we developed a systems pharmacology model for SGLT-mediated renal glucose reabsorption in humans with or without pharmacological modulation of SGLT2 activity. In general, this model adequately described four separate data sets from different study settings (Mogensen, 1971; Wolf et al., 2009; DeFronzo et al., 2013; Polidori et al., 2013) and replicated severe

glucosuria (79 g/day) in normoglycemic human subjects with homozygous SLGT2 mutations (Santer et al., 2003).

The prediction of UGE at the plasma glucose level of 400 mg/dL and higher in the subjects treated with dapagliflozin was lower than the observed data from DeFronzo et al., (**Figure 2**). This possibly results from compensatory effects in the renal tubules when glucose concentrations are drastically elevated. Bank and Aynedjian (1990) proposed that high glucose concentration in the proximal tubules would stimulate water reabsorption in the proximal portion and enhance compensatory water excretion in the more distal portion. In the DeFronzo et al. study (2013), an increase in urine volume was observed with escalation of plasma glucose level. This hydrodynamic change in response to glucose level may interfere with tubular glucose reabsorption. These processes, however, were not included in the model. Nevertheless, the unsatisfactory performance at high glucose levels (over 400 mg/dL) is unlikely to hamper the utility of the model because those glucose levels are irrelevant to most of normal or even diabetic conditions. Overall, the performance of the model suggests that the model is useful for mechanistically evaluating the roles of SGLT1 and SGLT2 in renal glucose reabsorption, and for predicting clinical pharmacodynamics of SGLT2 inhibitors.

## **CHARACTERIZATION OF SGLTs OPERATION WITHOUT THE PRESENCE OF SGLT2 INHIBITION**

The Vmax values of SGLT1 and SGLT2 were estimated to be 20 mmol/h and 94 (healthy)/110 (diabetic) mmol/h, respectively, and the glucose Km values for SGLT1 and SGLT2 were estimated to be 0.5 and 4 mM, respectively. The sum of the Vmax values (i.e., total reabsorption capacity) and the two Km estimates are similar to previously reported estimates (Mogensen, 1971; Diez-Sampedro et al., 2001; Chao and Henry, 2010; Hummel et al., 2011; DeFronzo et al., 2013). The Vmax2/Vmax1 ratio in the healthy subject (4.7) is consistent with that in the rat (5.4) (Yamaguchi et al., 2011). These estimates reinforce the concept of SGLT1 being a high affinity, low capacity transporter and SGLT2 being a low affinity, high capacity transporter for renal glucose reabsorption.

Under near normoglycemic conditions (average plasma glucose ∼80–120 mg/dL) in both healthy and diabetic subjects, SGLT2 and SGLT1 are operating at about 40 and 20% of their respective capacities, and contributing to 90 and 10% of total glucose reabsorption, respectively. With the increase in plasma glucose concentration, SGLT2 operation efficiency steadily increases to near 90% of its capacity, whereas SGLT1 operation efficiency jumps sharply to over 80% of its capacity

and then steadily approaches to the maximum. The relative contributions of 90%/10% gradually becomes 80%/20% for SGLT2 and SGLT1 as plasma glucose rises. These results solidify the current characterization of the relative contributions of the two transporters to renal glucose reabsorptoin without the presence of SGLT2 inhibition (Chao and Henry, 2010; DeFronzo et al., 2012).

## **CHARACTERIZATION OF SGLT's OPERATION IN THE PRESENCE OF SGLT2 INHIBITION**

With the treatment of dapagliflozin at its clinical dose (10 mg QD), the majority of SGLT2 is occupied by dapagliflozin molecules (occupancy up to 98% at the peak exposure). The total activity of SGLT2 in a healthy or T2DM subject is suppressed considerably, from ∼40% of operation efficiency without SGLT2 inhibition to only 10% with the treatment of dapagliflozin. Consequently, the contribution of SGLT2 to renal glucose reabsorption declines from 80 to 90% at baseline to less than 50% with SGLT2 inhibition. Meanwhile, the importance of SGLT1 to renal glucose reabsorption jumps sharply. The operation efficiency of SGLT1 reaches over 90%, up from 20% at baseline. As a result, SGLT1 accounts for over 50% of renal reabsorption when SGLT2 is inhibited, much higher than the 10–20% at baseline.

**FIGURE 5 | Model derived operation efficiency (defined as glucose reabsorption rate/Vmax × 100% for either SGLT1 or SGLT2) for both SGLTs at baseline and after dapagliflozin**

**treatment in the healthy subjects (A) and patients with diabetes (B) under the stepwise hyperglycemic clamp procedure in DeFronzo et al. (2013).**

## **THEORETICAL MAXIMUM INHIBITION OF RENAL GLUCOSE REABSORPTION**

The simulations of loss of function of SGLTs (reduction in Vmax) vs. renal glucose reabsorption provide a clean relationship for assessing the theoretical maximum inhibition of the reabsorption. In a healthy subject under physiological conditions, an 87.5– 100% loss of SGLT2 function results in a 32–49% of inhibition of renal glucose reabsorption. In a diabetic patient, the glucose reabsorption vs. loss of SGLT2 activity curve shifts downwards, i.e., somewhat greater inhibition of reabsorption. With a daily average plasma glucose level of 150 mg/dL, a complete loss of SGLT2 activity lowers the reabsorption by 70%. This greater extent of inhibition in diabetics is due to the up-regulated activity of SGLT2 in the disease state (Rahmoune et al., 2005).

The loss of SGLT1 function has only mild inhibitory effect on renal glucose reabsorption. An entire loss of SGLT1 function leads to only 10% of inhibition of glucose reabsorption in a normoglycemic healthy subject and up to 15% of inhibition in a diabetic patient with a daily average plasma glucose level of 150 mg/dL.

Our results of theoretical maximum inhibition of renal glucose reabsorption due to loss of activity of SGLT1 or SGLT2 are in general agreement with the findings in Sglt1/2 knock-out mice. In Sglt2−*/*<sup>−</sup> mice the renal glucose reabsorption is reduced to ∼50% of that in wild-type mice at euglycemia, and is further reduced with increase in filtered glucose load (Vallon et al., 2011). The knock-out of Sglt1−*/*<sup>−</sup> in mice causes a 2–3% decrease in total renal glucose reabsorption (Gorboulev et al., 2012; Powell et al., 2013). The numerical discrepancy in the maximum influence of SGLT1 loss (2–3% in mice vs. 10% in humans) is yet to be understood. It may reflect a real inter-species difference in the contribution of SGLT1, a result secondary to inter-species differences in other physiological factors, or an inter-study variation as well as random errors. Extension of our systems pharmacology model to mice with appropriate physiological parameters could shed light on this issue.

## **EXPLANATION TO THE PUZZLING MODERATE INHIBITION OF RENAL GLUCOSE REABSORPTION BY POTENT SGLT2 INHIBITORS**

Based on the modeling and simulations, it is likely that the apparently moderate inhibition of renal glucose reabsorption induced by potent SGLT2 inhibitors is a combined result of two physiological determinants:


Rieg et al. (2014) recently observed a 56% of lowering of renal glucose reabsorption in mice with complete SGLT2 inhibition, and an entire demolition of reabsorption in mice lacking both SGLT1 and SGLT2. This result confirms the hypothesis of SGLT1 compensation. The extrapolation of this finding to the clinic, however, is complicated by the differences in the experimental conditions in the Rieg et al. (2014) and clinical trials (Komoroski et al., 2009a,b; Devineni et al., 2013; Heise et al., 2013) (see Table S1). While a complete blockage of SGLT2 is likely in the Rieg et al. (2014) with drastically elevated concentration (free plasma concentration at least 10–15-fold higher than *in vitro* IC50) of empagliflozin over the duration of 30 min for UGE collection, in the clinical trials with once daily dosing, it is unlikely to maintain a 100% blockage of SGLT2 throughout a day over which luminal drug concentrations decline and 24 h UGE is collected. Thus, to explain the apparently moderate inhibition of renal glucose reabsorption by potent SGLT2 inhibitors in the clinic, the residual activity of SGLT2 should not be overlooked.

It is worth pointing out that dapagliflozin does severely suppress SGLT2 activity at its approved dose of 10 mg/day, as demonstrated by the simulations (Figure S2) at steady state in a hypothetical healthy subject with a constant plasma glucose level of 100 mg/dL treated with dapagliflozin at various doses. The SGLT2 activity decreases with increase in dose; from 20 mg onward, there is mild further decrease in SGLT2 activity. For SGLT1, its activity is nearly saturated at 10 mg. These results seem to be consistent with previous clinical observations that the UGE effect of dapagliflozin saturates at 20 mg (Komoroski et al., 2009a).

## **EFFECT OF AN SGLT1/2 DUAL INHIBITOR ON GLUCOSURIA IN COMPARISON WITH A SELECTIVE SGLT2 INHIBITOR**

It has been a question whether or not an SGLT1/2 dual inhibitor would induce greater glucosuria than a highly selective SGLT2 inhibitor (Chao and Henry, 2010; Abdul-Ghani et al., 2013). Abdul-Ghani et al. (2013) hypothesized that glucosuria induced by an SGLT2 inhibitor with a moderate selectivity over SGLT1 (e.g., capable of inhibiting SGLT1 activity by 30%) may be substantially greater than with a highly selective SGLT2 inhibitor. Using our model, we examined the sensitivity of UGE to Vmax1 and Ki1 in humans. We found that UGE was mildly to moderately sensitive to Vmax1 but not Ki1 in the presumably clinically tolerable ranges. The insensitivity to Ki1 is implied by Equation (2). With a treatment of 10 mg dapagliflozin, the glucose concentration in the PST rises to at least 20-fold of Km1. In order to moderately suppress SGLT1-mediated reabsorption through competitive inhibition, the luminal inhibitor exposure has to reach tens of fold of Ki1, a level that cannot be safely achieved in humans.

Therefore, without the consideration of its effect on intestinal SGLT1, whether or not a dual inhibitor will induce stronger glucosuria than a selective SGLT2 inhibitor is dependent on the mode of interaction between the dual inhibitor and SGLT1. A competitive inhibition of SGLT1 is unlikely to afford the dual inhibitor augmented effect on glucosuria. Other modes of inhibitions (non-competitive or uncompetitive) that attenuate Vmax1 may augment glucosuria mildly to moderately with a dual inhibitor.

In summary, to clarify mechanistically and quantitatively the operating characteristics of SGLT1 and SGLT2 in renal glucose reabsorption, we developed a systems pharmacology model with emphasis on renal glucose filtration, reabsorption, and transfer along the proximal tubules with or without SGLT1/2 inhibition. The model was calibrated using DeFronzo et al. (2013) and evaluated against three other data sets (Mogensen, 1971; Wolf et al., 2009; Polidori et al., 2013). Simulations using this model provided insights into the operating characteristics of SGLTs under normo- and hyperglycemic conditions in the healthy and diabetic state with or without SGLT2 inhibition. The simulations solidified the current concept of the relative contributions of SGLT1/2 to renal glucose reabsorption without the presence of SGLT2 inhibition. Moreover, the simulations elucidated quantitatively the operating characteristics of SGLTs when SGLT2 is inhibited. Further simulations clarified the relationships between SGLT1/2 capacity and renal glucose reabsorption in humans. Based on our modeling and simulations, we propose that the apparent moderate inhibition of renal glucose reabsorption observed clinically with SGLT2 inhibitors is a combined result of two physiological determinants, SGLT1 compensation and residual SGLT2 activity. This model will be valuable in mapping SGLT2 genotype to its functionality, and in predicting, through the incorporation of a plasma glucose-insulin model, the efficacy of an SGLT2 inhibitor in patients with diabetes, especially pediatric patients and patients with type 1 diabetes, for whom clinical data remain scarce.

## **AUTHOR CONTRIBUTIONS**

Yasong Lu: designed and executed the study, wrote and finalized the report; Steven C. Griffen: designed the study, critically reviewed and approved the report; David W. Boulton: designed the study, critically reviewed and approved the report; Tarek A. Leil: designed the study, reviewed the execution, critically reviewed and approved the report.

## **SUPPLEMENTARY MATERIAL**

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

### **REFERENCES**


malabsorption by trafficking defects. *Biochim. Biophys. Acta* 1453, 297–303. doi: 10.1016/S0925-4439(98)00109-4


**Conflict of Interest Statement:** This study was sponsored by Bristol-Myers Squibb and AstraZeneca. All authors were employees of Bristol-Myers Squibb when the study was being conducted.

*Received: 28 July 2014; accepted: 24 November 2014; published online: 10 December 2014.*

*Citation: Lu Y, Griffen SC, Boulton DW and Leil TA (2014) Use of systems pharmacology modeling to elucidate the operating characteristics of SGLT1 and SGLT2 in renal glucose reabsorption in humans. Front. Pharmacol. 5:274. doi: 10.3389/fphar. 2014.00274*

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

*Copyright © 2014 Lu, Griffen, Boulton and Leil. 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.*

## A computer-based quantitative systems pharmacology model of negative symptoms in schizophrenia: exploring glycine modulation of excitation-inhibition balance

## *Athan Spiros 1, Patrick Roberts 1,2 and Hugo Geerts 1,3\**

*<sup>1</sup> Computational Neuropharmacology, In Silico Biosciences, Inc., Berwyn, PA, USA*

*<sup>2</sup> Department of Biomedical Engineering, Oregon Health and Science University, Portland, OR, USA*

*<sup>3</sup> Department of Laboratory Pathology, Perelman School of Medicine, University of Pennsylvania, Philadelphia, PA, USA*

## *Edited by:*

*Tarek A. Leil, Bristol-Myers Squibb, USA*

#### *Reviewed by:*

*Ronald P. Hammer, University of Arizona College of Medicine, Phoenix, USA Barry Horwitz, National Institutes of Health, USA*

#### *\*Correspondence:*

*Hugo Geerts, Computational Neuropharmacology, In Silico Biosciences, Inc., 686 Westwind Dr., Berwyn, PA 19312, USA e-mail: hugo-geerts@ in-silico-biosciences.com*

Although many antipsychotics can reasonably control positive symptoms in schizophrenia, patients' return to society is often hindered by negative symptoms and cognitive deficits. As an alternative to animal rodent models that are often not very predictive for the clinical situation, we developed a new computer-based mechanistic modeling approach. This Quantitative Systems Pharmacology approach combines preclinical basic neurophysiology of a biophysically realistic neuronal ventromedial cortical-ventral striatal network identified from human imaging studies that are associated with negative symptoms. Calibration of a few biological coupling parameters using a retrospective clinical database of 34 drug-dose combinations resulted in correlation coefficients greater than 0.60, while a robust quantitative prediction of a number of independent trials was observed. We then simulated the effect of glycine modulation on the anticipated clinical outcomes. The quantitative biochemistry of glycine interaction with the different NMDA-NR2 subunits, neurodevelopmental trajectory of the NMDA-NR2B in the human schizophrenia pathology, their specific localization on excitatory vs. inhibitory interneurons and the electrogenic nature of the glycine transporter resulted in an inverse U-shape dose-response with an optimum in the low micromolar glycine concentration. Quantitative systems pharmacology based computer modeling of complex humanized brain circuits is a powerful alternative approach to explain the non-monotonic dose-response observed in past clinical trial outcomes with sarcosine, D-cycloserine, glycine, or D-serine or with glycine transporter inhibitors. In general it can be helpful to better understand the human neurophysiology of negative symptoms, especially with targets that show non-monotonic dose-responses.

**Keywords: schizophrenia, negative symptoms, glycine, computer model, dose-response relationship, drug**

## **INTRODUCTION**

Negative symptoms in schizophrenia are a major cause of functional deficit for patients wanting to return to professional life. While many of the approved antipsychotics can control the positive symptoms, negative symptom dysfunction is often not addressed properly by drug therapy alone (Rosenbaum et al., 2012). In addition, there are species differences for animal models that have large ramifications for drug development in schizophrenia (Peleg-Raibstein et al., 2012) and consequently psychiatric disorders have one of the lowest probabilities of clinical success, close to 7% (Hay et al., 2014). Because of these limitations, companies are de-emphasizing psychiatric diseases (Hyman, 2014), suggesting a need for completely novel technologies.

Negative symptoms can be divided in two moderately correlated factors (Horan et al., 2011): experiential impairments (diminished motivation and enjoyment of social, vocational, and recreational activities) and expressive impairments (diminished non-verbal and verbal communication). Experiential impairments are best represented by avolition and anhedonia, while expressive impairments are related to flat affect. Both these dimensions play an important role in the clinical phenotype.

Glutamate modulation through increased glycine mediated stimulation of the NMDA-R has been proposed as a strategy for addressing negative symptoms in schizophrenia. Consequently, a number of glycine modulators have been studied in clinical trials. In humans, the GlyT1 inhibitor ORG25935 reduced

**Abbreviations:** ASL-fMRI, Arterial Spin Labeling functional Magnetic Resonance Imaging; BOLD-fMRI, Blood Oxygen Level Dependent functional Magnetic Resonance Imaging; CBF, Cerebral Blood flow; DA, Dopamine; dlPFC, dorsolateral Prefrontal Cortex; EC50, concentration at which 50% of maximal effect is reached; Enk, encephalin; ER, endoplasmic reticulum; GlyT1, Glycine Transporter subtype 1; mAChR, muscarinic acetylcholine receptor; mGluR, metabotropic glutamate receptor; MSN, Medium Spiny Neuron (major cell type in striatum); NMDA, N-methyl-D-aspartate (glutamate receptor subtype); PANSS, Positive and Negative Symptoms Scale in schizophrenia; PDSP, Psycho-active Drug Screening Program; PET, Positron Emission Tomography; PFC, Prefrontal Cortex; PSD-95, postsynaptic density complex; QSP, Quantitative Systems Pharmacology; SP, Substance P; vmPFC, ventro-medial Prefrontal Cortex.

the ketamine-induced increases in measures of psychosis and perceptual alterations with an effect size of 0.71 and 0.98, respectively, but worsened some aspects of learning and delayed recall (D'Souza et al., 2012). Studies with the GlyT1 inhibitor GSK1018921 suggested that target engagements up to 80% were well tolerated (Ouellet et al., 2011). The Janssen GlyT1 inhibitor R213129 enhanced scopolamine-induced finger tapping impairment in healthy volunteers, while electroencephalography alpha power was increased and scopolamine-induced impairment of the Stroop test was partly reversed (Liem-Moolenaar et al., 2010). The Pfizer GlyT1 inhibitor PF03463275 was ineffective at the highest dose (NIH NCATS website http://www*.*ncats*.*nih*.* gov/research/reengineering/rescue-repurpose/therapeuticuses/ directory*.*html).

In a meta-analysis with 800 subjects from 26 studies, glycine, D-serine, and sarcosine had effect sizes ranging from 0.40 in negative symptoms to 0.28 for cognitive and 0.26 for positive symptoms, whereas D-cycloserine did not improve any symptom domain. Interestingly, patients on risperidone or olanzapine, but not clozapine, improved (Tsai et al., 2004).

Glycine directly activates the glycineB site on the NMDA-R, but needs to be given in large quantities; D-serine is another endogenous activator of the NMDA-R on a different binding site and sarcosine was found to be a GlyT1 inhibitor (Wolkenberg and Sur, 2010). The absence of target engagement data in these clinical trials makes it difficult to interpret the clinical outcome.

Preclinical data on cognitive effects together with target engagement studies in non-human primates for two GlyT1 inhibitors strongly suggest an inverse U-shape dose-response (Eddins et al., 2014); in this study the highest doses consistently failed to improve cognition and bitopertin was found to be effective only at the lowest and medium doses, but not at the highest dose. An inverse U-shape dose-response is a difficult property for any clinical trial development; although such a dose-response is often observed in preclinical animal models, it is often difficult to relate this to actual human target engagement levels. Therefore, exploring the neurophysiology of such complex dose-responses in a humanized translational model is of crucial importance to drug development. In this report we will use an *in silico* quantitative systems pharmacology model (Geerts et al., 2013b) that integrates preclinical information with clinical neuropathology, imaging, and clinical data and that has been successful for cognitive enhancements in schizophrenia (Geerts et al., 2013a) and Alzheimer's disease (Roberts et al., 2012; Nicholas et al., 2013) and for motor side-effects of new antipsychotics (Geerts et al., 2012).

The remainder of the introduction will be devoted to the biological rationale for identifying the brain regions and neurophysiological processes that play a role in the clinical phenotype of negative symptoms. Unlike preclinical animal models, we will use predominantly imaging studies from patients and their relationship to clinical scales.

## **BIOLOGICAL RATIONALE FOR COMPUTER MODEL OF NEGATIVE SYMPTOMS**

## *Brain regions/neurophysiology involved in negative symptoms*

The prefrontal cortex and ventral striatum are key brain regions involved in the processing of negative symptoms. From ASL-fMRI imaging studies to measure cerebral blood flow (CBF) in schizophrenic patients on antipsychotics medications (Pinkham et al., 2011), hypofrontality was most prominent in individuals with more severe negative symptoms. A large meta-analysis of 25 imaging studies (Goghari et al., 2010) suggests an inverse correlation between BOLD-fMRI activity of the ventromedial cortex and the degree of negative symptoms. Metabolic activity, measured by PET imaging, is reduced as negative symptoms increase in patients without antipsychotics (Wolkin et al., 1992) and physical anhedonia scale scores were negatively correlated with the hypoactive dorsomedial PFC metabolism (Park et al., 2009).

Another study suggests that activity of R. orbitofrontal cortex, but not anterior cingulate correlates with the self-reported Chapman Physical Anhedonia Scale (Harvey et al., 2010). As anhedonia together with avolition and apathy form the more "experiental" factor in negative symptoms, as opposed to flat affect that is more "expressive" (Horan et al., 2011); this suggests that lower activity of the R. orbitofrontal dysfunction might play a role in negative symptoms. Moreover, an inverse correlation of negative symptoms with R. anterior prefrontal cortex activity at rest (Mingoia et al., 2012) suggests that basal cortical activity is proportionally lower in patients with predominantly negative symptoms but the identity of the cortical region depends upon the task involved or the measurement condition.

This overview suggests that the cortical activity especially of the vmPFC and the right orbitofrontal cortex is lower in schizophrenia patients, and that increased activation might correspond to improved symptoms.

Imaging studies of ventral striatum pathology in schizophrenia (Menon et al., 2001; Harvey et al., 2010) suggest a profound and proportional dysfunction, with more negative symptoms associated with decreased activation level. In patients, lower ventral striatum activation in patients is proportional to the severity of negative symptoms, an effect that is independent of medication (whether medication-free, on typical or atypical antipsychotics) (Juckel et al., 2006a,b). In schizophrenia patients in psychotic remission (Sorg et al., 2013) basal activity of the ventral striatum is increased and this increase is correlated with improvements of negative symptoms such as emotional withdrawal and blunted affect.

## *Cellular localization of NMDA-NR2 subunits*

The activity of the cortical region is driven by pyramidal cell firing in general and by glutamatergic action in particular. Therefore, NMDA-R is an interesting target for negative symptoms. However, because the cortical activity is defined by the balance of excitation over inhibition, it is of interest to take into account the differential localization of NMDA-R on pyramidal cells and interneurons in cortical circuits. mRNA localization studies of different NMDA-NR2 subunits in the rat and mice hippocampus, suggest that NR2C*/*2D are localized on inhibitory interneurons while NR2A/NR2B seem to be more concentrated on pyramidal cells (Monyer et al., 1994). Functional evidence was provided by elimination of NR2C subunit having no effect on the strongly rectifying NMDA current in pyramidal cells of the prefrontal cortex (Zhang et al., 2012).

Reelin deficient heterozygous mice showed significantly enhanced MK-801-induced locomotor hyperactivity and startle, which was associated with significant up-regulation of NR1 subunits, but down-regulation of NR2C subunits in the frontal cortex (van den Buuse et al., 2012), suggesting that loss of activity on inhibitory neurons through reduced NMDA-NR2C synapses leads to a lower GABA tone, a functional disinhibition, and a higher locomotor activity. These and other preclinical data strongly suggests that while NR2A and NR2B are expressed on pyramidal excitatory cells, the NR2C subunit is localized on inhibitory neurons.

## *Change of NMDA subunits with schizophrenia pathology*

The NR2B subunit is upregulated during neurodevelopment of the brain and is likely to play a relatively larger role in schizophrenia, in line with the neurodevelopmental hypothesis of schizophrenia pathology. Indeed in postmortem dorso-lateral prefrontal cortex samples of schizophrenia patients vs. healthy control, increased phosphorylation of NR2B at Y1336 is found (Funk et al., 2012), probably leading to a higher functional activity by reducing endocytosis (Jiang et al., 2011; Li et al., 2011). In patients with schizophrenia, a significant effect of GRIN2B (human NMDA receptor 2B subunit gene, NR2B) genotype on habituation (Hokyo et al., 2010) suggests a bigger role for NR2B mediated processes.

Altered expression of mRNA for proteins involved in in microtubule-associated tracking complex of NR2B such as KIF17, APBA1, CASK, mLin7A, and mLin7C in cortical layers III and IV of schizophrenia patients, which overlapped with NR2B but not NR2A transcripts suggests that NR2Bcontaining NMDA receptor transport could be selectively compromised in schizophrenia (Kristiansen et al., 2010a,b). In a subcellular endoplasmic reticulum (ER)-enriched fraction from postmortem brain, ER expression of NR2B and PSD-95 was decreased in dorsolateral prefrontal cortex in schizophrenia. The data suggest that changes in NR2B processing in schizophrenia involve increased ER exit of NR2B containing NMDA receptors suggesting a higher membrane expression level (Kristiansen et al., 2010b).

Furthermore, a cross-sectional study of over 900 human brains from the publicly available BrainCloud website (http:// braincloud*.*jhmi*.*edu/) suggests an increase in cortical mRNA for the NR2B subunit during the adolescent period (10–20 years) that reverts for older brains. This suggests that during neurodevelopment the NR2B subunit is upregulated in the human brain but its expression tends to decrease with age.

In summary these data suggest that the NMDA- NR2B subunit is upregulated in schizophrenia patients.

## **GLYCINE TRANSPORTER PHYSIOLOGY**

In order to estimate the range of free glycine level that can be readily achieved in the living human brain, we need to consider the neurophysiology of the glycine transporter T1, found mostly on astrocytes but also on neuronal cells and is a co-transporter system driven by the electrogenic movement of 2 Na+ and 1 Cl− over the cell membrane at a slow turnover rate of 10/s (Cherubino et al., 2010). Kinetics follow Michaelis-Menten dynamics with Km in the range of 10–20 uM (Okamoto et al., 2009; Cherubino et al., 2010).

The astrocyte membrane potential is in the range of −75 mV (Ma et al., 2014) and does not share the same temporal dynamics as neuronal cells. The membrane can depolarize substantially in the case of ischemic and traumatic brain injury (Strong and Dardis, 2005), but we assume that the astrocyte membrane potential is close to the equilibirum value in schizophrenia.

## **METHODS**

## **RECEPTOR COMPETITION MODEL**

Many antipsychotic drugs on the market have different affinities for multiple receptors, therefore calculating the receptor change for a given exposure level of the drug at each of these receptors is important, because they will affect the membrane potential of key neuronal circuits and their emergent properties.

The receptor model simulates the competition between endogenous neurotransmitter and up to four agents, (for instance two drugs with their metabolites or a drug and radioactive tracer) at postsynaptic receptors with full presynaptic autoreceptor coupling to neurotransmitter release based on the affinities of the drug for all receptors in the synaptic cleft (Spiros et al., 2010). This is performed using a set of ordinary differential equations that takes into account different neurotransmitter release patterns and modulated by presynaptic autoreceptors, including presynaptic facilitation and depression processes. The dopaminergic synapse is further calibrated (Spiros et al., 2010) using data on dopamine dynamics measured with fast cyclic-voltammetry in monkey slices (Cragg et al., 2000) and human cortical imaging data (Slifstein et al., 2008), while the serotonin synapse with 5-HT1B as a presynaptic autoreceptor is calibrated using a combination of preclinical fast cyclic voltammetry constrained by human imaging data (Roberts et al., 2012).

The affinity parameters for each antipsychotic and neurotransmitter for human receptors were derived from the *in vitro* experiments performed at the Psychoactive Drug Screening Program (PDSP) and reported in the PDSP database (http://pdsp*.*med*.* unc*.*edu/indexR*.*html) with the advantage that the affinity values are derived under the same standardized assay conditions. For different values of target engagement (e.g., D2R occupancy), we then calculated the change of postsynaptic receptor activation for all the receptors involved in the computer model based on the affinities of the drug for different receptors.

## **CORTICAL-STRIATAL MODEL FOR NEGATIVE SYMPTOMS**

Based on the human imaging studies, we developed a dual cortical-striatal model for the neurobiology of negative symptoms (**Figure 1**). The cortical neuronal network consists of 20 excitatory neurons and 10 inhibitory interneurons and has been described before (Geerts et al., 2013a). This model has been calibrated from *in vivo* single-unit recordings in primates during a working memory task and reduces some of the problems associated with species difference in inhibitory tone. Synchronous firing of the target pyramidal cells is initiated by injecting a transient current at *t* = 2000 ms. The network then fires in a synchronized pattern before it gets degraded by the background noise and the interference of the distractor neurons. The simulated

neural activity represents the right orbitofrontal cortex or the vmPFC.

Functional representations, driven by preclinical experiments on the coupling between receptor activation and changes in voltage-gated ion channel conductance, of the dopamine (D1, D2, D3, D4, DAT, COMT), serotonin (5HT1A, 5HT1B, 5HT2A, 5HT2C, 5HT3, 5HT4, 5HT6, SERT), norepinephrine (alpha1A, alpha2A, NET), cholinergic (M1 mAChR, M2 mAChR, α<sup>7</sup> nAChR, α4β<sup>2</sup> nAChR, and AChE), glutamate NMDA (different subunits NR2A-NR2B-NR2C), AMPA, mGluR2, mGluR5, GlyT1, GABA-A α<sup>1</sup> and GABA-A α2, histamine H3 and PDE-10 targets are currently implemented in the model.

Although the intracellular pathways activated by receptor modulation are not modeled in full detail, we implement the effects as a transfer function on ion channel permeability or transporter functionality. For instance a change in dopamine D1R activation on cortical neurons is implemented by changing the slow K<sup>+</sup> channel *Iks* conductance (Yang and Seamans, 1996) and the High-voltage activated (Hva) Ca++ -channel, based on preclinical electrophysiological measurements (Law-Tho et al., 1994).

Schizophrenia pathology in the cortical network is introduced as a reduction in glutamate tone (Coyle, 2006), decreased dopamine tone (Meyer-Lindenberg et al., 2002; Weinberger, 2007) in the cortex, impaired GABA physiology through a decrease in GAD67 activity resulting in lower GABA release (Gonzalez-Burgos et al., 2010) and increased background noise level (Winterer et al., 2000). Such a pathology when implemented in the computer model leads to a deterioration of a marker for cognitive outcome of about 1.5 standard deviations (Geerts et al., 2013a). In this case, rather than the length of time a certain firing pattern can be independently held, the cortical readouts of the model for negative symptoms of schizophrenia is the average firing rate and BOLD-fMRI.

The ventral striatum model has been described in detail as part of the quantitative systems pharmacology platform for schizophrenia (Geerts et al., 2012; Spiros et al., 2012). Briefly, the model calculates the excitability of the medium spiny neuron (MSN), the major GABA-ergic cell type in the nucleus accumbens, when driven by afferent cortical projections and gated by both hippocampal and amygdala projections. Changes in membrane potential are calculated using partial differential equations that are solved in NEURON (Hines and Carnevale, 1997). If C is the membrane capacitance, then the time course of the membrane potential V can be determined from the following equation:

$$C\frac{\partial V}{\partial t} = I\_{\rm Ksi} + I\_{\rm KA} + I\_{\rm CI} + I\_{\rm Ca} + \dots \tag{1}$$

where *IX* is the current associated with channel X.

We simulate three types of neuronal MSN cells: SP<sup>+</sup> = D1R<sup>+</sup> cells that project to the direct pathway; Enk<sup>+</sup> = D2R<sup>+</sup> cells that project to the indirect pathway; and a small number of D+ <sup>1</sup> D<sup>+</sup> 2 cells that project to both pathways. In the SP<sup>+</sup> cells the D1R mostly affects the *Kir*<sup>2</sup> channel and increases the L-type Ca++ current (Hernandez-Lopez et al., 2000), while in Enk<sup>+</sup> cells, D2R activation affects the A-type K+ current (Falk et al., 2006). In addition, D2R activity modulates the presynaptic Glu release on the afferent cortical fibers (O'Donnell and Grace, 1994; Bamford et al., 2004).

For instance, the inward rectifying potassium current, *Kir*2, is modified by the dopamine D1R activation u (Kuzhikandathil and Oxford, 2002; Falk et al., 2008) so that the total current, *I* = *u* · *IKir*2. With a conductance, *gK*, and a reversal potential, *EK* = −90 mV, the current is given by *IKir*<sup>2</sup> = *gK (V* − *EK)* with a voltage dependent form

$$\log \epsilon = \bar{\mathbf{g}}\_{\mathbf{K}} \frac{1}{1 + \exp\left(-\frac{V - V\_h^K}{V\_\epsilon^K}\right)}\tag{2}$$

where *<sup>g</sup>*¯*<sup>K</sup>* <sup>=</sup> <sup>1</sup>*.*2 mS/cm<sup>2</sup> is the maximum conductance, *Vh* = −111 mV is the value of the membrane potential that causes half activation and *Vc* = −11 mV describes the sensitivity of the change (Mermelstein et al., 1998; Gruber et al., 2003).

The amount of DA released in the striatal dopaminergic synapse is increased by 5-HT2C receptor inhibition (Abdallah et al., 2009), while 5-HT3R antagonism decreases striatal DA (De Deurwaerdere et al., 1998; Porras et al., 2003). Cholinergic modulation affects the excitability of MSN through an effect on Cl- channel (Shen et al., 2005, 2007) through postsynaptic M1R mAChR. In addition, M2 mAChR located on corticostriatal terminals (Hersch et al., 1994) inhibit the glutamatergic input to MSNs (Malenka and Kocsis, 1988; Sugita et al., 1991; Calabresi et al., 1998; Hernandez-Echeagaray et al., 1998). Adrenergic alpha1A-R block decreases gating signal stimulation of the GABA spiny neuron (Braga et al., 2004; Aroniadou-Anderjaska et al., 2007). All these processes are implemented using the appropriate differential equations with a linear relationship between the increase of DA and normalized activation level.

#### **IMPLEMENTATION OF THE BOLD-fMRI READOUT**

In order to calculate a measure of the BOLD-fMRI outcome from the computer model, we implemented a series of biophysical relations between excitatory and inhibitory neuronal activity as determined by experimental studies (Sotero and Trujillo-Barreto, 2007, 2008). The relevant equations are implemented describing the relationships between excitatory and inhibitory neuronal activity, glucose consumed, oxygen consumed, and CBF changes to obtain a measure of the BOLD signal with the Balloon model (Buxton et al., 2004; Buxton, 2012) with the parameters provided from a review study (Sotero and Trujillo-Barreto, 2008).

With *v*(*t*) the normalized cerebral blood volume, *f*(*t*) the normalized CBF and *q*(*t*) the doxyhemoglobin content, the BOLDfMRI signal *y*(*t*) is described by

$$\mathcal{Y}(t) = V0(a\_1\left(1 - q\right) - a\_2(1 - \nu) \tag{3}$$

$$\text{With}\frac{d\nu(t)}{dt} = \frac{1}{t\_0} (f\left(t\right) - f\_{\text{out}}\left(\nu, t\right))\tag{4}$$

$$\frac{dq\,(t)}{dt} = \frac{1}{t\_0}(m\,(t) - f\_{\rm out}\,(\nu, t)\,\frac{q\,(t)}{\nu\,(t)})\tag{5}$$

$$fout\left(\nu, t\right) = \nu \exp\left(\frac{1}{a}\right) + \pi d\nu(t)/dt\tag{6}$$

Furthermore, with *me*(*t*) and *mi*(*t*) the metabolic rate of oxygen consumption from excitatory and inhibitory cells, respectively

$$m\left(t\right) = \left(\gamma m\_{\epsilon}\left(t\right) + m\_{\dot{\imath}}\left(t\right)\right) / \left(\gamma + 1\right);\tag{7}$$

$$m\_i(t) = \mathcal{g}\_i(t) \text{ and } m\_\varepsilon(t) = \mathcal{g}\_\varepsilon(t)(2 - \ge (t)\,)/(2 - \ge 0 \,) \tag{8}$$

$$\text{And}\,\mathbf{g}\,\left(\mathbf{t}\right) = \left(2\gamma\mathbf{g}\_e\left(\mathbf{t}\right) + \left(2 - \mathbf{x}o\right)\mathbf{g}\_i\left(\mathbf{t}\right)\right)\left(\left(2\gamma + 2 - \mathbf{x}o\right)\right) \quad \text{(9)}$$

With *ue*(*t*) and *ui*(*t*) the excitatory and inhibitory neuronal activity, both *ge*(*t*) and *gi*(*t*), the glucose level normalized to baseline consumption, are further defined by

$$\frac{d\mathbf{g}\_{\mathfrak{e}}(t)}{dt} = \mathbf{s}\_{\mathfrak{e}}(t) \tag{10}$$

$$\text{and } \frac{d\mathbf{g}\_i(t)}{dt} = s\_i(t) \text{ with}$$

$$\frac{ds\_{\varepsilon}\left(t\right)}{dt} = \frac{a\_{\varepsilon}\left(u\_{\varepsilon}\left(t-\delta\_{\varepsilon}\right)-1\right)}{\mathfrak{r}\_{\varepsilon}} - \frac{2s\_{\varepsilon}\left(t\right)}{\mathfrak{r}\_{\varepsilon}} - \frac{\mathfrak{g}\_{\varepsilon}\left(t\right)-1}{\mathfrak{r}\_{\varepsilon}\*\mathfrak{r}\_{\varepsilon}} \tag{11}$$

With an identical equation for *si*(*t*) with all indexes referring to inhibitory interneuron activity *ui*(*t*).

The CBF f(t) is defined by *df(t) dt* = *s*(*t*).

Where *ds(t) dt* = *ε ee t* − *δ<sup>f</sup>* − 1 <sup>−</sup> *<sup>s</sup>(t) <sup>τ</sup><sup>s</sup>* − (*f (t)* − 1*/τ<sup>f</sup>* .

Values for the different constants are given in Table 3 of Sotero and Trujillo-Barreto (2008). For instance, *ae* is the efficacy of glucose consumption response to excitation (1.2); *c* the steepness of the sigmoid function (2.5) and *d* the position of the threshold for the sigmoid function (1.6), τ the time constant that controls how fast Cerebral Blood volume adjusts to changes in CBF (10 s), *a*<sup>1</sup> weight for deoxyhemoglobin change (3.4) and a2 weight for blood volume change (Rosenbaum et al., 2012).

## **IMPLEMENTATION OF THE GLYCINE NEUROPHYSIOLOGY IN THE MODEL**

The ratio of NR2A/NR2B subunit on pyramidal excitatory synapses vs. NR2C/NR2D subunit on inhibitory cells is an important driver of glycine modulation. Glycine interacts differently with different NMDA-NR2 subtypes even if the binding site is on the NR1 subunit. The potentiation of NMDA current by glycine has been measured experimentally and can be described by a Hill equation. The NMDA-R conductance g can be described as follows.

$$\mathbf{g} = \mathbf{g}\_{\text{max}} \frac{[\mathbf{G} \mathbf{l} \mathbf{y}]^n}{[\mathbf{G} \mathbf{l} \mathbf{y}]^n + [\mathbf{E} \mathbf{C}\_{50}]^n} \tag{12}$$

With *g*max a maximal conductance value, [*Gly*] the concentration of extracellular glycine and with numerical values for EC50 and Hill slope determined from a number of experimental studies (Kutsuwada et al., 1992; Laurie and Seeburg, 1994; Matsui et al., 1995; Woodward et al., 1995; Chen et al., 2008).

The EC50 (concentration at which effect is 50% of maximum) and Hill slopes (*n*) for different experimental conditions are given in **Table 1**. It will be clear that with different values for the Hill equation and EC50 we get rich dynamics in terms of the ratio of the NMDA currents on pyramidal-pyramidal synapses (mostly NR2A-NR2B) vs. the NMDA currents on pyramidal-inhibitory synapses driven by the NR2C-NR2D subunits (see **Figure 3** for example).

**Table 1 | Experimentally determined values for EC50 and Hill slopes for glycine-glycine site interaction on the NMDA-NR2 subunit from different experimental conditions.**


*EC50 values are definitely higher and there is a trend for higher Hill slope for excitatory-excitatory NR2A/<sup>B</sup> subunits (especially when considering the NR2B subunits) over the excitatory-inhibitory NR2C/<sup>D</sup> subunits.*

## **CALIBRATION OF THE MODEL WITH CLINICAL DATA ON NEGATIVE SYMPTOMS**

The model is subsequently calibrated using historical clinical trials. Historical clinical data in schizophrenia patients were collected by querying PubMed with the keywords "drug X" and "schizophrenia" and trial" in the period since 1986. Restricting the data to clinical double-blind placebo-controlled studies on drug monotherapy using stable schizophrenia patients for a short duration (4–12 weeks), resulted in 91 papers and 71 drug-dose combinations. For each drug dose combination, we calculated the change in postsynaptic receptor activation using the receptor competition model using the appropriate affinities of the neurotransmitter, the drug and its metabolite. "

We assume a linear normalized relationship between receptor activation and biological effect on physiological responses such as *Xeff <sup>Y</sup>* <sup>=</sup> *<sup>X</sup><sup>A</sup> <sup>Y</sup>* <sup>−</sup>*X<sup>C</sup> Y XC Y* ; (Equation 13) where *XA <sup>Y</sup>* and *<sup>X</sup><sup>C</sup> <sup>Y</sup>* are the actual activation levels of receptor X subtype Y (for instance D1) after treatment (A) and the untreated (placebo) control levels (C).

Such short clinical studies are common in the clinical testing of antipsychotics and motor side effects can arise very early with treatment. For each study, the average outcome of a patient group on the reported clinical trial was entered into a database. In the case of multiply reported results for the same drug-dose, the weighted average outcome based on number of patients was calculated. The list of clinical studies can be requested from the corresponding author.

## **RESULTS**

## **CALIBRATION OF THE MODEL FOR NEGATIVE SYMPTOMS OUTCOME**

Extraction of the relevant information from the clinical database results in 34 drug-dose combinations of short-term clinical trials

BOLDfMRI readout of the cortical part of the computer model and activity in the ventral striatum computer model. A weighted sum of these two parameters (40% cortical input and 60% ventral input) was defined as a proxy for negative symptoms. In addition, a few biological coupling parameters in the ventral striatum model were adjusted to achieve a robust correlation between model results and clinically reported outcomes.

**FIGURE 3 | (A)** Relative currents through NMDA- NR2C*/*<sup>D</sup> subunits in pyramidal-inhibitory glutamatergic synapses (e-i) over NMDA currents and NMDA- NR2A*/*<sup>B</sup> subunits in pyramidal-pyramidal glutamatergic synapses (e-e) and the ratio of NR2C*/*<sup>D</sup> over NR2A*/*<sup>B</sup> currents as a function of glycine concentration. This particular graph corresponds to a value of 0.35 and 0.25 uM for EC50 of NR2A*/*<sup>B</sup> and NR2C*/*<sup>D</sup> respectively and a Hill coefficient of 1.84 and 1.1 for NR2A*/*<sup>B</sup> and NR2C*/*<sup>D</sup> respectively. This particular case shows a minimum of 0.93 on the ratio of NR2*<sup>C</sup>/<sup>D</sup>* over NR2*<sup>A</sup>/<sup>B</sup>* at a glycine concentration of 1.1 uM. For certain combinations of the EC50 and Hill coefficients of the interaction between glycine and its binding site on NMDA receptors, the ratio shows a U-shape dose-response. For glycine concentrations where the e-i over e-e ratio drops below 1 (in this case between 0.5 and 30 uM), this leads to a relative dominance of the excitatory tone over the inhibitory tone and a greater firing and BOLD-fMRI outcome. **(B)** Dose-response of the model outcome for negative symptoms as a function of glycine concentration for different combinations of parameters for glycine-NR2C*/*<sup>D</sup> vs. glycine-NR2A*/*<sup>B</sup> interactions after calibration with clinical data. The parameter settings are EC50 for NR2C*/*D, Hill coefficient for NR2C*/*D, EC50 for NR2A*/*B, Hill coefficient for NR2A*/*B, respectively so that the conditions shown are A (0.2 uM, 0.7, 0.3 uM, 1.7); B (0.2 uM, 1, 0.22 uM, 1.5); C (0.21 uM, 1, 0.25 uM, 2); D (0.2 uM, 1, 0.35 uM, 2); and E (0.20 uM, 1.46, 0.35 uM, 1.78). In a clinical trial, glycine concentration can be altered using either a glycine transporter inhibitor or substitution with glycine. Note the inverse U-shape dose-response that parallels the ratio of current through the NMDA-NR2 subunits located at the excitatory-excitatory synapses over the current through the NMDA-NR2 subunits located at the excitatory-inhibitory synapses. Maximal effect is in the two-point range but is only obtainable if the typical glycine concentration is around 1 uM. With higher glycine concentration, in many cases the effect decreases, with a worsening of negative symptoms in some cases for glycine concentrations greater than 5 uM.

(4–12 weeks) with outcomes on the PANSS negative subscale. Each of these drug-dose combinations was then first simulated with the appropriate PET tracer displacement (with 11Craclopride) studies in the dopamine receptor competition model to yield the functional intrasynaptic drug concentration that corresponded to the observed tracer displacement. These functional drug doses were then applied to the receptor competition model for all postsynaptic receptors (other dopamine receptors, in addition to serotonergic, cholinergic, and norepinephrine receptors) in the computer model using the appropriate pharmacological interaction between that particular drug and the human receptor subtype. With the resulting changes in postsynaptic receptor activation for all synapses in the negative symptoms QSP computer model, *in silico* computer model results were obtained for each drug-dose combination, in particular changes in BOLDfMRI readout of the cortical part of the computer model and activity in the ventral striatum computer model. A weighted combination of these two parameters was then defined as a proxy for negative symptoms.

Those results were then compared to the corresponding actual clinical readouts on PANSS Negative. **Figure 2** shows the correlation between the model outcome of these drug-dose combinations and the actual reported changes in PANSS Negative. A few biological coupling parameters in the ventral striatum model were adjusted to achieve a robust correlation. The observed correlation suggests that the model captures a substantial part of the variance.

## **DOSE-RESPONSE OF GLYCINE MODULATION**

In order to simulate the effect of modulating glycine concentration on the network outcome, we proceeded by calculating the Hill equations for the interaction between glycine and the NR2 subunit on the excitatory-excitatory (e-e) glutamatergic synapses with NMDA subunits (NR2A-NR2B) and the excitatory-inhibitory (e-i) glutamatergic synapses with NMDA subunits NR2C-NR2D. **Figure 3A** shows the currents through the two synapse types and their ratios as a function of glycine concentration. Due to the complex non-linear effects, for a number of EC50 and Hill coefficient parameter settings, the ratio of inhibitory to excitatory NMDA response decreases, reaches a minimum and then increases again, allowing for a complex non-monotonic doseresponse. The range of EC50 and Hill coefficient parameters for which this can be observed is examined in the following section.

When entering these changes for the effect of glycine on ee and e-i glutamatergic synapses as corresponding changes in NMDA maximum conductances on the respective synapses, the impact on the output of the network can be calculated. **Figure 3B** shows the anticipated clinical outcome on the PANSS negative scale based on the calibration (**Figure 2**) for a number of combinations of the interaction parameters between glycine and its binding site on the NMDA receptor. The clinical benefit follows an inverse U-shape dose-response relationship with the free glycine concentration provided that the typical glycine concentrations are less than 0.5 uM, very similar to the dose-response of the ratio of the currents, NR2C*/*2D to NR2A*/*2B NMDA subtypes. In fact using a number of different values for the EC50 and Hill coefficient parameters, we can confirm the correspondence between these two dose-responses (data not shown).

Although the dose-responses in general are similar for the different glycine interaction parameters, the glycine concentration for which there is a clinical benefit varies substantially. As no experimental data are available for the interaction of glycine with the glycine binding site in the living human brain, in principle accordance between the model outcome and the clinical results will help to narrow down the range of likely parameters. In the following section we will address these issues in more detail.

## **SENSITIVITY OF THE MODEL TO BIOLOGICAL COUPLING AND PHYSIOLOGY PARAMETERS**

We then studied the sensitivity of the model outcome as a function of small changes on all possible parameters, including the changes that reflect the implementation of the schizophrenia pathology. This could be envisaged as the intrinsic variability in a large set of patients possibly driven by the appearance of genotypic changes in neurophysiological and neuropathological pathways. For this we allowed the parameter to change by a certain fraction around their calibrated value. A general effect size is calculated by dividing the difference between the maximum and minimum relative effects by the fractional range in parameter settings. From the outcomes reported in **Table 2**, it is clear that the changes in pathology implementation at the level of the cortical network lead to the largest changes in the results. The effect of the parameters on the Hill equations for the interaction between glycine and its binding site on the NMDA-R will be explored in detail in the following section.

## **SENSITIVITY ANALYSIS ON GLYCINE HILL EQUATION**

A crucial set of parameters is the relative values of the EC50 and Hill coefficients for the glycine-NMDA current effect through the

**Table 2 | Sensitivity of the model outcome for different parameters and coupling constants that are changed between 20 and 50% in both directions for the glycine dose-response.**


*Maximal and minimal effect on the dose-response are derived. The general effect size is calculated as the difference between maximal and minimal outcome divided by the input range. The data suggest that GABA and NMDA reduction have a big impact on the outcome, likely because they are driving the baseline excitation-inhibition balance upon which glycine modulation will act.* different NR2 subunits. As shown in the previous section, the inverse U-shape dose-response in the cortical network outcome corresponding to a clinical benefit on PANSS negative is associated with a U-shape dose-response of the ratio of inhibitory over excitatory effect (compare **Figure 3A** with **Figure 3B**). Note that the glycine concentration for maximal effect on the network corresponds to the glycine concentration of the minimum in the ratio of NR2C*/*<sup>D</sup> over NR2A*/*B.

Because the experimentally determined values reported in **Table 1** for the interaction of glycine with the human NMDA NR2 subunits were performed in an artificial *in vitro* system and could be quite different from the actual human *in vivo* situation, we systematically studied the effect of changing the values for EC50 on the network outcome.

**Figure 4** shows the sensitivity analysis when probing different parameter ranges for EC50 and Hill coefficients of the interaction between glycine and the NMDA-NR2 subunits. It is clear that higher Hill coefficients for the NR2A*/*<sup>B</sup> subunits compared to the NR2C*/*<sup>D</sup> subunits is necessary for a beneficial effect on the network outcome with the range increasing with larger differences between the Hill coefficients for NR2A*/*2B vs. NR2C*/*2D subunits.

## **PHYSIOLOGICAL RANGE OF GLYCINE CONCENTRATION**

The free glycine concentration in the human brain is regulated by a 2Na+-Cl−-Gly co-transporter system and its value is constrained by the Nernst-Goldman equation. Therefore, the functional free glycine concentrations in steady-state equilibrium conditions is dependent upon the range of concentrations for intracellular and extracellular Na+ and Cl−. Assuming the glycine transporter is located on astrocytes, all calculations are done for astrocyte membrane potential and intracellular ion concentrations. Astrocyte membrane potential, while not changing on the same time scale as neuronal membrane potential is supposed to be in the −50 to −70 mV range for steady state conditions.

With the exception of extreme pathological situations such as in stroke or neurotrauma, Na+ and Cl− ion concentrations in the human brain are tightly regulated. **Figure 5** shows a number of solutions to the Nernst-Goldman equation for different ranges of intracellular Na+ and intracellular glycine concentrations. All figures are derived for constant values of intracellular Cl− of 6 mM and extracellular Cl− concentration of 120 mM. It is clear that in the absence of extreme pathology, the range of free extracellular glycine is limited and is unlikely to exceed 10 uM.

### **INTERNALIZATION OF NMDA-R AT HIGH GLYCINE CONCENTRATIONS**

High glycine exposure in principle could overstimulate the NMDA-R and lead to epileptic seizures. In preclinical slice work, NMDA-R internalization has been observed at very high glycine concentrations, typically with an EC50 value in the range of 40 uM (Nong et al., 2003). While the Nernst-Goldman equations that regulate the free glycine concentration as a function of free Na and Cl strongly suggest that glycine concentrations beyond 5 uM are highly unlikely (see **Figure 5**), we nevertheless simulated the effect of NMDA-NR2 subunit internalization on the computer model outcome.

**FIGURE 4 | Sensitivity analysis of the inverse U-shape dose-response for glycine-dependent outcome of the computer model as a function of the Hill slope and the EC50 values for both classes of NMDA-NR2 subunits.** The x-axis represents EC50 values for the e-i NR2C*/*<sup>D</sup> subunits and the y-axis represents the EC50 values for the NR2A*/*<sup>B</sup> subunits. Shown are ratios of different Hill coefficients on NR2A*/*Bvs. NR2C*/*<sup>D</sup> subunits: **(A)** 1:1, **(B)** 1.25:1, **(C)** 1.5:1, **(D)** 1.75:1, and **(E)** 2:1. Combinations of values in green color are associated with a U-shape dose-response for NMDA effects of e-i over e-e with a minimum between 0.5 and 10 uM of glycine

concentration and therefore a beneficial effect on the network outcome. Combinations in yellow are associated with values resulting in monotonic decrease of the ratio and combinations in blue are associated with a non-biological minimum much smaller than 0.5 uM. The results indicate that a beneficial inverse U-shape dose-response on negative symptoms is associated with large differences between the Hill coefficients for NR2A*/*2B vs. NR2C*/*2D subunits. In addition, in most cases, a beneficial inverse U-shape dose-response is associated with greater EC50 values for NR2A*/*2B than for NR2C*/*2D.

We simulated two conditions of internalization (**Figure 6**). The first condition assumes the same EC50 value (40 uM) for both NR2A/B subunits as for NR2C/D subunits. As expected, for very low levels of glycine where the internalization process has a very limited effect, the dose-response shows a peak for glycine levels in the low uM range. The results further suggest that internalization of the NMDA-R assuming the same EC50 values for the two types of NR2 subunits leads to a substantial decrease (corresponding

to a collapse of network activity at glycine levels beyond 20 uM). The second assumption assumes that the internalization process would have the same affinity as glycine itself, i.e., an EC50 of 32 uM for the NR2A*/*<sup>B</sup> subunits and an EC50 of 50 uM for the NR2C*/*<sup>D</sup> subunits. In that case, the simulations suggest that the model outcome first improves beyond the no internalization case before collapsing at glycine levels beyond 40 uM.

## **DISCUSSION**

This report describes a quantitative systems pharmacology computer model based on physiologically realistic interactions in models of a cortical network and the ventral striatum. The major result of this simulation is the prediction of an inverse U-shape dose-response with glycine that is a consequence of the shifting balance between excitation and inhibition in the cortical network, secondary to an interesting difference in pharmacological properties of glycine for the different NMDA subunits regulating excitatory and inhibitory tone.

The sensitivity analysis suggests that there are a substantial number of parameter combinations that result in such an inverse

U-shape dose-response with glycine. The exact values for the interaction in the human brain is unknown and probably is different for each subject, but in general an inverse U-shape dose-response can be achieved when both the EC50 and the Hill coefficient for the glycine effect on the NR2*C/<sup>D</sup>* subunit is lower than for the NR2A*/*<sup>B</sup> subunit. It is of interest to note from different experimental data that on average the interaction of glycine with the NMDA-NR2B subunit indeed suggests a higher value for the Hill coefficient. As noted in the Biological Introduction, the neurodevelopmental trajectory of schizophrenia tends to delay the appearance of the mature NMDA-NR2A subunit, so that there is a relatively higher contribution of the NMDA-NR2B subunit to the excitatory tone in schizophrenics. This allows the interaction of glycine with the e-e NMDA receptor to have a higher Hill slope in combination with a lower EC50 concentration, promoting an inverse U-shape dose-response.

The clinical Phase II data with the glycine inhibitor bitopertin suggest indeed that the clinical outcome follows such as doseresponse. Some studies with D-serine, D-sarcosine or glycine have often reported mixed effects (Singh and Singh, 2011), with some but not all showing a clinical benefit and the interpretation is hampered by the lack of data on proper target engagement. It is conceivable that this could be a consequence of the nonlinear dose-response with patients on different points of the dose-response. It is also of interest to note that negative symptoms seem to be most improved when glycine or D-serine levels are increased (Singh and Singh, 2011). The QSP platform when calibrated suggests a rather limited effect of glycine modulation on clinical changes in PANSS negative in the range of 2–2.5 points. Note that the patient population used for calibration was not selected for extremely high negative symptomatology, resulting in a baseline PANSS negative between 18 and 24. Although this observed change is in line with published data on glutamatergic modulation (Singh and Singh, 2011), it suggests that such a limited effect might be difficult to be detected in clinical situations.

These results also lead to the important observation that because of the different contribution of the NMDA-NR2 subunits, the interaction between glycine and the NMDA receptors likely is different in healthy volunteers as compared with schizophrenia patients. Such a difference is extremely difficult to achieve in preclinical animal models. Therefore, great care needs to be taken to extrapolate positive or negative findings from a Phase I proof-of-concept study in healthy volunteers to actual schizophrenia patients.

Similarly, rodent models often lack the right mix of NMDA receptor subtypes to simulate very well the actual pathology mediated and often do not show the inverse U-shape doseresponse (Alberati et al., 2011). Not realizing this could substantially hamper the clinical development and often can lead to failed clinical trials. In addition, clinical trials have been performed as augmentation strategy, i.e., the glycine modulator intervention is given to patients on stable antipsychotic medication. Such comedications can have a direct effect through affecting the metabolism of the active compound which is dependent upon the genotype of the specific Cytochrome P450 enzyme but can also be modified by other comedications such as nicotine (Tsuda et al., 2014). The comedications can also have an indirect effect on the dose-response of glycine level modulation through non-linear interactions on the excitation-inhibition balance that affect the emergent properties such as the BOLD-fMRI signal. This paper does not address the issue of comedication, but we are planning to perform such an in-depth analysis in a follow-up paper.

## **LIMITATIONS OF THE MODEL**

Firstly, different reports suggest that D-serine plays a more important role as co-agonist on the NMDA-R in the cortex (Fossat et al., 2012) while other studies suggest a role for both D-serine and glycine in regulating neuronal morphology in rodent somatosensory cortex (Balu et al., 2012). To a certain extent, the interaction of D-serine with the co-agonist site on the NMDA-R is quite similar to glycine's interaction (Chen et al., 2008), but the level of free D-serine is regulated by serine racemase localized in neurons (Balu et al., 2014) and by a Na+-independent alanine-cysteineserine transporter system (Maucler et al., 2013). This suggests that most of the conclusions for glycine, with the exception of the limited range of glycine driven by its specific Na+-dependent co-transporter system can be applied to the modulation of D-serine.

There is also some discussion about the nature of the NMDA NR2 subunits on the inhibitory cell types in cortical networks. mRNA studies in the human brain localize the NR2C subunit predominantly to the cerebellum (Monyer et al., 1992) although there are lower levels present in the cortex (Allen Brain institute data http://human*.*brain-map*.*org/). However, NR2D subunits are likely present in cortical areas on excitatory-inhibitory synapses and could play a predominant role in the generation of the inhibitory tone. The exact interaction parameters between glycine and its co-agonist binding site on the NMDA receptor in the human brain are unknown, but the computer model suggests a range of interaction parameters that would correspond to an inverse U-shape dose-response. Although adding a greater contribution of NR2D to the inhibitory tone will increase the EC50 of the e-i interaction, it will also decrease the Hill coefficient as compared to the interaction of glycine with the e-e synapses, which has been shown to be favorable for an inverse U-shape dose-response.

The model presented here does not address the other modulatory agents such as extracellular proteins, zinc, polyamines, and neurosteroids. All these molecules can influence the dynamics of glycine-mediated amplification of NMDA-currents and we assumed that these modulators do not change in schizophrenia. In principle, if new data become available suggesting a change in these neuromodulators as a consequence of schizophrenia pathology, detailed biochemical data could be incorporated in this platform to estimate their impact.

In summary, this report simulates the anticipated doseresponse of glycine level modulation on an emergent property (BOLD-fMRI) of a computer-based neuronal circuit that has been calibrated against clinical outcomes for negative symptoms. It provides a physiological explanation for the appearance of an inverse U-shape dose-response based on a biologically constrained set of interaction parameters between glycine and the co-agonist site on different types of NMDA-NR2 subunits and the electrogenic character of the Gly/2Na/Cl co-transport system. A notable limitation is that this study deals with the effect of glycine modulation in the absence of any antipsychotic medication and therefore does not reflect the real clinical study design.

## **REFERENCES**


phase-locking during information processing. *Clin. Neurophysiol.* 111, 837–849. doi: 10.1016/S1388-2457(99)00322-3


**Conflict of Interest Statement:** All authors are employees of In Silico Biosciences. Hugo Geerts and Athan Spiros are the inventors on two patents related to Quantitative Systems Pharmacology.

*Received: 10 June 2014; accepted: 23 September 2014; published online: 21 October 2014.*

*Citation: Spiros A, Roberts P and Geerts H (2014) A computer-based quantitative systems pharmacology model of negative symptoms in schizophrenia: exploring glycine modulation of excitation-inhibition balance. Front. Pharmacol. 5:229. doi: 10.3389/ fphar.2014.00229*

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

*Copyright © 2014 Spiros, Roberts and Geerts. 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.*

## Exploring BSEP inhibition-mediated toxicity with a mechanistic model of drug-induced liver injury

#### *Jeffrey L. Woodhead1 \*, Kyunghee Yang1, Scott Q. Siler 1, Paul B. Watkins 1, Kim L. R. Brouwer 2, Hugh A. Barton3 and Brett A. Howell <sup>1</sup>*

*<sup>1</sup> The Hamner-UNC Institute for Drug Safety Sciences, The Hamner Institutes for Health Sciences, Research Triangle Park, NC, USA*

*<sup>2</sup> Division of Pharmacotherapy and Experimental Therapeutics, UNC-Eshelman School of Pharmacy, University of North Carolina at Chapel Hill, Chapel Hill, NC, USA <sup>3</sup> Pharmacokinetics, Dynamics, and Metabolism, Worldwide Research and Development, Pfizer, Inc. Groton, CT, USA*

#### *Edited by:*

*Sergey Ermakov, Bristol-Myers Squibb, USA*

#### *Reviewed by:*

*Adriana Maggi, University of Milan, Italy Rink-Jan Lohman, The University of*

*Queensland, Australia*

#### *\*Correspondence:*

*Jeffrey L. Woodhead, The Hamner-UNC Institute for Drug Safety Sciences, The Hamner Institutes for Health Sciences, 6 Davis Drive, Research Triangle Park, NC, USA*

*e-mail: jwoodhead@thehamner.org*

Inhibition of the bile salt export pump (BSEP) has been linked to incidence of drug-induced liver injury (DILI), presumably by the accumulation of toxic bile acids in the liver. We have previously constructed and validated a model of bile acid disposition within DILIsym®, a mechanistic model of DILI. In this paper, we use DILIsym® to simulate the DILI response of the hepatotoxic BSEP inhibitors bosentan and CP-724,714 and the non-hepatotoxic BSEP inhibitor telmisartan in humans in order to explore whether we can predict that hepatotoxic BSEP inhibitors can cause bile acid accumulation to reach toxic levels. We also simulate bosentan in rats in order to illuminate potential reasons behind the lack of toxicity in rats compared to the toxicity observed in humans. DILIsym® predicts that bosentan, but not telmisartan, will cause mild hepatocellular ATP decline and serum ALT elevation in a simulated population of humans. The difference in hepatotoxic potential between bosentan and telmisartan is consistent with clinical observations. However, DILIsym® underpredicts the incidence of bosentan toxicity. DILIsym® also predicts that bosentan will not cause toxicity in a simulated population of rats, and that the difference between the response to bosentan in rats and in humans is primarily due to the less toxic bile acid pool in rats. Our simulations also suggest a potential synergistic role for bile acid accumulation and mitochondrial electron transport chain (ETC) inhibition in producing the observed toxicity in CP-724,714, and suggest that CP-724,714 metabolites may also play a role in the observed toxicity. Our work also compares the impact of competitive and noncompetitive BSEP inhibition for CP-724,714 and demonstrates that noncompetitive inhibition leads to much greater bile acid accumulation and potential toxicity. Our research demonstrates the potential for mechanistic modeling to contribute to the understanding of how bile acid transport inhibitors cause DILI.

**Keywords: BSEP inhibition, bile acids and salts, drug-induced liver injury, mechanistic modeling, bosentan, CP-724,714, population variability**

## **INTRODUCTION**

Inhibition of the bile salt export pump (BSEP) by a drug has been implicated as a risk factor for the drug's potential to cause drug-induced liver injury (DILI) (Dawson et al., 2011; Morgan et al., 2013). Several high-profile DILI-causing drugs, for example troglitazone (Funk et al., 2001; Smith, 2003), bosentan (Fattinger et al., 2001; Eriksson et al., 2011), and nefazodone (Kostrubsky et al., 2006), have been shown to be BSEP inhibitors. Furthermore, prediction of toxicity by these molecules has been uneven; for example, neither bosentan nor troglitazone displayed toxicity in animal models (Leslie et al., 2007; Lauer et al., 2009). Predictions involving hepatobiliary transporter IC50 values in *in vitro* assays have shown better predictive ability (Dawson et al., 2011; Morgan et al., 2013). Improving the ability to predict the frequency and severity of DILI with BSEP inhibitors will allow those involved in drug development to better judge the risk involved with moving a drug in development into the clinic or beyond early-stage clinical trials.

DILIsym® is a multi-scale mechanistic model incorporating numerous functions of the liver and disruptions of the function with the goal of predicting the DILI potential of drugs at various stages in the development process (Howell et al., 2012, 2014; Woodhead et al., 2012; Shoda et al., 2014). Previously, we have constructed and validated a model of bile acid homeostasis and transporter inhibition within DILIsym® (Woodhead et al., 2014). We found that inhibiting BSEP could lead to significant increases in bile acid concentrations in the liver, and that the effects of bile acid transporter inhibitors should be considered on a simulated population as well as on a single baseline individual. We have expanded that model to include a representation of bile acidmediated toxicity based on experiments performed by Yang et al. (2013). In that work, we constructed a relationship between intracellular bile acid concentration and cellular ATP. We used this relationship to predict cellular necrosis using the existing relationship between ATP and cell death in DILIsym®. This model has been used previously to effectively predict the frequency and timing of ALT elevations observed with troglitazone in clinical trials, and has predicted the difference between troglitazone and the similar non-toxic drug pioglitazone (Yang et al., 2014).

In the present work, we use this model of bile acid-mediated cytotoxicity to model the observed hepatotoxicity of bosentan, telmisartan, and CP-724,714. Bosentan is a currently marketed medication for pulmonary arterial hypertension that carries a black-box warning for hepatotoxicity. In clinical trials a dose of 1000 mg/day of bosentan caused between 8 and 18% of individuals to experience increases in serum ALT greater than 3-fold (Fattinger et al., 2001). Bosentan has also been shown to be a relatively potent BSEP inhibitor. Telmisartan, an angiotensin II receptor agonist used in hypertension treatment, is also a relatively potent BSEP inhibitor, but has not been reported to cause hepatotoxicity in humans (Morgan et al., 2013). CP-724,714 is an anti-cancer drug that was terminated from development after Phase II clinical trials revealed liver signals (Guo et al., 2008). CP-724,714 has been shown to inhibit multiple transporters, including BSEP, in addition to being a mitochondrial toxin (Feng et al., 2009). In this report we explore the DILIsym® software's (version 2C) predictions for the toxicity of bosentan in humans, and the lack of toxicity of bosentan in rats and telmisartan in humans. We will also test several theories about the toxicity of CP-724,714 in order to suggest potential viable explanations for the observed toxicity in early clinical trials.

## **METHODS**

The construction and validation of the bile acid homeostasis model in DILIsym® is described in a previous paper (Woodhead et al., 2014). DILIsym® models the synthesis and enterohepatic recirculation of two main potentially toxic bile acids, chenodeoxycholic acid (CDCA) and lithocholic acid (LCA), and their amide (and sulfate in the case of LCA) conjugates. DILIsym® describes the intrahepatic accumulation of these toxic bile acids as well as the concentrations in the gallbladder, portal blood, gut lumen, and systemic blood. Concentrations of bile acid transporter inhibitors are modeled using a physiologically-based pharmacokinetic model (PBPK) described in depth in previous papers (Howell et al., 2012; Woodhead et al., 2012). The bile acid concentrations are linked to ATP decline as described below; the effect of ATP decline on eventual hepatocyte necrosis is described by a model of the hepatocyte life cycle also described in previous papers (Howell et al., 2012; Woodhead et al., 2012; Shoda et al., 2014). A diagram of the interaction between the various submodels of DILIsym® employed in this paper can be seen in **Figure 1**.

The DILIsym® model of bile acid toxicity is based on *in vitro* experiments comparing the intracellular level of LCA and CDCA to cellular ATP levels (Yang et al., 2013). In order to construct the connection between bile acid accumulation and toxicity, a small model of the DILIsym® ATP turnover model was constructed (Yang et al., 2013); a diagram of this model is shown in **Figure 2**. The relationship between intracellular ATP and intracellular bile acids was modeled by the following equation:

$$\frac{d[ATP]}{dt} = k\_{usage} - k\_{synth}S$$

where *kusage* is the rate of *ATP* usage in the cell, k*synth* is the rate of ATP synthesis in the cell, and *S* is the signal for ATP synthesis inhibition by bile acids. *S* has Hill type behavior and is given by the following equation:

$$S = \frac{1}{1 + \frac{V\_{\text{max},S}[BA]\_{delay}^H}{K\_{m,S}^H + [BA]\_{delay}^H}}$$

where *V*max*,<sup>S</sup>* is the maximum denominator for the signal, *Km,<sup>S</sup>* is the Michaelis constant for the signal, and *H* is the Hill coefficient for the signal. [*BA*]*delay* is the delayed intracellular bile acid concentration, given by the standard delay equation:

$$\frac{d[BA]\_{delay}}{dt} = \text{tr}\left( [BA] - [BA]\_{delay} \right).$$

The delay constant *τ* , *V*max*,S*, *Km,S*, and *H* were the parameters that were fitted to the ATP time course from the experiment. These parameters were then applied to the intracellular bile acid concentrations modeled in DILIsym®. We applied the parameters for the unconjugated bile acids measured in the experiments to the conjugated bile acids in DILIsym®; for example, we used the unconjugated LCA parameters to describe the toxicity of LCA amide and sulfate conjugates. While some data exist describing the relationship between intracellular amide- and sulfateconjugated bile acids and toxicity (Chatterjee et al., 2013), these are not enough to justify using different toxicity parameters for conjugated and unconjugated bile acids; this is an area of potential refinement for DILIsym®.

Physiologically-based pharmacokinetic (PBPK) models were constructed for bosentan, telmisartan, and CP-724,714 in humans using available *in vivo* time course data (Weber et al., 1996; Stangier et al., 2000; Munster et al., 2007). For simulations where bosentan metabolism was impaired, the *V*max for the metabolism of both metabolites was represented as 10% of their baseline value. Both bosentan and telmisartan are active transport substrates; the active transport parameters for both were adapted for use in DILIsym® from Ménochet et al. (2012). The bosentan PBPK model includes the major and minor metabolite, with metabolism parameters based on published *in vitro* metabolism data (Ubeaud et al., 1995); these were included because the minor metabolite of bosentan also inhibits BSEP (Fattinger et al., 2001). The major metabolite of telmisartan, a glucuronide, has not been shown to inhibit BSEP and so was not included in the PBPK model. CP-724,714's complex metabolism was not modeled explicitly; however, a single metabolite compartment in the liver was included in order to test the hypothesis of whether metabolite accumulation could explain CP-724,714 toxicity. PBPK model results for human bosentan are shown in **Figure 3**; those for human telmisartan are shown in **Figure 4**; and those for human CP-724,714 in **Figure 5**. For bosentan in rats, data did not exist in the literature. Data provided to us from Amgen, Inc. (Thousand Oaks, CA) were used to parameterize the bosentan rat PBPK model shown in **Figure 6**.

Inhibition constants for bosentan were taken from published research that found *Ki* values for both bosentan (12µM) and the minor metabolite of bosentan (8.5 µM) (Fattinger et al., 2001).

The inhibition of BSEP by bosentan and its metabolite were modeled as noncompetitive in agreement with the data. These data are from rat Bsep vesicles; we assigned the same *Ki* value to humans, as there is evidence suggesting that rat Bsep and human BSEP are inhibited with similar potency by bosentan (Mano et al., 2007). Bile acid uptake inhibition constants for bosentan were taken from work done by Leslie et al. (2007). The telmisartan *Ki* was approximated using the IC50 value reported in the literature (Morgan et al., 2013). A list of *Ki* values used in our simulations is included in **Table 1**. In addition, bosentan is a potent inducer of its own uptake into the liver and metabolism (Dingemanse and van Giersbergen, 2004); these effects are included in our PBPK model. It is important to note that for each drug, only the parameters related to drug pharmacokinetics and transporter inhibition

were changed; all other remaining parameters were not changed from drug to drug.

For CP-724,714 IC50 values, experiments were run in rat Bsep, dog Bsep, and human BSEP-expressing vesicles. CP-724,714 was synthesized by Pfizer, Inc. (Groton, CT), while radiochemicals were purchased from Perkin Elmer (Perkin Elmer, Waltham MA). Other chemicals were purchased from Sigma-Aldrich (Sigma-Aldrich, St Louis, MO) or were of analytical grade. Membrane vesicles were obtained from Sf9 cells not expressing and expressing human BSEP (Solvo Biotechnology, Szeged, Hungary), rat Bsep (AB Life Technologies, Waltham, MA), or dog Bsep (AB Life Technologies).

CP-724,714 was incubated (2 or 5 min) with membrane vesicle preparations (total protein: 50µg/well) and the probe substrate, taurocholate (2µM). Serial dilutions of CP-724,714 (30 mM stock, 10 mM, 3.16 mM, 1 mM, 316µM, 100µM, 31.6µM, 10µM) were prepared in DMSO. Incubations in duplicate were carried out in the absence or presence of 4 mM ATP to distinguish between transporter mediated uptake and passive diffusion into the vesicles. CP-724,714 was added to the reaction mixture in 0.75µl of solvent (1% of the final incubation volume). Glyburide (100 – 0.1µM) was used as the positive control inhibitor. Reaction mixtures were pre-incubated for 10 min at 37◦C. Reactions were started by the addition of 25µl of 12 mM MgATP (or AMP, as disodium salt, for background controls), preincubated separately. Reactions were stopped by the addition of 200µl of ice-cold washing buffer and immediate filtration via glass fiber filters mounted to a 96-well plate (filter plate). The filters were washed, dried and the amount of substrate inside the filtered vesicles determined by liquid scintillation. Maximal observed inhibition ranged between 94.8 and 99.8% for CP-724,714. Further information on this method can be found in the literature (Kis et al., 2009).

Bosentan dosing in humans was simulated as a 500-mg twicedaily dose for 30 days. The doses were given once every 12 h; this is one of the dosing regimens in the clinical trials that demonstrated toxicity as reported by Fattinger et al. (2001). Telmisartan dosing was simulated as a 50, 3000, or 12,000-mg dose once daily for 30

days; common clinical dosing regimens range from 40 to 80 mg per day (Meredith, 1999). CP-724,714 was dosed three times daily, or once every 8 h. A range of CP-724,714 doses were simulated in order to cover the range of exposure levels reported by Guo et al. (2008). Three simulated meals per day were included in all human simulations, and simulated to occur concurrent with drug dosing. Bosentan dosing in rats was simulated as a 50 mg/kg dose once daily.

DILIsym® contains simulated populations, or SimPops™, that represent a plausible range of variability in key model parameters. The human SimPops™ and rat SimPops™ used in the course of these simulations contain variability in several bile acid transport parameters and are described in our previous work (Woodhead et al., 2014; Yang et al., 2014). There are 331 individuals in the human SimPops™ and 191 individuals in the rat SimPops™; although the numbers of individuals are different, both the human and the rat SimPops™ are designed to account for the entire plausible range in variability in bile acid transport parameters and bile acid profiles observed in a sample population (Woodhead et al., 2014). Pharmacokinetic variability in bosentan disposition was also included in one SimPops™ simulation; for this population, a normal distribution around four variables (oral absorption coefficient, liver uptake *V*max, and major and minor metabolite formation *V*max) was superimposed upon the existing small human SimPops™. The range of these parameter values was based on a known range of variability for hepatobiliary transporters (Meier et al., 2006) and upon previous simulations where pharmacokinetic variability was considered (Woodhead et al., 2012).

## **RESULTS**

Results from the SimPops™ analysis with bosentan are presented in **Table 2**. When we simulated bosentan in the human SimPops™ using the baseline assumptions (noncompetitive inhibition, no basolateral inhibition, normal metabolism), 1 individual out of 331 developed ALT elevations greater than 3-fold above the baseline value. When potential variability in pharmacokinetics was incorporated into the SimPops™, the predicted incidence rate rose to 3 out of 331. While this successfully predicts the potential toxicity of bosentan, the incidence rate is well below the incidence rate of 8–18% observed during the clinical trials (Fattinger et al., 2001).

A genetic polymorphism that limits CYP metabolism has been shown to correlate with an increased rate of toxicity from bosentan exposure (Markova et al., 2013). We have simulated this case using the human SimPops™ by decreasing the bosentan metabolism *V*max values for both minor and major metabolites by 10-fold. This change led to eight simulated individuals out of 331 (2.42%) developing ALT elevations. DILIsym® correctly predicts the increased risk of toxicity due to this potential genetic polymorphism, though the predicted incidence rate is still less than the overall rate in the general population from the clinical study.

When bosentan was also simulated in the rat SimPops™, zero individual rats out of 191 developed ALT elevations. This is consistent with preclinical observations which have reported no toxicity in the rat (Leslie et al., 2007). It has been suggested that the difference in bile acid uptake inhibition by bosentan between rats and humans contributes to the species difference in hepatotoxicity (Ansede et al., 2010). In order to test this theory, we ran a simulation in the rat SimPops™ with uptake inhibition eliminated (*Ki* set to 1 <sup>×</sup> <sup>10</sup>10); though bile acid accumulation in the liver was predicted (results not shown), no toxicity was observed. This suggests that the difference in the bile acid pool between humans and rats is likely a more significant contributor to the species difference in toxicity.

In order to ensure that the model can differentiate non-toxic BSEP inhibitors from toxic BSEP inhibitors, we modeled telmisartan in the human SimPops™. No toxicity was predicted in the human SimPops™, which is consistent with clinical observations. This was true whether telmisartan was modeled as a competitive or a noncompetitive inhibitor.

While there is scant difference between predicting one individual developing toxicity and predicting zero individuals developing toxicity (out of 331), a fuller understanding of the risk of a given compound can be gained by investigating more mechanistic data within DILIsym®. **Figure 7** displays the minimum hepatic ATP for each individual in the human and rat SimPops™ for bosentan and the human SimPops™ for telmisartan (modeled as both a competitive and non-competitive inhibitor of BSEP). We can see that while only one simulated individual in the human SimPops™ had an ATP decline after bosentan dosing that ultimately led to toxicity, many more had ATP reductions visibly below the baseline value. This was not true of the rat, or of telmisartan in either case; the simulated individuals all have hepatic ATP values very close to the baseline value. This places the difference between bosentan and telmisartan, and the difference between rat and human, in sharper relief and demonstrates the potential of the model for predicting species differences in toxicity and differentiating between toxic and non-toxic BSEP inhibitors.

The simulation results for CP-724,714 in the baseline human individual, and a comparison of these results to the clinical observations in Guo et al. (2008), are shown in **Figure 8**. The graph compares the normalized liver function test (LFT) elevation

## **Table 2 | SimPops™ simulation results for simulated ALT elevations caused by bosentan and telmisartan.**


**Table 1 | Inhibition constants used for the compounds in the simulations conducted for this paper.**


*\*Denotes that the Ki was approximated using the listed IC50 value for which mode of inhibition was not determined.*

**FIGURE 7 | ATP levels in individuals after simulated dosing of bosentan to the human SimPops™ (A) and the rat SimPops™ (B); and with simulated telmisartan dosing to the human SimPops™ when represented as a**

**competitive (C) and noncompetitive (D) inhibitor.** Each point represents an individual within the human SimPops™; the baseline ATP concentration in humans is 4.2 mM, while in rats the baseline ATP concentration is 2.0 mM.

**measured by AUC0-24 on Day 1 of Cycle 2 (Day 22 overall) in the baseline human in DILIsym® and the clinical trial patients from Guo et al. (2008).** The lines refer to the dose response simulated in DILIsym® under various assumed conditions listed in the table legend. "Baseline assumptions" refers to the case where only the parent CP-724,714 inhibits BSEP competitively and does not inhibit the electron transport chain (ETC). "Metabolite inhibition" refers to the case where the metabolite of

"Poor clearance" refers to the case where the metabolite biliary clearance value is set to a value 10-fold lower than the parent biliary clearance value. "Noncompetitive" refers to the case where all BSEP inhibition is modeled as noncompetitive. "Added ETC inhibition" refers to the case where inhibition of the electron transport chain was simulated for both parent and metabolite. The points are individual patients from Guo et al. (2008). The purple dashes are the individuals in a human SimPops™.

against the AUC of the drug in the patient's bloodstream on Day 22 (Day 1 of Cycle 2) of drug dosing, a measure of steady-state drug exposure. The normalized LFT elevation used by the Guo et al. (2008) paper to which our simulations are compared is given by the following expression (Guo et al., 2008):

$$\max\left(\frac{\text{ALT}\,\text{fold}\,\text{increase}}{5}, \frac{\text{blind}\,\text{in}\,\text{fold}\,\text{in}\,\text{increase}}{3}\right)$$

The simulation results are the maximum normalized LFT elevation at each dose simulated for each individual case. The graph shows that bile acid inhibition alone cannot explain the clinical toxicity, as liver function tests did not elevate in the simulations when bile acid inhibition alone was simulated. This was true for both competitive and noncompetitive inhibition. However, the inclusion of mitochondrial toxicity did, in fact, lead to a toxic response that was augmented when bile acid accumulation also occurred. The data are best described here by this combination of competitive inhibition and electron transport chain (ETC) inhibition.

The simulations also predict that both BSEP and mitochondrial ETC inhibition by CP-724,714's major metabolite are necessary to explain the toxicity, since CP-724,714 is rather rapidly metabolized by the liver and so parent residence time within the liver is somewhat limited. When the potential activity of the metabolite is not included, no toxicity is shown no matter the mechanism or combination of mechanisms selected. The accumulation of this metabolite is also necessary to explain the toxicity; a biliary clearance of 10% of the value for the parent compound was used in order to reproduce the clinical toxicity.

SimPops™ results for CP-724,714 in humans are also shown in **Figure 8**, with each individual in the SimPops™ represented by a purple dash. We found that while the simulated baseline individual did not display the toxicity that would have been expected from the clinical data, the population sample did contain several individuals who developed clinically-relevant ALT elevations, but only if the BSEP inhibition was noncompetitive. Furthermore, the range of injury in the SimPops™ is far wider than the range reported in the clinical study; the most severe normalized LFT elevation in our simulated population was 6.5, while the largest clinical normalized LFT in that exposure range was about 2. Twelve individuals in our 331-individual SimPops™, or 3.62%, developed toxicity; this is lower than the 36% of individuals in the exposure range near the simulated dose who developed LFT elevations. No individuals in the SimPops™ developed toxicity if the inhibition was competitive, demonstrating the importance of differentiating between modes of inhibition when determining BSEP inhibition constants.

Uncertainty about the biliary clearance of the CP-724,714 metabolite has an impact on simulated ALT elevations. **Figure 9** shows simulated ALT elevations at a constant dose of CP-724,714 when the biliary clearance of the metabolite is modulated. This variable does not affect the plasma pharmacokinetics of parent CP-724,714; it is thus a degree of freedom in DILIsym®. However, as **Figure 9** shows, it has a significant effect on the predicted toxicity of CP-724,714; DILIsym® suggests, therefore, that experiments clarifying the amount of hepatic accumulation and clearance of CP-724,714's major metabolite should be conducted in order to fully elucidate the molecule's toxicity.

### **DISCUSSION**

We have used DILIsym® to model bosentan and found that, using the mechanistic simulation results calculating ATP decline, DILIsym® suggests that the potential for toxicity in a human population is greater than that in a rat population. Furthermore, the same ATP decline simulation results show a clear difference between the human population response to bosentan and

to telmisartan. While the prediction of the toxic potential of bosentan was well below the actual clinical incidence rate, these results nevertheless show promise for the ability to use a mechanistic mathematical model of DILI, rather than small-animal models, to predict the human toxicity of a BSEP inhibitor. More revealingly, the modeling exercise suggested a potential reason for this discrepancy in toxicity between rats and humans. Previous work in this area demonstrated the difference in uptake transporter inhibition between rats and humans (Leslie et al., 2007; Ansede et al., 2010) and proposed that this was a contributing factor to the species difference. However, our simulations suggest that removing the uptake inhibition from the rat model altogether does not lead to liver toxicity in the rat. This suggests that the difference in bile acid pools and metabolic pathways are most likely the strongest contributor to the species difference in toxicity. The rat bile acid pool has more non-toxic bile acids, such as cholic acid and muricholic acid, and less of the toxic bile acids CDCA and LCA than the human (Hofmann, 2009; García-Cañaveras et al., 2012). Furthermore, the rat can hydroxylate LCA into the less-toxic hyodeoxycholic acid, which is a pathway that does not occur in the human (Hofmann, 2004). Because of these species differences in bile acid metabolism, our modeling questions the utility of current small-animal models in the prediction of BSEP inhibitor-mediated toxicity in humans.

Our results comparing the ATP decrement in the liver caused by bosentan and telmisartan demonstrate that DILIsym® can differentiate between a toxic and a non-toxic compound with a similar IC50. In the case of bosentan and telmisartan, the difference is mostly a pharmacokinetic one; telmisartan is dosed at a far lower dose than is bosentan. However, our results with increasing telmisartan dose suggest that a much higher dose of telmisartan, compared to bosentan, is necessary to cause toxicity in a human SimPops™. The suggestion is that dose is not necessarily predictive on its own; differences in metabolism and in liver uptake transporter affinity and capacity both are likely contributors to the differences in toxicity between bosentan and telmisartan.

While we qualitatively predict the species difference in bosentan toxicity and the toxicity of bosentan compared with the safety of telmisartan, we substantially underpredict the incidence rate of the toxicity of bosentan in the general population. Clinical trials predict 8–18% toxicity at the 1000 mg/day dosing level; DILIsym® predicted a much lower incidence rate even when individuals with metabolic polymorphisms were considered.

There are several possible reasons for this discrepancy, mostly related to the possibility of other mechanisms of toxicity not currently included in the DILIsym® model of bosentan. First, we do not represent the inhibition of basolateral bile acid transport by bosentan; this could prevent a potential clearance pathway for toxic bile acid species and lead to more toxicity. Recent research has shown that bosentan has the potential to inhibit MRP4, a basolateral bile acid transporter (Morgan et al., 2013). Exploratory simulations representing bosentan as a basolateral transporter inhibitor (reported in **Table 2**) suggest that this explanation is likely only part of the problem; if the bosentan basolateral transporter *Ki* is particularly low, however, this could account for some of the discrepancy in predicted incidence rates.

Second, DILIsym® does not represent bile acid toxicity in a mechanistic manner. Recent research suggests that the bile acids affect the mitochondria and potentially lead to the mitochondrial membrane permeability transition (Schulz et al., 2013); representing this effect in DILIsym® could introduce extra sources of variability and could lead to a higher predicted toxicity incidence rate. Indeed, work is underway in this area for DILIsym® and preliminary results suggest that the mechanistic model increases the predicted rate of bosentan toxicity without predicting toxicity in the rat; this work will be the focus of a subsequent paper.

Third, bosentan could cause toxicity through mechanisms not currently included in the model for bosentan or not currently represented in DILIsym®. Bosentan has been shown to inhibit phospholipid transport in rodents (Fouassier et al., 2002). Phospholipids protect the bile duct from the cytotoxic effect of bile acids; inhibition of phospholipid transport has been shown to lead to liver toxicity in the case of itraconazole (Yoshikado et al., 2011). While bosentan's mechanism of action is dissimilar to that of itraconazole, this shows that the ability of bosentan to interfere with phospholipid transport could be at least partially responsible for the observed toxicity of bosentan. Phospholipid transport and the toxic effect of its inhibition are not currently represented in DILIsym®; this is an area for potential future research and refinement.

Also, bosentan's effects on the mitochondria have not been elucidated. Previous research has shown that there is a correlation between the ability to inhibit BSEP and toxic effects on the mitochondria, and that this combined effect is itself correlated with DILI risk (Aleo et al., 2014). Furthermore, our simulations with CP-724,714 demonstrated that the toxicity seen in the clinic could not be explained without representing the drug's mitochondrial effects; while our research suggests that it is plausible that a compound could cause toxicity through bile acid transporter inhibition alone, it is also plausible that the toxicity of bosentan is due to a combined mitochondria/bile acid toxicity mechanism. Indeed, our modeling showed that this was most likely the case with CP-724,714. Research in this area is underway in our group and will focus on potential synergy between mitochondrial toxicity and bile acid buildup in the liver.

Fourth, DILIsym® does not represent the intracellular trafficking and potential localization of the drug within the hepatocyte. Further research in this area could help us improve DILIsym® and thus our prediction of toxicity incidence rates.

Our simulations with CP-724,714 demonstrate the ability of mechanistic modeling to consider and prioritize multiple hypotheses when modeling compounds where the cause of toxicity is not fully established. Potential variability in the degree of hepatic accumulation of a compound that inhibits BSEP is particularly important, as shown by **Figure 9**. This is an especially salient point to consider given that the toxicity of compounds such as troglitazone are suspected to be due to the hepatic accumulation of a BSEP-inhibiting metabolite (Masubuchi, 2006). Furthermore, the difference between competitive and non-competitive inhibition, the potential contribution of ETC inhibition to the observed toxicity, and the toxic potential of CP-724,714's metabolites are all sources of uncertainty that required a hypothesis-based approach to modeling. While this is of limited value for predicting toxicity on its own, the hypothesisbased modeling approach is potentially valuable in determining which experiments would be most impactful in using the model to properly predict toxicity. In the case of CP-724,714, these experiments include hepatic accumulation studies, BSEP inhibition studies, and ETC inhibition studies on CP-724,714's metabolites.

## **ACKNOWLEDGMENTS**

We would like to acknowledge Yurong Lai (Pfizer) for assistance with the BSEP inhibition studies, Michael D. Aleo and Denise Robinson-Gravatt for discussions around DILIsym® implementation for CP-724,714, and the product data management and project team scientists whose CP-724,714 data facilitated this modeling effort. We would also like to acknowledge Ryan Morgan (Amgen, Inc.) for many helpful discussions and assistance with the direction of this research. We also acknowledge Ryan Morgan, Chuck Qualls, Hisham Hamadeh, Sabina Buntich, Nataraj Kalyanaraman and Dean Hickman (Amgen) and Andy Mould and Jim Howard (Covance Laboratories) for their work on the bosentan pharmacokinetic rat studies that have been used to optimize the rat bosentan model.

## **REFERENCES**


toxicity in rat sandwich-cultured hepatocytes: incorporation into a mechanistic model of drug-induced liver injury. *Toxicologist* 132:226.

Yoshikado, T., Takada, T., Yamamoto, T., Yamaji, H., Ito, K., Santa, T., et al. (2011). Itraconazole-induced cholestasis: involvement of the inhibition of bile canalicular phospholipid translocator MDR3/ABCB4. *Mol. Pharmacol*. 79, 241–250. doi: 10.1124/mol.110.067256

**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.

*Received: 13 June 2014; accepted: 22 October 2014; published online: 07 November 2014.*

*Citation: Woodhead JL, Yang K, Siler SQ, Watkins PB, Brouwer KLR, Barton HA and Howell BA (2014) Exploring BSEP inhibition-mediated toxicity with a mechanistic model of drug-induced liver injury. Front. Pharmacol. 5:240. doi: 10.3389/fphar. 2014.00240*

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

*Copyright © 2014 Woodhead, Yang, Siler, Watkins, Brouwer, Barton and Howell. 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.*

## Applications of linking PBPK and PD models to predict the impact of genotypic variability, formulation differences, differences in target binding capacity and target site drug concentrations on drug responses and variability

## *Manoranjenni Chetty1\*, Rachel H. Rose1, Khaled Abduljalil 1, Nikunjkumar Patel 1, Gaohua Lu1, Theresa Cain1, Masoud Jamei <sup>1</sup> and Amin Rostami-Hodjegan1,2*

*<sup>1</sup> Simcyp Limited (a Certara Company), Blades Enterprise Centre, Sheffield, UK*

*<sup>2</sup> Manchester Pharmacy School, University of Manchester, Manchester, UK*

#### *Edited by:*

*Sergey Ermakov, Bristol-Myers Squibb, USA*

#### *Reviewed by:*

*Valeria Bruno, University of Rome "Sapienza," Italy Georg Hempel, University of Münster, Germany*

#### *\*Correspondence:*

*Manoranjenni Chetty, Simcyp Limited (a Certara Company), Blades Enterprise Centre, John Street, South Yorkshire, Sheffield S2 4SU, UK e-mail: m.chetty@simcyp.com*

knowledge-rich physiologically based pharmacokinetic (PBPK) models with pharmacodynamics (PDs) models. Four distinct applications that were developed and tested are presented here. PBPK models were developed for metoprolol using different CYP2D6 genotypes based on *in vitro* data. Application of the models for prediction of phenotypic differences in the pharmacokinetics (PKs) and PD compared favorably with clinical data, demonstrating that these differences can be predicted prior to the availability of such data from clinical trials. In the second case, PK and PD data for an immediate release formulation of nifedipine together with *in vitro* dissolution data for a controlled release (CR) formulation were used to predict the PK and PD of the CR. This approach can be useful to pharmaceutical scientists during formulation development. The operational model of agonism was used in the third application to describe the hypnotic effects of triazolam, and this was successfully extrapolated to zolpidem by changing only the drug related parameters from *in vitro* experiments. This PBPK modeling approach can be useful to developmental scientists who which to compare several drug candidates in the same therapeutic class. Finally, differences in QTc prolongation due to quinidine in Caucasian and Korean females were successfully predicted by the model using free heart concentrations as an input to the PD models. This PBPK linked PD model was used to demonstrate a higher sensitivity to free heart concentrations of quinidine in Caucasian females, thereby providing a mechanistic understanding of a clinical observation. In general, permutations of certain conditions which potentially change PK and hence PD may not be amenable to the conduct of clinical studies but linking PBPK with PD provides an alternative method of investigating the potential impact of PK changes on PD.

This study aimed to demonstrate the added value of integrating prior *in vitro* data and

**Keywords: PBPK linked PD models, CYP P450 genotypes and response, heart drug concentration and QTc, formulation effects on drug response, target site concentrations and response**

## **INTRODUCTION**

Physiologically based pharmacokinetic (PBPK) modeling provides a mechanistic platform for the integration of the concentrationtime profile of the drug with realistic physiological and biological processes in the body. This modeling approach offers an advantage over traditional compartmental modeling approaches since it potentially allows for extrapolation and further investigation into conditions for which pharmacokinetic (PK) studies have not been conducted, thereby informing and accelerating the drug development process. Predictions on drug–drug interactions, first in man dosing, optimal clinical study designs, dosage requirements for drugs that are metabolized by polymorphic enzymes and dosage adjustments in disease states are some of the PBPK applications that could potentially be used during the drug development and regulatory submission processes (Chen et al., 2012; Huang and Rowland, 2012; Rostami-Hodjegan, 2012; Sinha et al., 2012; Zhao et al., 2012; Rowland, 2013; Vieira et al., 2014).

Since the primary concern in drug development is the efficacy and safety of the drug, PBPK linked pharmacodynamic (PD) models can be very valuable, offering a platform for exploring the effect of variability in various physiological, biochemical, and formulation factors on the response to the drug, especially where clinical studies have not or cannot be conducted. An important advantage of PBPK linked PD models is the ability to link the drug concentration at the probable site of action with toxicological and/or therapeutic effect. This is especially important when the plasma concentration is not a good surrogate for the concentration at the site of drug action (Rostami-Hodjegan, 2013), as demonstrated in a recent study with rosuvastatin, a cholesterol lowering drug which is a substrate of the OATP1B1 influx transporter (Rose

et al., 2014). A PBPK/PD model that used liver concentrations of rosuvastatin demonstrated improved ability to capture the effect of the OATP1B1 c.521T > C single nucleotide polymorphism on the change in cholesterol synthesis rate in response to rosuvastatin, compared to a model using plasma concentrations of the drug. Based on significant plasma concentration differences, a dosage adjustment in rosuvastatin may have been considered in patients with the OATP1B1 c.521T > C polymorphism, but a clinical study had shown no significant differences in response in the two OATP1B1 phenotypes. However, using PBPK modeling with a PD model driven by free liver concentrations of rosuvastatin, no net difference in liver concentrations of the drug in patients with the polymorphism was observed and their clinical response was similar, suggesting that dosage adjustments were unnecessary.

Many other scenarios in drug development and clinical practice can benefit from incorporating the prior systems knowledge into PBPK models when these are linked to PD modeling. Four such examples from distinct areas are presented in this paper to demonstrate the wide range of applications of this approach. The first case study considers CYP genotypes which can have a significant effect on drug PK. This case study is a quantitative prediction of the variation in the clinical response (measured as heart rate) to standard doses of metoprolol in ultrarapid metabolizers (UMs), extensive metabolizers (EMs), and poor metabolizers (PMs) of CYP2D6. The PBPK model used differences in CYP2D6 abundance obtained from *in vitro* studies to simulate phenotypic differences in metoprolol PK. The second case study explores the potential for the application of PBPK/PD modeling to predict the response to a controlled release (CR) formulation using *in vitro* data for the CR formulation. A PBPK/PD model was developed and verified for a nifedipine immediate release (IR) formulation. It was then used for prediction of the PK and PD profiles of the CR formulation using only the dissolution profile of the CR formulation. The third case study investigates the application of a PBPK model linked to a semi-mechanistic PD model developed for one drug to predict the response to a second drug based on clinical data for the first drug that acts on the same target. A PBPK/PD model was developed and verified for triazolam using an operational agonism PD model. A zolpidem specific target binding parameter obtained by *in vitro* studies was then used with the triazolam model to predict response to zolpidem. Such models are useful for comparing several potential drug compounds that belong to the same therapeutic class, when clinical data is available for just one of the compounds. In the fourth case study, the PD model is driven by target site drug concentrations and used to gain a mechanistic understanding of a clinical observation. The differences in the potential for QT prolongation by drugs such as quinidine in Caucasian and Asian females are well known, although the reason for this difference has not been established. Using PBPK models with heart concentrations of quinidine as the input to PD models, a higher sensitivity to heart concentrations of quinidine was demonstrated in Caucasian females.

## **MATERIALS AND METHODS**

The Simcyp population-based simulator (V12 www.simcyp.com; Jamei et al., 2009) was used in the development of the PBPK/PD models and the simulations. Simcyp compound files (parameters

listed in **Tables 1–4**) that were validated previously (Howgate et al., 2004; McGinnity et al., 2008; Polasek et al., 2010; Rowland Yeo et al., 2004; Patel et al., 2014) and population data (Jamei et al., 2014) available in Simcyp (V12) were used. Clinical data used in the case studies were digitized from published clinical studies using the Getdata software. Simulations using the developed models were verified by comparison with clinical data prior to further predictive applications. Simulations were found to be acceptable if the predicted parameters were within two fold of the observed data (Guest et al., 2011) and visual predictive checks showed observed data within the 5 and 95% percentiles of the predicted data. Methodological details relevant to the individual case studies are described below.

## **CASE STUDY 1: ASSESSING THE IMPACT OF GENOTYPICALLY CONTROLLED ELIMINATION**

Although plasma concentrations of metoprolol and effects on heart rate have been shown to correlate significantly with CYP2D6 metabolic phenotype in clinical studies (Kirchheiner et al., 2004;



#### **Table 2 | Parameters used for the Nifedipine PBPK model (Simcyp V12).**

#### **Table 3 | Parameters used for the triazolam PBPK model (Simcyp V12).**



Sharma et al., 2005), the prevalence of some phenotypes may not be adequately high in a study population to discern the differences in PK and PD. Therefore, it would be of value to use the prior *in vitro* information on metabolism together with PK and PD information in prevalent phenotypes of CYP2D6 to conduct virtual clinical studies with a view to assess the potential pharmacological differences in various less frequent phenotypes, prior to the conduct of clinical studies or in lieu of such studies when the studies are not feasible and yet providing

a recommendation is more prudent than leaving a void in prescribing information.

The reduction in heart rate due to a standard 100 mg dose of metoprolol in virtual healthy Caucasian populations was simulated and stratified for their CYP2D6 phenotypes. The

#### **Table 4 | Parameters used for the zolpidem PBPK model (Simcyp V12).**

#### **Table 5 | Parameters used for the Quinidine PBPK model (Simcyp V12).**


Simcyp metoprolol compound file (**Table 1**) was used with a minimal PBPK model, first order absorption and elimination by enzyme kinetics. The study design was matched to that of


Kirchheiner et al. (2004). Simulated contribution of the CYP2D6 phenotypes (EM, PM, and UM) to metoprolol PK within Simcyp is based on the propagation of the differences in CYP2D6 abundance, obtained from *in vitro* data. Concentration-time profiles published by Kirchheiner et al. (2004) and Sharma et al. (2005) were compared with the predicted profiles and PK parameters. PD differences in the different phenotypes were assessed by an Emax model published by Kirchheiner et al. (2004) and assumed to be the same regardless of CYP2D6 genotype. Based on this direct effects PD model, response was propagated via changes in the plasma concentration profile. PD simulations were compared with clinical observations from the above two studies, to verify that the model predicted the PD corresponding to the different phenotypes adequately.

## **CASE STUDY 2: ASSESSING THE CONSEQUENCES OF MODIFYING THE DRUG FORMULATION**

Nifedipine is a dihydropyridine calcium channel blocker commonly used in the treatment of hypertension and exerts its hypotensive effect primarily through arterial dilation. CR formulations are now recommended in the treatment of hypertension as they have been shown to offer a number of clinical benefits over IR nifedipine (Reitberg et al., 1987; Schug et al., 2002; Meredith and Elliott, 2004; Wonnemann et al., 2006). This study aimed to integrate the PBPK models describing the plasma profile of IR and that for the CR nifedipine (Shimada et al., 1996; Brown and Toal, 2008) with the PD model available for IR nifedipine (Shimada et al., 1996), to identify whether it could be extrapolated to predict the response to nifedipine GITS, a CR formulation that is reported to achieve a zero order release rate sustained over 24 h, through an osmotic release mechanism (Brown and Toal, 2008).

The PK and PD profiles for nifedipine in the treatment of hypertension were simulated using the Simcyp nifedipine compound file (**Table 2**), a minimal PBPK distribution model and elimination by enzyme kinetics. To simulate the PK profile of IR nifedipine the first order absorption model was used, while for the nifedipine GITS [a CR formulation reported to achieve a zero order release rate sustained over 24 h through an osmotic release mechanism (Brown and Toal, 2008)] formulation effects were described by a mechanistic absorption model within Simcyp (ADAM) using *in vitro* data (dissolution data) for the CR profile and intrinsic solubility of nifedipine (Janssen Therapeutics, 2013). The PD model relating the nifedipine plasma concentration to the change in systolic blood pressure was a dynamic binding PKPD model, as described by Shimada et al. (1996). Parameters used in the PD model are shown in **Table 2**. Simulated study design was matched to that reported for clinical studies, including age, proportion of females and fasted or fed state dosing. Ethnicity was also matched to the clinical study using the built in Simcyp Japanese and North European Caucasian populations. Where ethnicity of study subjects was not reported it was assumed based on the location of the approved study site or the country of residence of the study authors. Using the developed PBPK/PD models, concentration and response profiles were simulated for two different doses of the GITS formulation and compared with clinical data.

## **CASE STUDY 3:** *IN VITRO IN VIVO* **EXTRAPOLATION OF DIFFERENCES IN PD**

The third case study investigates the application of PBPK linked to a semi-mechanistic PD model to predict the response to a drug based on clinical data for a different drug that acts on the same target. Such models are useful for comparing several compounds that belong to the same therapeutic class. Semi-mechanistic PD models combine mechanistic aspects of the PD relationship with empirical features and are commonly used where the mechanism of drug action is not fully understood or when there is insufficient data available to develop a fully mechanistic model. The operational model of agonism (equation 1) was developed from receptor theory to describe *in vitro* pharmacology (Black and Leff, 1983) and has also previously been applied in PKPD modeling (Van der Graaf et al., 1997, 1999; Cox et al., 1998; Cleton et al., 1999, 2000; Zuideveld et al., 2004; Jonker et al., 2005).

$$\mathbf{E} = \frac{\mathbf{E}\_m \cdot \mathbf{r}^n \cdot [\mathbf{A}]^n}{(\mathbf{K}\_\mathbf{A} + [\mathbf{A}])^n + \mathbf{r}^n \cdot [\mathbf{A}]^n} \tag{1}$$

Mechanistic features are incorporated in terms of drug binding affinity (KA<sup>−</sup> which represents the binding affinity of drug A to the receptor) and intrinsic efficacy (ε; proportional to the transducer ratio τ), both drug-dependent parameters for

which information can be measured *in vitro*. The conversion of receptor activation to the PD response is described empirically by the system-dependent parameters Em, the maximum effect achievable in the system, n, the slope of the occupancy effect relationship and τ, which is related to the receptor concentration transduction properties of the tissue. System-dependent parameters are shared for drugs with the same mechanism of action in the same system. The operational model of agonism has previously been used to describe the PD effect of benzodiazepines in animal models (Cleton et al., 1999, 2000).

In this example, the operational model of agonism was used to describe the hypnotic effects of triazolam, as measured by change in beta-EEG amplitude, and to extrapolate the model to zolpidem by changing only the drug related parameters of the PBPK-PD model. The hypnotic effects of both zolpidem and triazolam are mediated via the same binding site on α1 subunit containing GABAA receptors. Triazolam and zolpidem were selected for two reasons. Firstly, the PD effect is related to the concentration of the parent compound only; zolpidem has no active metabolites, while, the active metabolite of triazolam is rapidly metabolized and thus does not contribute significantly to activity. Garzone and Kroboth (1989) Secondly, several clinical PK and PD studies have been reported for both compounds by the same group, which is important for the PD response since there is no standardization of the measurement and analysis of EEG recordings used as the PD effect measure, making it difficult to pool and compare data collected by different research groups. In this example, clinical data is used to establish the model for triazolam only, thus zolpidem is treated as a compound in pre-clinical development, with published clinical data used only to confirm the accuracy of the modeling approach.

Simulations of triazolam and zolpidem PK and PD were performed in virtual Caucasian healthy volunteers (HVs) and the Simcyp Triazolam (**Table 3**) and Zolpidem (**Table 4**) compound files with the first order absorption model, minimal PBPK distribution model and clearance described by enzyme kinetics. Unbound plasma concentration was used as the input to the PD model.

Equilibrium dissociation constant (KA) values for triazolam (1 nM) and zolpidem (53 nM) and relative intrinsic efficacy of the two compounds (Equation 1) were identified from published *in vitro* data (Garzone and Kroboth, 1989; Hadingham et al., 1993; Van der Graaf et al., 1999; Zuideveld et al., 2004). Since triazolam and zolpidem have comparable efficacy, the transduction ratio (τ) was assumed to be the same for both compounds. Simcyp Parameter Estimation module (using weighted least square objective function and Nelder–Mead optimization methods) was used to determine the values of τ, Em and n in addition to the effect compartment elimination rate (keo) to account for hysteresis in the response, using published data (Smith et al., 2001; Sancar et al., 2007). The weighted mean of the estimated values of each parameter was used in the simulations.

The quality of these parameter estimates to predict the PD response to triazolam was tested by the ability of the model to predict the PD response following interaction with the CYP3A inhibitor ketoconazole (von Moltke et al., 1996; Greenblatt et al., 1998). Thereafter, the dose of zolpidem predicted to produce the equivalent response to 0.25 mg oral triazolam (based on the maximal response and the area under the effect curve between 0 and 12 h) using the PBPK-PD model developed for triazolam and the KA and intrinsic efficacy for zolpidem, was identified using the Simcyp automated sensitivity analysis module for the dose range 0.01–1000 mg oral zolpidem. This was then compared with the clinically applicable dose for verification of the PBPK/PD model.

## **CASE STUDY 4: UNDERSTANDING THE COVARIATES DETERMINING PD VARIABILITY**

The fourth case study explores the ethnic differences in the QTc prolongation by quinidine using the target tissue (heart) drug concentrations. It has been suggested that in specific situations, PBPK/PD models are more likely to allow a better understanding of true PD variability versus variability resulting from drug disposition alone (Rostami-Hodjegan, 2013), which is usually reflected by plasma concentrations. Quinidine is known to cause lengthening of the QT interval in the electrocardiogram (ECG),with greater potential for QT prolongation in females (Benton et al., 2000; El-Eraky and Thomas, 2003; Shin et al., 2007). Ethnic differences in QT prolongation have also been demonstrated (Shin et al., 2007), with greater QT prolongation observed in Caucasian females than in Korean females, despite no significant differences in plasma concentrations. These differences in QT prolongation may be of significance clinically since lengthening of the QT interval corrected for heart rate (QTc) that is >500 ms is believed to be a contributory factor to the life-threatening side effect of Torsades de pointes observed with some drugs (Bednar et al., 2001).

Traditional PK/PD models linking plasma concentrations to QT changes in Caucasians and Koreans have reported a higher Emax (the maximum value of QTc changes) values in Caucasian females with similar EC50 (concentration of quinidine required to produce 50% of the maximum response) values in both ethnic groups, suggesting similar sensitivity to quinidine concentrations in the two groups (Shin et al., 2007). PBPK/PD modeling using free heart concentrations of quinidine that may be more relevant to the QT prolongation effect of the drug may have a greater potential to provide an understanding of the ethnic differences in the observed QTc changes.

Data from the study by Shin et al. (2007) were used to develop the PBPK/PD model, with virtual Caucasian HV and virtual Chinese HV (to represent Korean). The Simcyp compound file for quinidine (**Table 5**), a full PBPK distribution model with first order absorption and clearance of quinidine of 19.4 (CV 38%) L/h in Caucasians and 18.16L/h (34%) in Koreans was used. This PBPK model was verified by comparison of the plasma concentration versus time profile with clinical data. The Emax model used the measured mean baseline QTc of 443 ms for Koreans and 445 ms for Caucasians (Shin et al., 2007). Input to the PD model was predicted free heart concentrations and parameter estimation was used to estimate Emax and EC50. EC50 was used as a marker of sensitivity and compared in the two groups of virtual subjects.

## **RESULTS**

**CASE STUDY 1**

In general both PK and PD profiles were predicted successfully, as is evident from **Table 6** that summarizes the PK and PD parameters and **Figure 1** where the simulated data has been superimposed on observed data (Shimada et al., 1996; Kirchheiner et al., 2004). These models successfully simulated PK and PD profiles of metoprolol and support the potential for prediction of genetic differences in PD once the PKPD relationship is established in wild-type genotypes.

The simulated CL (Dose/AUC) of the UM group was found to be 16- and 2-fold higher than that of PM and EM groups, respectively, suggesting that UMs may not achieve adequate therapeutic response on a standard dose of 100 mg metoprolol. Simulated mean PD profiles showed that the area under the effect curve in PMs was 4-fold higher than that in UMs, and 2-fold higher than that in EMs. The simulated/observed ratios for the maximum reduction in heart rate and absolute area under effect curve are 0.94 and 1.2 for PMs, 1.0 and 0.94 for EMs, and 0.96 and 0.73 for UMs groups, respectively.

It is clear from these results that the status of CYP2D6 phenotype has an impact on the reduction in heart rate. PMs are of particular interest as the PD effect is higher and takes longer to return to the initial point. In comparison with EMs, and UMs, the longer action of metoprolol in PMs is a result of residence of drug in the body (see plasma concentration profile for PMs), which is caused by the lower clearance of metoprolol in PMs group. These differences indicate significant effects on metoprolol dosing in the corresponding groups of patients which could have been predicted *a priori.*

Simulation results showed consistency with clinical observations in terms of significant differences of metoprolol PK/PD profiles between PMs and UMs with a marginal change between EMs and UMs. UMs may not achieve optimal target concentrations of metoprolol, which can lead to a lower benefit from the standard 100 mg dose of the drug compared with PMs.

## **CASE STUDY 2**

Predicted PK and PD profiles for IR nifedipine in Japanese hypertensive subjects suggested that the model was successful in recovering the clinical data (Kikuchi et al., 1982). Comparison of the PK and PD parameters (Cmax and Rmax respectively) showed that the predicted Cmax/observed Cmax and predicted Rmax/observed Rmax are within the 2-fold acceptability criteria (**Table 7**).

Both the magnitude and sustained plateau (>24 h) of the PK and PD profiles were well captured for 60 mg nifedipine GITS formulation, with mean clinical data falling within the range of the mean values of simulated trials (**Figure 2**). The comparative PK and PD ratios in **Table 7** also confirm the successful prediction of the PK/PD profile of the 60 mg GITS formulation, for which a rich *in vitro* data set was available.

However, for a 30 mg multi-dose study of nifedipine GITS, visual inspection suggests that PK and PD is overpredicted (**Figure 2**), although the majority of observed values are within the 5 and 95% percentiles of the predicted profiles. Based on


values reported by Sharma et al. (2005).

#### **Table 6 | Observed vs. predicted "PRED" Metoprolol PK/PD parameters in healthy volunteers by CYP 2D6 metabolizer status.**

**and UMs (F).** Simulations are presented as the mean of 10 trials (bold black

**Table 7 | Comparison of the predicted and observed Cmax and maximum reduction in systolic blood pressure (Rmax) for the different nifedipine formulations and doses.**


*Data are reported as the mean* ± *SD (where reported). Observed values are from (a) Shimada et al. (1996), (b) Meredith and Elliott (2004), and (c) Brown and Toal (2008).*

**FIGURE 2 | Predicted and observed (A) plasma concentration profile and (B) change in systolic blood pressure after a single dose of nifedipine 60 mg GITS in North European hypertensive subjects (Meredith and Elliott, 2004).** Predicted and observed **(C)** plasma

concentration profile and **(D)** change in systolic blood pressure after the initial dose and daily dosing of nifedipine 30 mg GITS for 15 days in North European hypertensive subjects (Brown and Toal, 2008).

comparative ratios (**Table 7**) Cmax was marginally overestimated after the first dose, with a Cmax (predicted)/Cmax (observed) ratio of 2.25 and the other parameters within 2-fold of the observations. *in vitro* data for propagation and the prediction of the PK and PD profiles.

## **CASE STUDY 3**

It is notable that in this study no parameter fitting based on clinical data was used, with the aim of mimicking a situation in which prediction of the formulation effect is based on the use of

Fitted values of Em, τ, n, and keo were 20.6, 1.0, 0.93, and 2.0 respectively. The resulting model was able to predict the

**ketoconazole.** Observed data are from **(A,B)** von Moltke et al. (1996) and **(C,D)** Greenblatt et al. (2000) and the simulated study designs were matched to these studies. Predicted maximal observed response (PD

measure of PD response. **(G)** Predicted mean PD response to 10 mg oral zolpidem. Observed data points are from Greenblatt et al. (2000; closed

circles) and Greenblatt et al. (2006; open circles).

PD response to 0.125 and 0.25 mg triazolam with and without ketoconazole DDI reasonably well as seen in **Figure 3**, although the PD response was underestimated at the highest plasma concentrations of triazolam (**Figure 3D**).

A dose of ∼10–13 mg zolpidem was predicted to result in the same maximal response (Rmax) and area under the effect curve as a 0.25 mg dose of triazolam (**Figures 3E,F**). Visual inspection of the concentration – effect curve shows a good prediction of the maximum effect of zolpidem but the duration of the effect is overestimated (**Figure 3G**) compared with the clinical data (Greenblatt et al., 2000, 2006).

## **CASE STUDY 4**

The PBPK model predicted clinically observed plasma PK profiles of quinidine in Caucasian and Korean (represented by Chinese HVs) females (Shin et al., 2007) adequately as verified by visual predictive checks (**Figures 4A,B**). Simulations of free heart concentrations of quinidine over time for both groups are shown in **Figures 4C,D**.

Estimated Emax and EC50 values were 190.0 ms and 1.53 μM respectively in Caucasian females and 175.19 ms and 1.80 μM respectively in Korean females. Visual predictive checks of the simulations suggested that these PD models recovered the greater QTc prolongation observed clinically (Kim et al., 2007; Shin et al., 2007) in Caucasian females adequately (**Figure 5**).

The estimated sensitivity parameters (EC50) showed a Caucasian:Korean ratio of 0.85, indicating a greater sensitivity to heart quinidine concentrations in Caucasian females. This suggests that a standard dose of quinidine has the potential to produce QTc in more Caucasian females than in Korean females because of this difference in sensitivity.

## **DISCUSSION**

Recent advances in IVIVE coupled PBPK models have facilitated informed covariate recognition of the observed PK variability. Further, these models allow connecting response to the unbound drug concentrations at the site of action which in turn improves our ability to link the concentration-response relationships

**Caucasian females (B).** Solid lines represent mean values; dotted lines represent the upper and lower confidence intervals and solid circles represent observed data (Shin et al., 2007). Predicted free heart concentrations in Asian females **(C)** and Caucasian females **(D)**. Solid lines represent mean values; dotted lines represent the upper and lower confidence intervals.

beyond merely relying on the plasma concentration as the response's driving force. Connecting PD models, even empirical or semi-mechanistic ones, to PBPK models is the natural progression after developing predictive PK models.

To illustrate the added value in utilizing PBPK/PD models four case studies are presented here. These have demonstrated that linking a PD model to a PBPK model allows the prediction of the effects of a change in metabolizing enzyme phenotype, drug formulation, drug receptor binding, or ethnic differences in sensitivity to the drug on the PD response through propagation of the change in PK. Such models may assume that the

concentration-response relationship remains unchanged when the PK changes occur. It is recognized that this is not always the case, and any mismatch between predictions and observations may provide additional information about the mechanism of action of a drug and covariates relevant to the PD responses. Such factors can then be investigated and built into the models.

The first case study illustrates the potential for prediction of genetic differences in PD once the PKPD relationship is established in wild-type genotypes. Although population pharmacokinetic (POPPK) studies have been valuable in informing investigators of PKPD differences associated with different phenotypes, these studies need to be powered adequately to recognize such differences. Clinical trial simulations similar to the one shown in this case study can also be used to investigate the design of studies and their power to ensure that less frequently occurring phenotypes that are predicted to be relevant to dose evaluation are included. It is noteworthy that such simulations can be used in lieu of clinical studies to inform drug labels. (Janssen Therapeutics, 2013).

The second case study demonstrates that integration of a PBPK model that accounts for differences in formulation effects with a dynamic PKPD binding model predicts differences in clinical PD observations reasonably well, a feature that is very challenging to implement in classical compartmental approach. Prediction of the formulation effect based on *in vitro* data is propagated to the prediction of the PK and PD profiles, without the use of parameter fitting to observed data. The marginal over-prediction of the plasma profile can, in a large part, explain the over-prediction of the PD response that was observed in the study by Brown and Toal (2008). Overestimation of the plasma profile and PD response to 30 mg nifedipine GITS in this study may relate to differences in the dissolution profile of the batch of 30 mg GITS tablets used in the clinical study (which was not reported), compared to the dissolution profile used for simulation. In the absence of a published dissolution profile for the 30 mg GITS tablets, a dissolution profile proportional to the 60 mg GITS tablets were assumed for the simulations. It might have been expected that a PD model for the hypotensive response to nifedipine developed for and IR would underestimate the response to CR of nifedipine. Rate of increase in the plasma nifedipine concentration has been shown to influence the haemodynamic response, with a more rapid increase associated with increased sympathetic nervous system activation, increased heart rate and a diminished reduction in blood pressure response (Kleinbloesem et al., 1987; Meredith and Elliott, 2004; Brown and Toal, 2008). This may not have been the case since no increase in heart rate in hypertensive subjects was observed in the study from which Shimada et al. (1996) took the nifedipine PK and PD data to develop the PKPD model that was used in this case study (Kikuchi et al., 1982).

In the third example, by changing only the PBPK model input and the KA value for zolpidem, a good estimate of the dose requirement of zolpidem was obtained in HVs. The dose estimated using this PBPK/PD model is in agreement with the recommended dose of 10 mg zolpidem in adults. However, when compared to clinical data (Greenblatt et al., 2000, 2006), although the maximal response to zolpidem was predicted reasonably well, the duration of the response was overestimated. A possible explanation could be the clockwise hysteresis observed in the clinical data for zolpidem, suggesting acute tolerance effects to zolpidem, as has previously been proposed (de Haas et al., 2010). However, clockwise hysteresis is not observed for triazolam, suggesting that differences in the mechanism of action between zolpidem and triazolam may not be fully accommodated by the model. This example demonstrates the application of a combined PBPK and semi-mechanistic PD model in predicting the response to a compound based on clinical data for a different compound that acts at the same target.

Results in the fourth example demonstrate that the PBPK/PD model that used unbound heart concentrations of quinidine to drive the changes in QTc prolongation was effective in recovering the clinically observed ethnic difference in QTc prolongation. This model, which was of greater physiological relevance than previously published models, enabled us to gain a plausible mechanistic explanation for the observed ethnic differences in QTc prolongation, despite the similarities in the measured plasma concentrations in the two population groups. The higher EC50 in Asian females illustrate that this group is less sensitive to the QTc prolongation effects of quinidine and require higher free concentrations of quinidine at the target site to produce an equivalent change in QTc prolongation. Further studies to elucidate the mechanistic basis for the differences in sensitivities and also to investigate the potential contribution of 3-hydroxy quinidine (a primary metabolite that may contribute to pharmacological activity) are warranted.

Knowledge of inter-patient variability in response to drugs is crucial during drug development and clinical practice. PBPK/PD models such as the ones presented above provide a seamless framework to assess the propagation of key PK variables resulting from differences in physiology, genetics, demographics, concurrent medications, different formulations, etc. through to PD effects. Predicting such variability using the 'bottom up' approach prior to planning clinical trials enables researchers to optimize study design and predict results that are likely to be more reflective of the general population using the drug. Furthermore, when measured plasma concentrations cannot be reliably correlated with PD effects, PBPK/PD offer a valuable alternative to traditional compartmental modeling.

## **AUTHOR CONTRIBUTIONS**

Manoranjenni Chetty, Rachel H. Rose, Khaled Abduljalil, Theresa Cain, Gaohua Lu, and Nikunjkumar Patel designed and performed the research and analyzed the data. Manoranjenni Chetty, Rachel H. Rose, Khaled Abduljalil, Masoud Jamei, and Amin Rostami-Hodjegan wrote the manuscript.

## **ACKNOWLEDGMENTS**

This work was funded by Simcyp Limited (a Certara Company). The Simcyp Simulator is freely available, following completion of the training workshop, to approved members of academic institutions and other not-for-profit organizations for research and teaching purposes. The help of Eleanor Savill in preparing the manuscript is appreciated.

## **REFERENCES**


Zuideveld, K. P.,Van Der Graaf, P. H., Newgreen, D., Thurlow, R., Petty, N., Jordan, P., et al. (2004). Mechanism-based pharmacokinetic-pharmacodynamic modeling of 5-HT1A receptor agonists: estimation of in vivo affinity and intrinsic efficacy on body temperature in rats. *J. Pharmacol. Exp. Ther.* 308, 1012–1020. doi: 10.1124/jpet.103.059030

**Conflict of Interest Statement:** Manoranjenni Chetty, Rachel H. Rose, Gaohua Lu, Theresa Cain, Khaled Abduljalil, Nikunjkumar Patel and Masoud Jamei are employees of Simcyp Limited (a Certara company). Amin Rostami-Hodjegan is an employee of the University of Manchester and part-time secondee to Simcyp Limited (a Certara Company).

*Received: 18 August 2014; accepted: 04 November 2014; published online: 26 November 2014.*

*Citation: Chetty M, Rose RH, Abduljalil K, Patel N, Lu G, Cain T, Jamei M and Rostami-Hodjegan A (2014) Applications of linking PBPK and PD models to predict the impact of genotypic variability, formulation differences, differences in target binding capacity and target site drug concentrations on drug responses and variability. Front. Pharmacol. 5:258. doi: 10.3389/fphar.2014.00258*

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

*Copyright © 2014 Chetty, Rose, Abduljalil, Patel, Lu, Cain, Jamei and Rostami-Hodjegan. 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.*

## Virtual Systems Pharmacology (ViSP) software for simulation from mechanistic systems-level models

#### *Sergey Ermakov1 \*, Peter Forster 2, Jyotsna Pagidala3, Marko Miladinov3, Albert Wang3, Rebecca Baillie4, Derek Bartlett 4, Mike Reed4 and Tarek A. Leil <sup>1</sup>*

*<sup>1</sup> Exploratory Clinical and Translational Research, Bristol-Myers Squibb, Princeton, NJ, USA*

*<sup>2</sup> Forster Solutions, LLC, Wilmington, DE, USA*

*<sup>3</sup> Research IT and Automation, Bristol-Myers Squibb, Princeton, NJ, USA*

*<sup>4</sup> Rosa & Co. LLC, San Carlos, CA, USA*

#### *Edited by:*

*Chiranjib Chakraborty, Galgotias University, India*

#### *Reviewed by:*

*Emilio Benfenati, Istituto Mario Negri, Italy Gian Marco Leggio, University of Catania, Italy*

#### *\*Correspondence:*

*Sergey Ermakov, Bristol-Myers Squibb Company, PO Box 4000, Princeton, NJ 08540-4000, USA e-mail: sergey.ermakov@bms.com* Multiple software programs are available for designing and running large scale system-level pharmacology models used in the drug development process. Depending on the problem, scientists may be forced to use several modeling tools that could increase model development time, IT costs and so on. Therefore, it is desirable to have a single platform that allows setting up and running large-scale simulations for the models that have been developed with different modeling tools. We developed a workflow and a software platform in which a model file is compiled into a self-contained executable that is no longer dependent on the software that was used to create the model. At the same time the full model specifics is preserved by presenting all model parameters as input parameters for the executable. This platform was implemented as a model agnostic, therapeutic area agnostic and web-based application with a database back-end that can be used to configure, manage and execute large-scale simulations for multiple models by multiple users. The user interface is designed to be easily configurable to reflect the specifics of the model and the user's particular needs and the back-end database has been implemented to store and manage all aspects of the systems, such as Models, Virtual Patients, User Interface Settings, and Results. The platform can be adapted and deployed on an existing cluster or cloud computing environment. Its use was demonstrated with a metabolic disease systems pharmacology model that simulates the effects of two antidiabetic drugs, metformin and fasiglifam, in type 2 diabetes mellitus patients.

**Keywords: quantitative systems pharmacology, system-level mechanistic models, virtual patient, model development software, simulation experiment, metabolic diseases model**

## **INTRODUCTION**

In recent years pharmaceutical R&D has seen an increase in the development and application of mechanistic, systems-level models to inform decision making. These models are better at describing the disease biology and drug pharmacology than the more traditional and empirical pharmacokinetic/pharmacodynamic (PK/PD) models (Lalonde et al., 2007; Milligan et al., 2013; Visser et al., 2013). They are typically called quantitative systems pharmacology (QSP) models to distinguish them from other systems biology models that do not incorporate drug pharmacology or pharmacokinetics and often do not account for the biology of disease and its progression. With a detailed representation of physiology and pharmacology QSP models include a significantly larger number of equations and parameters compared to what is

**Abbreviations:** ER, extended release; FFA, free fatty acids; G6-P, glucose-6 phosphate; GI, gastrointestinal; GNG, gluconeogenesis; GIP, gastric inhibitory peptide; GLP-1, glucagon like peptide 1; GPR40, G-protein coupled receptor 40; IIGU, insulin independent glucose utilization; MDSP, metabolic diseases systems pharmacology; PD, Pharmacodynamic; PK, pharmacokinetic; RE, rapidly equilibrating; RR, readily releasable; SE, slowly equilibrating; TAG, triacylglycerol; UI, User Interface; ViSP, virtual systems pharmacology.

normally seen in traditional PK/PD models (Mager et al., 2003; Danhof et al., 2007). Despite the challenges of accurately determining all of the model parameters, QSP models can nevertheless be very informative by allowing the generation of quantitative hypotheses about the efficacy and/or safety of drugs prior to testing them in humans, or when testing in new patient populations (De Graaf et al., 2009; Kuepfer et al., 2012). Examples of mechanistic system-level models include, the first attempt to mathematically model the circulatory system in the human body (Guyton et al., 1972) and more recently, the HumMod model (Hester et al., 2011), and models of glucose homeostasis (Schaller et al., 2013), rheumatoid arthritis (Rullmann et al., 2005), hypertension (Hallow et al., 2014), and drug induced liver injury (DILI) (Shoda et al., 2014). A recent review by Schmidt et al. (2013) describes the process of how these models can be built and used.

With our incomplete knowledge of disease biology, QSP models can be used to make and test assumptions about the intrinsic variability in biological pathways. Because of the deterministic nature of the current approaches to the development of QSP models, it is expected that one set of initial conditions will normally produce a single set of outcomes with a unique solution trajectory from its initial state to the final one. In order to represent variability observed in a given population, multiple simulations with different initial conditions must be generated, each simulation implementing a virtual experiment or a virtual patient. Combined, the results from a sufficient number of simulations will provide estimates on the expected degree of variability in a population of patients. However, this comes at a price, by involving a large number of simulations and a large number of varied parameters the modeling process becomes computationally intensive. This challenge can be efficiently tackled only by employing the state of the art high-performance computing technology.

Likewise, the creation of large QSP models is not a trivial task; it requires the use of sophisticated and specialized software applications. Available tools range from complex sets of distributed software packages connected through a common portal to smaller yet versatile software programs capable of producing detailed mechanistic models. Examples of the former include the Garuda Alliance (Ghosh et al., 2011) and Physiome project (Thomas et al., 2008; Randall Thomas, 2009), while examples of the latter are JDesigner (Vallabhajosyula and Sauro, 2007), Entelos PhysioLab (Shoda et al., 2010), Mathworks SimBiology (MathWorks)1, Bayer's PK-Sim and MoBi (Eissing et al., 2011), and ISB's DBSolve Optimum (Gizzatkulov et al., 2010). For a given problem the choice of a proper modeling tool could become a difficult task by itself. In addition to purely scientific considerations dictated by the scope of the model, the software should meet multiple criteria to be considered optimal: an intuitive user interface, numerous differential equation solvers and library functions, a convenient way of storing and handling large number of parameters, ease of setting up multiple simulations and executing them in parallel, multiple-format import-export capabilities, reasonable cost and technical support, and an existing base of trained users. In this paper we present a simple and userfriendly Virtual Systems Pharmacology (ViSP) platform designed to quickly set up, run, and handle multiple simulation tasks in a flexible and scalable hardware/software environment. The platform is neither model nor software specific and can utilize existing cluster or cloud computing infrastructure for large-scale simulations. The ViSP platform was successfully used with a Metabolic Diseases Systems Pharmacology (MDSP) model to simulate multiple antidiabetic therapies in healthy and Type 2 diabetes mellitus (T2DM) patients.

## **METHODS**

In order to create flexible and versatile QSP software for setting up and running large scale simulations the following requirements were formulated:


The main innovation that enables the first three requirements is to separate the process of constructing a QSP model from the process of running simulations originated by that model. During the design step the model is created and then saved in what could be software proprietary file format. Afterwards it is the operating system (OS) that runs simulations, a process that could be implemented as software independent. One way to achieve this is to compile the model file(s) and convert into an executable code in which high-level model instructions are translated into low-level machine commands. Features proprietary to the modeling software will be removed during this translation and the executable will rely only on OS instructions. In order to keep the simulation process maximally flexible all model parameters need to be external to the executable file (i.e., their values not being hardcoded into the model). This means all model parameters should be presented as input parameters. Once input values are provided, together with the executable they will define a unique simulation task. Multiple simulation instances can be created by combining different input parameter values with executables, which can be run on a grid of processors in a cluster, in a cloud, or in a mixed environment (see **Figure 1**). As soon as the executable is compiled for the OS which runs the hardware, all computational resources integrated by the grid management software will be available for simulations.

Another important concept of the QSP software architecture is making it a web-based client-server application with a relational database back-end. With many advantages emphasized below it helps to address the requirements formulated in points 4 through 7. For example, a web-based platform is inexpensive to deploy and maintain, since it does not require installing software on every computer and web browsers are now omnipresent. Modifications and upgrades to the software could be done on the server side with no user intervention, thus reducing IT cost and time. Webbased software is easily accessible over networks inside or outside the company by multiple users whose access privileges can be regulated by IT departments. More advantages offered by webbased architecture come from the rich selection of software tools available for UI, front, and back end programming. Additionally, the database allows for a robust, reliable, dynamic, structured and complex relational data model where items such as models, users, virtual patients, project information, UI configuration settings, simulation inputs and results can be stored, managed and displayed by the UI components or other server-side modules.

Since QSP software is required to be capable of handling multiple models, its UI needs to be dynamic and configurable. This

<sup>1</sup>http://www*.*mathworks*.*com/products/simbiology/.

means the UI should expose model parameters in the way specific to each particular model. Also, some parameters should be immediately available via UI while other should be hidden. The latter could be necessary because physiological models often have hundreds of parameters, of which only a relatively small subset may be of interest for performing simulations. Such flexibility could be achieved by dividing all model parameters into meaningful groups, e.g., parameters that describe caloric value and composition of meals consumed by patients, or parameters specifying drug regimen, and so on (see **Figure 1**). Then, only the groups that are of interest will be selected during the UI configuration process; they will show up as sub-sections in the UI with specific parameters inside. An example of such an interface is given in the Results section.

Grouping provides additional benefits for handling and storing parameter values in a structured way. For instance, the same parameter group may get assigned different value sets corresponding to different individuals, here called virtual patients. Similar manipulations can be done with groups describing therapies, meals, and so on. Each set of values can be given a name and stored in the database with options for search and reuse. Once a sufficient number of such value sets is accumulated in the database, the end user's task of setting up simulations will be reduced to simply finding and selecting appropriate value sets. Again, the use of relational database enables and empowers this process, making it another key concept implemented in the QSP software architecture. In addition to parameter value sets, practically all other information about the model, UI configuration, and simulation results is stored in the database that is searchable and that preserves the relationships between these pieces of information.

The software architecture with the features discussed above was implemented in the ViSP platform, a flexible tool for setting up and running large-scale QSP simulation tasks. Together with attempts to make simulation process less dependent on specific modeling tools and proprietary model formats we tried to establish a more universal workflow for simulations (**Figure 2**). In this article we describe the implementation of this simulation workflow using the MDSP model designed to study the effects of anti-diabetic drugs. The MDSP model itself was generated using JDesigner software (Sauro et al., 2003) and then exported as an m-file used by Matlab® software by Mathworks (MathWorks) (**Figure 2**). Afterwards, the "main" program controlling the simulation process was added to the model, and the code was modified such that all model parameters became input parameters as per the requirements described above. Using the Matlab Compiler Toolkit®, the code was then compiled into a standalone executable. The executable was subsequently uploaded into the ViSP database while a designated Power User (experienced user with highest privileges) configured the user interface (UI) to reflect the specifics of the model. Once the UI is in place all users may set up, submit and run simulations on any number of nodes, for instance via Amazon Web Services™. When simulations are complete the results are saved to the database in a text format for further processing. The same process can be repeated with any model as soon as it can be converted into an executable file.

ViSP specific details along with its application to the MDSP model are illustrated in the Results section. The MDSP model high-level organization is outlined in the next sub-section, however a complete description is out of this paper's scope. The mathematical basis of the MDSP model is provided in the Appendix.

model development software. It is then converted into an executable file with the help of 3-d party software. The Power User also configures the user interface (UI) that is specific to the model by using ViSP. Executable, text file with full set of parameter baseline values, and Virtual Patients are uploaded by ViSP into the database (DB). Power and Regular Users can

## grid (cluster, cloud) and retrieve the results after simulations are completed. The latter will be available for analysis to users through the ViSP UI. The section of the workflow that deals with model conversion into an executable file (inside gray rectangle) is currently implemented outside of ViSP.

## **MDSP MODEL**

The MDSP mathematical model was developed to mechanistically describe the basic physiological and pathophysiological processes involved in T2DM. It represents essential systems and mechanisms regulating glucose and lipid metabolism and describes pathophysiological changes related to T2DM together with the PKPD effects for several classes of antidiabetic drugs (for recent review of mathematical models of diabetes please see Ajmera et al., 2013). The core of the model simulates intake and processing of nutrients, and their distribution and utilization by different body tissues and organs as schematically represented by a block diagram on **Figure 3**. The nutrients enter in the form of meals (up to three per day) with a specific percentage of carbohydrates, fats and proteins and with a given caloric content (all these can be modified through the ViSP UI). From the GI tract the nutrients are absorbed into the bloodstream and the model further tracks glucose and lipid metabolism by the brain, liver, muscle and adipose tissues (**Figure 3A**) (Zierler, 1999). In

the liver, glucose is phosphorylated to become glucose-6 phosphate (G6-P) to be afterwards converted into glycogen (Agius, 2008). Both reactions have their counterparts working in the opposite direction such that the net glucose flux into/out of the liver maintains plasma glucose concentration within a specific range depending on the feeding condition. The liver produces glucose from the three-carbon substrates through the process of gluconeogenesis (Radziuk and Pye, 2001). The above processes are subject to insulin and glucagon regulation, and are disrupted in T2DM, resulting in increased hepatic glucose output in the postprandial and fasting states compared to healthy subjects.

Muscle also stores glucose in the form of glycogen but unlike the liver it does not use glycogen to release glucose back into circulation (Laurent et al., 2000). Insulin regulates muscle glucose uptake, storage and utilization, and T2DM have decreased sensitivity to these insulin effects. In regard to the other tissues involved in glucose metabolism, adipose tissue uptakes glucose for either storage or oxidation (**Figure 3A**), brain consumes glucose in a constant, insulin independent fashion, while kidneys normally reabsorb all filtered glucose unless its concentration exceeds a threshold value (Rave et al., 2006; Marsenic, 2009). Lipids in the model are represented as pools of triacylglycerols (TAGs) and free fatty acids (FFA) stored and transported between several compartments (**Figure 3A**).

As mentioned earlier, nutrient disposal by tissues is tightly regulated by multiple hormones, insulin being the most important one. In the MDSP model the secretion and action of insulin is described by a multistep process (**Figure 3B**) that is coupled to plasma glucose concentration. Other factors included in the model that affect levels of insulin are beta-cell mass and betacell function (Bouwens and Rooman, 2005), activation of cAMP pathways (Fridlyand et al., 2007) and activation of Ca+ pathways (Bertuzzi et al., 2007). Insulin is degraded primarily by the liver and partially by peripheral tissues, with C-peptide being an important by-product and biomarker of insulin secretion tracked by the model. Insulin's counterpart glucagon is described by a simpler two-compartment dynamic model (**Figure 3B**). Its regulatory effects are implemented as stimulating gluconeogenesis and glycogenolysis in the liver. Two other metabolic incretin hormones, glucagon like peptide-1 (GLP-1) and gastric inhibitory peptide (GIP) are also implemented in the model (**Figure 3B**), since they represent potential targets for therapies. GLP-1 affects glucose uptake and oxidation by adipose tissue and both hormones influence insulin and glucagon secretion in response to glucose.

The ViSP platform was used to calibrate (see Appendix) and then run the MDSP model to simulate the effects of meals, glucose and meal tolerance tests, and several antidiabetic drugs in different patient phenotypes. Examples of simulation results for two of such drugs are presented in the following Results sections. The first example illustrates the simulated effects of metformin, considered by many as a standard of care for T2DM patients, compared to literature data. The other example presents results with a relatively new class of drugs, GPR40 agonists (GPR40a), with simulations reproducing the effects of fasiglifam (TAK-875). The last example presents simulation results for metformin + TAK-875 combination therapy.

### **RESULTS**

#### **ViSP PLATFORM**

The ViSP software features several primary user-interface components, the first of which is the Explorer (see **Figure 4**, left side). It organizes user's models and data in a tree-like hierarchical structure in which top elements are projects. A ViSP project typically comprises all information related to simulation tasks that pertain to the specific research topic. Each project can contain one or several models, for instance, different versions created in the course of the model development. The model is represented by an executable file which is uploaded into the ViSP databases every time a new model is created. The executable is accompanied by a text file which contains a list of model input parameters and their baseline values. Every model can be associated with one or several user interfaces (UI) that are configured to fit particular project needs or user preferences. The next level down in the hierarchy comprises groups of parameters that are subsets of input parameters in the sense explained in the Methods. Each group may further contain multiple value sets, reflecting, for instance, settings for different drug regimens (**Figure 4**). All elements of the structure are stored in a database that facilitates handling of the model, data, and results.

Another important feature in ViSP is the Simulation Manager. It provides a means to customize the UI to the content of the model and then prepare and launch simulation tasks. The model UI can be configured by a Power User in a simple setup by creating sections that deal with particular aspects of the model. For example, in the section of the MDSP model specifying meal regimen, out of all parameters related to meals only the parameters defining the regimen are selected (see **Figure 4**, right side). Consequently, only these parameters will be presented in the UI through a series of controls. The Power User assigns meaningful text labels to these controls and specifies how they should be displayed in the UI, as a check box to turn a parameter on/off, as a text entry field, or as a dropdown selection. The configuration table helps to arrange the controls in a simple grid by specifying row and column numbers (see **Figure 5**). By creating various sections in this manner, the Power User has full control over which of the hundreds of parameters of the model to display, and how. The Power User can create several UIs targeting different groups of users which require particular aspects of the model to be exposed.

Once the sections of the UI have been configured, the Simulation Manager presents an option of selecting some or all of the Virtual Patients (VPs) known to the system for this model (see **Figure 6**). For convenience VPs are classified according to their phenotypes, thus facilitating proper VP selection required for simulations. Even though VPs come with all parameters defined, there is an option to change some of them if a user finds that necessary. After VP selection is complete, the Simulation Manager will create a set of single simulations for every combination of the settings and VPs. For example, if three VPs were selected and the parameters were set to apply one therapy, then three simulations will be generated and submitted, one for each patient. However, if two therapies were selected (the same drug with different dose, or two different drugs, etc.), then the Simulation Manager will generate six simulations accordingly. Simulation Manager also allows the therapies (dropdown selections) to be applied as "Combination" treatments, which in the above example means both therapies get applied to each patient, resulting in three simulations (see **Figure 6**).

The final settings that are specified through an additional UI window (not shown) are duration of simulation, output variables and time intervals between outputs. Once those are provided the

simulations are fully defined and can be submitted to the computational grid. The user will get notified by e-mail upon task submission and when simulations are completed. The results that are saved in a series of text files can be retrieved afterwards via the Results Manager for further analysis and graphical visualizations. ViSP itself is capable of generating graphical plots which can be viewed directly as part of the results.

ViSP's Administrator tool provides means to register and grant access only to users who are authorized to use the software and its data, thus preventing any proprietary information from disclosure. Additionally all users are divided into Power Users and Regular Users based on their privileges. As was described above Power Users are allowed to create and modify projects, configure model UIs, and set up and run simulations, while Regular Users can perform only the last two functions.

## **SIMULATIONS METFORMIN**

The pharmacokinetics of metformin was simulated by using a three-compartment model (**Figure 7A**) derived from the Pentikainen et al. (1979). The model was calibrated to fit PK characteristics for a 500 mg single dose (Pentikainen et al., 1979) and multiple 500 mg twice daily doses (Graham et al., 2011) obtained with healthy individuals. An adequate fit has been achieved for both data sets as evidenced by **Figure 8**. It was deemed acceptable to apply the same calibration for simulating metformin pharmacodynamic (PD) effects in T2DM patients. This assumption is supported by the data from the Tucker et al. (1981), which found little difference in metformin PK between healthy and T2DM individuals.

Metformin has multiple sites of action, including liver, muscle, adipose tissue, GI tract, and pancreas. Despite the fact that metformin is perhaps the most widely used antidiabetic therapy its exact mechanism of action remains unclear (Kirpichnikov et al., 2002). Among its reported primary PD effects are decreased hepatic glucose production (Stumvoll et al., 1995; Campbell et al., 1996), increased peripheral tissue sensitivity to insulin (Bailey and Turner, 1996), and increased glucose utilization. Other less commonly described effects include lowering FFA levels and increasing lipid oxidation (Perriello et al., 1994), increased glucose utilization by the GI tract, and a delayed, more distal GI glucose absorption (Bailey et al., 2008). There are different opinions on whether metformin directly affects β-cells (DeFronzo, 1999), however some evidence exists that it improves the function of β-cells (Patane et al., 2000; Bi et al., 2013) and their response

to glucose. In the MDSP model the above mentioned metformin PD effects were implemented as multipliers in the rate equations. Functionally they are expressed as Hill equations (Appendix, Equation A3) representing either metformin-mediated activation

parameter is associated with a particular control type (e.g., TEXTBOX) that is supplied with a text label. Each control location is defined by UI

or inhibition. A study on the short-term effects of metformin in T2DM patients (Eriksson et al., 2007) was selected for simulations demonstrating the model's ability to reproduce metformin therapeutic outcomes. In this study an escalating metformin dose (500 mg *qd* for 7 days followed by 500 mg *bid* for 7 days and then by 1000 mg *bid* for 14 days) was applied to a group of T2DM patients with fasting plasma glucose concentration between 7 and 12 mM. At the beginning of each subsequent dose an oral glucose tolerance test (OGTT) was performed to check the effects of the previous dose on glucose and other metabolic characteristics, (for further details see paper by Eriksson et al., 2007). In simulations the same treatment regimen was reproduced for a representative virtual patient that matched the study enrollment criteria including body weight, fasting plasma glucose (FPG), age, etc. **Table 1** provides a comparison between clinical and simulation data for key parameters, including FPG, area under the glucose concentration curve for OGTT, and percent change in fasting plasma insulin (FPI) concentration from day 0 before treatment and after 7, 14, and 28 days of metformin. Overall there is good agreement between simulations and data, with simulations slightly under-predicting the decrease in FPG especially at higher doses.

## **GPR40 AGONIST (TAK-875)**

the inset (front).

TAK-875 is a selective GPR40 agonist that improves glycemic control in T2DM patients by potentiating postprandial insulin secretion in a glucose dependent manner with a minimal risk for hypoglycemia (Kaku, 2013; Yabuki et al., 2013). A single dose TAK-875 PK study with healthy volunteers (Naik et al., 2012) and a multiple dose study with T2DM patients (Leifke et al., 2012) were used to establish and calibrate the PK section of the MDSP model (**Figure 7B**). An enterohepatic recirculation (EHRC) was included in order to better fit the TAK-875 clinical data (**Figure 9**).

if necessary. With the above configuration the section will look like on

The mechanism by which TAK-875 potentiates insulin secretion involves activation of GPR40 in pancreatic β-cells followed by a cascade of reactions increasing the levels of secondary intracellular messengers. This eventually results in increased Ca2<sup>+</sup> release that enhances the movement of insulin granules and their fusion with the plasma membrane, leading to subsequent insulin release (Burant, 2013). In the MDSP model all these events are simplified into one mechanism representing the net TAK-875 amplification of the Ca2<sup>+</sup> effect on insulin secretion. GPR40 is also expressed in enteroendocrine cells of the intestine, and it has been hypothesized that GPR40 activation may potentially lead to increased secretion of GLP-1 and GIP hormones (Luo et al., 2012; Mancini and Poitout, 2013). These pathways are represented in the model as hypotheses, so that their potential impact on efficacy could be evaluated. However, the results of the TAK-875 clinical study in T2DM patients did


not demonstrate increases in GIP or GLP-1 following an OGTT (Leifke et al., 2012). Therefore, the secretion of GIP and GLP-1 via the intestine was disabled for simulations with TAK-875. In choosing a representative VP for simulations, as in the case of metformin, we selected one with steady-state characteristics comparable to the mean values found in the study enrollment criteria.

to be run in simulations. Parameters for selected VPs can be modified,

**Figure 10** compares simulation results with clinical data from a multiple ascending dose study of TAK-875 in T2DM subjects that received either placebo or one of the 25, 50, 100, 200 mg daily doses (Leifke et al., 2012). Data shown in the figure illustrate the short-term TAK-875 effects on steady state responses (FPG concentration, **Figure 10A**), and dynamic responses (2 h post-OGTT glucose concentration **Figure 10B**) after 14 days with different levels of drug exposure. Simulations provide adequate predictions in both occasions although simulated glucose post-OGTT values seem to follow clinical data more closely.

#### **METFORMIN—TAK875 COMBO**

combination therapy.

Since metformin is used as the standard of care for treating hyperglycemia in T2DM patients, we repeated the above simulation of TAK-875 in combination with 500 mg metformin twice daily as a background therapy. Simulation results suggest that additional therapeutic benefits could be achieved by combination therapy by further lowering FPG and post-prandial glucose excursions (**Figure 11**). Interestingly, the effect on post-prandial glucose appears to be more pronounced, with the response at TAK-875 doses of 50 mg and higher approaching a plateau (**Figure 11B**). In contrast the decline in FPG over the same dose range

**FIGURE 8 | Metformin simulated steady-state concentration profile (line) is plotted against clinical data (diamonds) (Graham et al., 2011) for 500 mg bid dose regimen.** On the inset single 500 mg dose simulation results are compared with data from Pentikainen et al. (1979) (mean ± SE) for metformin peak concentration *C*max, time to peak *T*max, and area under the curve (AUC).

(**Figure 11A**) did not appear to have reached saturation.Without metformin (**Figure 10**) this plateau in OGTT response is observed in both clinical and simulation data but at higher (*>*100 mg) doses than with combination therapy.

## **DISCUSSION**

QSP models bring new insights into our understanding of the mechanism of action of drugs and they help in optimizing decision making in pharmaceutical R&D (Schmidt et al., 2013). However, multiple obstacles need to be overcome in order to increase recognition of the value of QSP within the pharmaceutical industry. Two challenges are worth mentioning in the context of this article. First, QSP models rely on large-scale simulations requiring high-performance parallel computing infrastructure that is expensive to run and maintain. Second, there is no industry standard software, such as NONMEM® for non-linear mixed effects PK/PD modeling that satisfies the diverse needs of the modeling community. Utilizing multiple tools increases the cost of model development and limits the exchange of models between scientists, thus creating additional barriers for model acceptance and application. A solution for the first challenge could be cloudbased computing, as the burden of creating and maintaining an up to date computational environment is outsourced to vendors of high-performance computing clusters (e.g., Amazon). By developing the ViSP platform we attempted to address the second challenge, i.e., making the simulation process less dependent on the modeling tools and creating a more universal workflow for simulations (**Figure 2**). The central idea behind the ViSP platform is to work with the model file once it satisfies two conditions; first, it is compiled into an executable file, and second, all model parameters are presented as input parameters. Combination of an executable binary file with an input text file fully defines a single simulation task that has the following benefits. On the one hand, it is no longer dependent on the file format or the specifics of the modeling tool that created the model. On the other hand,


**Table 1 | Effects of metformin in 28 day study in T2DM patients, comparison between clinical data by Eriksson et al. (geometric mean ± 95% conf. interval) and simulation.**

**FIGURE 9 | TAK-875 concentration profiles for 25, 50, 100, and 200 mg once daily doses administered to T2DM patients.** Clinical data points are mean values ±SD (Leifke et al., 2012) shown by black color markers connected with lines overlaid with simulated data shown by gray color markers only. Panel **(A)** presents results at day 1, panel **(B)** at day 14.

the model preserves all possibilities for its customization since all its parameters are available through the input file. Additionally, when launched it runs as a single computer process that provides flexibility in choosing the hardware (multicore processors, cluster, or cloud) that can be used for computations. The only requirement here is to use the proper compiler when creating the executable file.

Large-scale simulation tasks in which ViSP could be useful originate from numerous applications. We employed ViSP for calibrating the MDSP model and for simulating clinical studies. In

the first case, multiple virtual patients (VPs) have been created in which only the parameters of interest were changed. Then a series of simulations equal to the number of virtual patients has been run and the process was repeated until the desired model behavior was achieved. The latter meant checking that sets of output parameters lay within the observed ranges derived from clinical or preclinical data. A similar procedure was used to create new VPs representing different phenotypes. In the future we are planning to automate this process, when parameter variations and model response analysis will be done without user intervention. The process just described could be applied to perform sensitivity analysis, for instance, to search for the pathway that responds the most (or the least) to the drug, or to characterize the drug response based on patient phenotype. This process could also be used to model a clinical trial, when different cohorts of VPs that satisfy the enrollment criteria are simulated and their responses are analyzed to provide suggestions for patient stratification.

One aspect of simulation workflow that remains outside the capabilities of the ViSP platform is how to convert model files saved in proprietary formats into an executable code (see **Figure 2**). Currently there is no universal mechanism inside ViSP allowing this to be done with an arbitrary file format. The proprietary nature of the model files prevents seeing model details, such as equations and parameters, making compiling such files into an executable impossible. Normally the modeling software itself does not offer this option either. The solution, however, exists if the modeling tool allows the export of models into a file format that can be read by other software. One such format is Systems Biology Markup Language (SBML), a computer language that is gaining ground inside the Systems Biology community for saving and exchanging models between users. Currently several model development tools offer SBML export capabilities, among them JDesigner (part of the Systems Biology Workbench, SBW) (Sauro et al., 2003), SimBiology by Mathworks (MathWorks), CellDesigner by the Garuda Alliance (Kitano et al., 2005), DBSolve Optimum by ISB (Gizzatkulov et al., 2010) and others. Once saved in SBML format, a model file can be translated into a different computer language that afterwards can be compiled into an executable (see **Figure 2**). We utilized SBW capabilities to export an SBML file into a MathWorks Matlab® file that later was compiled into a binary executable file using the MathWorks compiler. Since the MDSP model was originally developed in JDesigner, which "natively" saves models in SBML (XML) format, there were no issues in exporting it to a Matlab file. However, if the SBML model file is produced by a different modeling tool, for example PhysioLab®, it may require some editing before saving it as Matlab code. Such modifications may be necessary since the level of SBML support varies in different modeling tools.

In conclusion, we developed a versatile web based software platform that provides capabilities for setting up and running massive simulation tasks originating from system-level mechanistic models. It is designed to conveniently handle diverse modeling projects with large number of parameters while being flexible with respect to model structure. Its utility was demonstrated with metabolic diseases model by simulating pharmacological effects of antidiabetic drugs, metformin and fasiglifam in healthy and diabetic individuals.

## **REFERENCES**


**Conflict of Interest Statement:** Sergey Ermakov, Jyotsna Pagidala, Marko Miladinov, Albert Wang, and Tarek A. Leil are employees of Bristol-Myers Squibb, P.Forster is an employee of Forster Solutions, LLC, Rebecca Baillie, Derek Bartlett, and Mike Reed are employees of Rosa and Co. LLC.

*Received: 18 August 2014; accepted: 30 September 2014; published online: 22 October 2014.*

*Citation: Ermakov S, Forster P, Pagidala J, Miladinov M, Wang A, Baillie R, Bartlett D, Reed M and Leil TA (2014) Virtual Systems Pharmacology (ViSP) software for simulation from mechanistic systems-level models. Front. Pharmacol. 5:232. doi: 10.3389/ fphar.2014.00232*

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

*Copyright © 2014 Ermakov, Forster, Pagidala, Miladinov, Wang, Baillie, Bartlett, Reed and Leil. 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.*

## **APPENDIX**

## **BASIC EQUATIONS**

Mathematically the MDSP model is organized as follows. Individual blocks *i* (*i* = 1*,..., n*) shown on **Figure 3** correspond to quantities of a material *Qi* (e.g., amount of nutrient, number of cells, etc.) distributed within a particular compartment. Blocks may also stand for a state or a signal attributed to an organ or body system, like e.g., vagal stimulation of the stomach. The quantities are expressed as time-dependent state variables described by ordinary differential equations (ODEs) (Equation A1) where the right side term *Fi,<sup>j</sup>* signifies the flux of *Qj* to compartment *i* (shown by arrows on the diagram). These fluxes may depend on multiple variables, for instance, concentration of hormones and/or drugs. Most of them follow similar functional form (Equation A2) in which basal rate *ki,<sup>j</sup>* is modified by activation, represented by the term *αi,<sup>j</sup>* and/or by inhibition, represented by the term *βi,j*. Activation and inhibition terms are given by either the Hill-type equation (Equation A3) or by a sigmoid shape function (Equation A4) where coefficients *γi,<sup>j</sup>* and *σi,<sup>j</sup>* regulate the slope and the shift of the curve. Here functions are shown for activation term *αi,j*, they are identical for *βi,j*. Parameters *αmax i,j* may assume only non-negative values. In the case of inhibition, *βmax <sup>i</sup>,<sup>j</sup>* are further restricted to be no greater than 1. The variable *Pj* may stand for either the quantity *Qj* or its concentration *Cj* in case when *Qj* represents the amount of substance. Parameters *ni,j*, *γi,j*, *σi,<sup>j</sup>* are constant values. In each particular instance the choice of the function, either Equation A3 or Equation A4, is dictated by a prior knowledge about the mechanism of action. Alternatively it is selected based on which function provides better fitting to available data. For details refer also to **Figure 3**.

$$\frac{dQ\_i}{dt} = \sum\_j F\_{i,j} \tag{A1}$$

$$F\_{i,j} = k\_{i,j} \cdot Q\_j \cdot \left(1 + \alpha\_{i,j}\right) \left(1 - \beta\_{i,j}\right) \tag{A2}$$

$$\alpha\_{i,j} = \alpha\_{i,j}^{\max} \cdot \frac{P\_j^{n\_{i,j}}}{P\_j^{n\_{i,j}} + K\_j^{n\_{i,j}}} \tag{A3}$$

$$a\_{i,j} = \alpha\_{\max} \cdot \frac{1}{2} \left( 1 + \tanh\left(\wp\_{i,j}\left(P\_j - \sigma\_{i,j}\right)\right) \right) \qquad (\text{A4})$$

In addition to ODEs, the MDSP model contains a number of algebraic equations that calculate additional quantities and parameters, for instance, insulin sensitivity index QUICKI (Katz et al., 2000), homeostatic model assessment index HOMA (Wallace et al., 2004) and others. Overall the model comprises more than 100 ODEs and more than 50 algebraic equations resulting in a large number of parameters (more than 800) associated with them.

## **MODEL CALIBRATION**

As explained in the Methods section, parameters in the model could be approximately divided into those describing the characteristics of an individual subject [virtual patient (VP)], parameters that represent drugs, parameters for clinical interventions, e.g., OGTT, parameters that describe meal regimen, and the like. All of them are input data defining a particular simulation. Whenever possible the values for parameters are taken from the literature, however, when the values are not available they are estimated by matching model behavior to known clinical and preclinical data. For a valid VP, values for plasma concentration of glucose, insulin, C-protein, glucagon, GLP-1, GIP, FFA, TAG obtained after simulated overnight fast are required to match clinical data characteristic of the phenotype this VP represents (healthy or T2DM). In addition a valid VP should correctly reproduce simulated interventions such as a standard OGTT and meal tolerance test (MTT) with glucose and insulin concentration profiles being similar to the ones observed clinically. Since multiple sets of parameters could potentially match the same data, these sets constitute alternative VPs that reflect clinical variability. Parameters defining VPs can be varied intentionally to produce subjects representing different behaviors (phenotypes). For the MDSP model we have created a number of VPs that we classified as belonging to the following 5 phenotypes: normal healthy, obese non-diabetic, type 2 diabetic patients with mild, moderate, and severe degree of the disease. Each VP is defined by more than 800 parameters and ViSP provides the convenience of storing and handling them in a structured way.

In general the simulations were performed as follows. At first all VPs are simulated for a period of time until they reach a steady state with a meal regimen that matches VP energy expenditure. The latter is calculated based on VP age, height, body weight, and activity level. After the steady state is achieved the simulation of interest is initiated by applying additional sets of parameters corresponding to drugs and interventions that reproduce the conditions of the clinical trial or other experiments of interest.

## ADVANTAGES OF PUBLISHING IN FRONTIERS

FAST PUBLICATION Average 90 days from submission to publication

## COLLABORATIVE PEER-REVIEW

Designed to be rigorous – yet also collaborative, fair and constructive

RESEARCH NETWORK Our network increases readership for your article

## OPEN ACCESS

Articles are free to read, for greatest visibility

## TRANSPARENT

Editors and reviewers acknowledged by name on published articles

## GLOBAL SPREAD Six million monthly page views worldwide

## COPYRIGHT TO AUTHORS

No limit to article distribution and re-use

IMPACT METRICS Advanced metrics track your article's impact

## SUPPORT By our Swiss-based editorial team

EPFL Innovation Park · Building I · 1015 Lausanne · Switzerland T +41 21 510 17 00 · info@frontiersin.org · frontiersin.org