gPKPDSim: a SimBiology®-based GUI application for PKPD modeling in drug development

Modeling and simulation (M&S) is increasingly used in drug development to characterize pharmacokinetic-pharmacodynamic (PKPD) relationships and support various efforts such as target feasibility assessment, molecule selection, human PK projection, and preclinical and clinical dose and schedule determination. While model development typically require mathematical modeling expertise, model exploration and simulations could in many cases be performed by scientists in various disciplines to support the design, analysis and interpretation of experimental studies. To this end, we have developed a versatile graphical user interface (GUI) application to enable easy use of any model constructed in SimBiology® to execute various common PKPD analyses. The MATLAB®-based GUI application, called gPKPDSim, has a single screen interface and provides functionalities including simulation, data fitting (parameter estimation), population simulation (exploring the impact of parameter variability on the outputs of interest), and non-compartmental PK analysis. Further, gPKPDSim is a user-friendly tool with capabilities including interactive visualization, exporting of results and generation of presentation-ready figures. gPKPDSim was designed primarily for use in preclinical and translational drug development, although broader applications exist. gPKPDSim is a MATLAB®-based open-source application and is publicly available to download from MATLAB® Central™. We illustrate the use and features of gPKPDSim using multiple PKPD models to demonstrate the wide applications of this tool in pharmaceutical sciences. Overall, gPKPDSim provides an integrated, multi-purpose user-friendly GUI application to enable efficient use of PKPD models by scientists from various disciplines, regardless of their modeling expertise. Electronic supplementary material The online version of this article (10.1007/s10928-017-9562-9) contains supplementary material, which is available to authorized users.


Introduction
Characterization of drug pharmacokinetics (PK) and its relationship with pharmacodynamic (PD), efficacy and safety endpoints are critical in pharmaceutical research and development and are used to inform efforts such as target and molecule selection, human PK projection, and clinical dose and schedule determination. Quantitative analysis and modeling and simulation (M&S) tools are frequently employed to characterize and predict PKPD relationships, understand drug mechanism of action, support optimal design of studies and reduce animal usage in preclinical studies, and guide clinical translation and clinical study design. Commonly used M&S approaches include noncompartmental analysis and compartmental PK models with empirical or mechanistic pharmacodynamics (PD) and toxicology and/or efficacy models. Further, exploring the variability within these models is a critical component of model-based analyses.
Typically, these tasks are performed by scientists with training and expertise in PKPD modeling and simulation. However, while model development might require such expertise, certain model-based exploration and applications could be conducted by scientists from other disciplines (e.g., pharmacology, toxicology, biomarker discovery & development, bioanalytical sciences, and clinical sciences) if provided with validated models in user-friendly, graphical-user-interface (GUI) based tools. These scientists could ideally directly utilize such modeling tools in their study design, data analysis, interpretation and decisionmakings.
A variety of software and applications are available to perform PKPD analysis, modeling, and simulation. Some of the most commonly used PKPD modeling tools include Phoenix Ò WinNonlin Ò , NONMEM Ò , Adapt V, SAAM II, R and MATLAB Ò SimBiology Ò . These modeling tools can be script-based or have a graphical user interface to help construct and simulate models. Pharmacometricians and PKPD M&S scientists are typically trained in one or more of these tools. Although some of the tools are more userfriendly than others, these are rarely used by non-specialists even when the mathematical models are already available. Hence a simpler multi-purpose tool with a userfriendly interface can promote the broader use of PKPD models in various stages of drug development, by scientists from diverse disciplines.
Wojciechowski et al. [1] have developed an interactive application using Shiny for the R programming language that can help with sharing and communication of pharmacometric model results. Shiny allows for customization of the application's user-interface to provide an environment for displaying user-input controls and simulation output that can be simultaneously updated with changing input. Ermakov et al. [2] have developed the Virtual Systems Pharmacology software platform for simulations of models in a software free environment. The platformenabled configuration of models developed in any software of choice for simulations on a web based interface with a back-end database capability. The platform, however, is not publicly available. Germani et al. [3] introduced the A4S Simulator, developed in MATLAB Ò , that uses a library of models and provides different dosing options; however, features such as non-compartmental analysis, interactive visualization and project management are limited in this program. Krause et al. [4] demonstrated the visualization of pharmacometric models using the Berkeley Madonna software with a particular focus on interactive sessions to facilitate communications with non-modelers.
We have developed a GUI-based application called gPKPDSim (Genentech PKPD Simulator) to broaden utilization of preclinical and translational PKPD models in drug development. The tool is made available with this manuscript and allows deployment of models developed by expert modelers to a broader group of scientists and researchers for hands-on use. Any model developed within SimBiology Ò can be configured into gPKPDSim. The application and features of gPKPDSim are demonstrated using a library of standard PKPD models for typical modelrelated drug development tasks.

Overview of software development
The gPKPDSim application was developed using MATLAB Ò , SimBiology Ò and the Statistics and Machine Learning Toolbox and the user requires a MATLAB Ò license with the mentioned toolboxes to be able to launch the application. If the Optimization or Global Optimization Toolbox is installed, gPKPDSim can also use the parameter estimation methods from either of those toolboxes. The default choice for the optimization method is nlinfit (which is suitable for nonlinear least-square problems), however, the modeler can change the optimization method when the Session file is configured (see Supplementary Method S1 for how to create a session file) The implementation of the application is based on the model-view-controller or MVC architecture. The ''model'' contains the analysis or configuration. This can be executed from the MATLAB Ò command-line independently from the front-end viewer. The model is represented by the back-end objects with the ''?PKPD'' MATLAB Ò package folder. The ''viewer'' designates the graphical user interface (GUI), which is mostly contained within the ''?PKPDViewer'' package. Implementation of the viewer is based on the ''GUI Layout Toolbox'', a programmatic layout manager, from MATLAB Ò Central TM . The ''controller'', represented by the viewer's callback functions, responds to inputs by the user by updating the model and subsequently the viewer. The controller is also contained with the ''?PKPDViewer'' package folder.
The application has been tested in MATLAB Ò R2016a and R2016b. The software is not compatible with earlier MATLAB Ò releases and currently has not been tested for releases after MATLAB Ò R2016b. The application requires MATLAB Ò and is not supported for distribution via MATLAB Ò Compiler. Both Windows and Mac operating systems are supported. The software is available as Supplementary Material and on MATLAB Ò Central TM (The gPKPDSim and NCA toolboxes are available at: https://www.mathworks.com/matlabcentral/fileexchange/ 65399-gpkpdsimtoolbox and https://www.mathworks.com/ matlabcentral/fileexchange/65303-simbiologynca)

Library of models
With this manuscript, we have provided a set of PKPD models that are frequently used in the pharmaceutical industry; however, any other model built in SimBiology Ò can be easily configured for gPKPDSim. The models in our library include (1) two-compartment pharmacokinetics model with IV and extravascular dosing. The model includes both non-specific clearance (typically captured by a linear term) and specific clearance (typically captured by a nonlinear Michaelis-Menten term to model target-mediated drug disposition). (2) Target-mediated drug disposition (TMDD) model [5], (3) physiologic indirect response models (with inhibitory or stimulatory effects on synthesis or degradation of some endogenous target or mediator) [6] and (4) minimal physiologically based pharmacokinetic (PBPK) model with target represented in the central compartment, leaky and tight tissues [7]. Each SimBiology Ò model is encapsulated in a ''Session'' file, that can be loaded in gPKPDSim by end-users for performing PKPD tasks. The schematic diagram of each model in the library is shown in Fig. 1.

GUI application overview
In the framework developed in this work, model development performed by the expert modeler is separated from utilization of the model, a task that can be performed by a broader range of scientists and researchers. Thus, using the framework shown schematically in Fig. 2, standard PKPD models can be developed once by the expert and used repeatedly by end-users in diverse projects.
The expert modeler builds the model in SimBiology Ò and configures it into a ''Session'' file, which contains specification of the SimBiology Ò model Variants, Doses, Species, Parameters and Settings for use within the user interface (see Supplementary Method S1 for how to create a session file). The Session file, once loaded in gPKPDSim, enables functionalities that include (1) simulation, (2) data fitting (parameter estimation), (3) population simulation, and (4) NCA. All these functionalities are supported by some common features such as interactive visualization, export of results as presentation-ready figures and Excel datasets, and saving the working Session files for future use. The end-user works exclusively in the environment of gPKPDSim and can send the saved session files to communicate with collaborators and share the results. Table 1 shows the description of different components of the gPKPDSim framework, and Fig. 3 shows a snapshot of gPKPDSim GUI.

gPKPDSim general settings
This section includes the common settings between the different functionalities. Table 2 lists a summary of the general settings and menus in the GUI. A more detailed description of variants and dosing are provided below (see Supplementary Method S3 for more details on the remaining gPKPDSim settings).

Variant
Displays a list of built-in variants from the SimBiology Ò model. A variant stores the names and values of parameters of the SimBiology model and allows the user to use values stored in a variant object as the alternate value to be applied during a simulation. For example, if you have different values for PK parameters in animals and humans, those can be differentiated through variants (see Fig. S1 for a visual example and https://www.mathworks.com/help/ simbio/ref/variantobject.html for more information). Each session file can select from the variants available in the SimBiology Ò sbproj file; not all variants from the sbproj file need to be included in the session file. If a single variant is activated, the values of the parameters in gPKPDSim will be overwritten by the values stored in the variant. If multiple variants are activated, the variants will overwrite the parameter values in an ascending order. The user then can change parameter values in the functionalityspecific section after activating the variants. At any point, the user can reset the value to default by clicking ''Reset to defaults'' button. To summarize, here is the order in which the parameter values are assigned:

Simulation
Enables simulation of the model encapsulated in the session file with different parameter values, variants, dosing regimens and simulation time. Once the simulation is completed, the user can visualize results using the interactive plotting features. The user can choose to run multiple simulations and retain the results of previous simulations for the purpose of visualization or results export to Excel files. The profile notes section of gPKPDSim includes a detailed summary of each simulation run. In addition to simulation, the user can also import data for visualization and comparison with the simulation results (Supplementary Method S3).

Non-compartmental analysis (NCA)
Estimate the PK parameters for the imported dataset. The dataset can include any combination of dosing route (IV and/or extravascular dosing), and schedule (single-dose and/or multi-dose scenarios). gPKPDSim can compute NCA for sparse (a PK sampling scheme in which not all subjects were sampled at all the time points) and serial (a PK sampling scheme in which all subjects were sampled at all the time points) schemes with automatic detection of the dosing values from the dataset (see Supplementary Method S2 for the dataset formats). It also allows the user to provide an array of time ranges to estimate C max and partial AUC values. For calculating the terminal half-life, the • Group: Shows the identifier for the Group column header in the dataset. As mentioned before, a Group is defined as a number of Subjects receiving the same dose material, through the same dose route, and according to the same dose schedule. • Dose schedule: Frequency of administration. This could be single or multiple. • Administration route: Intravenous or extravascular administration. • Lambda_Z: First-order rate constant of elimination (commonly referred to as the terminal slope). • R2: A coefficient of determination that is a statistical measure of how well the regression line for half-life determination approximates the real data points.  The expert modeler develops the model in SimBiology Ò and encapsulates it into a ''Session'' file, which can be launched in gPKPDSim by end-users. The Session file, once loaded in gPKPDSim, enables the end-user to perform different functionalities, including simulation, data fitting (parameter estimation), population simulation, and NCA. For some of these functionalities, the user must provide a dataset. All the functionalities are supported by features such as interactive visualization, and export of results as presentationready figures and Excel datasets. When the user completed their tasks, the session file can be saved for future use. Note that the end-user works exclusively in the environment of gPKPDSim • AUC_extrap_percent: Percentage of the AUC_infinity value that is from extrapolation.   This functionality allows the end-user to estimate model parameter values, providing a best fit of model simulation to data. The user imports the dataset in gPKPDSim and maps all the data outputs to be used in data fitting as well as the dosing information for each group of subjects. The user can also select the error model for the objective function, specify if pooled data should be used to perform the data fitting and select parameters to be estimated. After the data fitting algorithm is complete, the user can save estimated parameter values as a new variant (which will appear in the list of variants and the tag column will show ''Custom'') within the current session, enabling use of the parameter estimates for subsequent simulation or population simulation functionality (Supplementary Method S3). Note that the data fitting functionality in gPKPDSim does not provide parameter estimation for mixed-effect models, which is available in other software such as SimBiology Ò , NONMEM Ò and Monolix.

Population simulation
The purpose of this functionality is to explore the impact of parameter variability on the model outputs of interest. The user specifies the distribution of parameters to be used; this can be based on expert knowledge or based on the standard errors in parameter estimates from the fitting task. Note that the population simulation functionality is not based on parameter distributions estimated from mixed-effect models, but relies on the end-user to provide the variability of each parameter. And example of how covariates are handled in gPKPDSim is provided in Supplementary Method S3. The user also specifies the number of simulations, i.e., the number of times the parameter distributions are sampled and the model is run. Once the population simulation is complete, the visualization plots display the median curve as well as a shaded region corresponding to the 5-95% boundaries for the output selected. Similar to other functionalities, the user can save the plots and/or export the results to an Excel file.

Results
Case Study #1: NCA, two-compartment model fit, and multi-dose projection for PK of a large molecule A common scenario encountered in pharmaceutical R&D is the use of data obtained with one dose regimen to characterize PK and to then project PK for other dosing scenarios. In this example, gPKPDSim is used to perform NCA [8] for PK data for single dose administration of a monoclonal antibody [9], followed by fitting of the data to a two-compartment PK model, which is then used to simulate and predict multi-dose PK and population variability. Thus, this example uses NCA, fitting, simulation, and population simulation functionalities of gPKPDSim. The instructions in Supplementary Instruction S1 can be followed for working through this exercise (see Table S5 for initial conditions). The dataset used in the analyses is provided in the supplement and includes systemic PK measurements for the study design in Table 3. NCA outputs are calculated for two scenarios: (1) sparse (or group-level) analysis resulting in two sets of NCA parameters, one for each dosing group (Table 4) and (2) serial (or sub-group-level) analysis resulting in six sets of NCA parameters, one for each animal and ( Table 5). The NCA results are computed to NCA results estimated by WinNonlin Ò for both cases (Table 4 and Supplementary  Table S3). We also included thirteen additional datasets for different dosing regimens including IV and extravascular, single and multi-dose scenarios (see Table S4 and Supplementary Files). In all cases, the absolute difference between the NCA results of gPKPDSim and those of WinNonlin Ò was less than 0.1%, when the number of data points used for calculation of parameters are matched. As seen in Table 5, the estimated value of clearance (CL) varies between 6.01 and 8.02 mL/kg/day across animals and the range for estimated volume of distribution (V_z) is 74.49-81.84 mL/kg. Next, we use the data fitting functionality to fit the PK data to a two-compartment model to estimate the PK parameters. The standard two-compartment antibody PK model with non-specific (captured by a linear term) and specific (captured by a nonlinear term) clearance used here is described by the following set of differential equations: where A C , A P , and A SC represent the amount of drug in the central, peripheral, and extravascular (subcutaneous) compartments; X C and X P denote the drug concentration in the central and peripheral compartments; and AUC represents the area under the curve for X C . The units for antibody amounts and concentrations are lg/kg and lg/mL. Parameters V 1 and V 2 represent the volume of central and peripheral compartments, CL represents the clearance from the central compartment, CL d denotes the distribution clearance between the central and peripheral compartments, k abs represents the absorption rate from the extravascular compartment and f bio is the bioavailability, that is the fraction of drug available in the central compartment after extravascular dosing (see Supplementary File casestudy1_TwoCompPK_equations.pdf for the units of parameters and species in the model).
The quality of the parameter estimation can be evaluated by the value of standard error on the fitted parameter values and the visual predictive checks, included in the Fitting Summary report (see Supplementary File casestudy1_fitting_summary_combined_pooled.pdf). For this dataset, the non-linear clearance term is not required to explain the data, as the drug exposure is dose-proportional for the low and high dose groups. The results for the pooled fitting are shown in Fig. 4a. The parameter estimates for the pooled and group-specific fittings, shown in Table 6, have a relatively low standard error. The estimated value of CL is 6.89 mL/Kg/day, which is consistent with the NCA results. The central and peripheral compartment volumes are estimated at 49.1 and 34.6 mL/Kg, respectively. Note that the performance of data fitting in gPKPDSim is essentially that of ''Fit Data'' task in SimBiology Ò .
Finally, we use the model and parameter estimates together to project PK for an alternate dosing regimen and explore potential PK variability. The parameter estimates from the pooled-fitting were saved in the form of a ''variant'', and then used for projecting the PK profile of a 10 mg/kg 4qw (weekly dosing for 4 weeks) IV dosing  regimen under different conditions listed in Table 7 to explore the effect of the PK parameters on drug exposures; this was done using the simulation and population simulation functionalities. The PK profiles are shown in Fig. 4b and d and the AUC curves are depicted in Fig. 4c and e (AUC d0-28 for the three scenarios are reported in Table 7). The simulation results suggest that decreasing the clearance by half results in a 38% increase in AUC for this molecule. Additionally, assuming 10% population variability on V 1 , and CL results in a median AUC value of 4065 lg 9 day/mL with the 5-95% boundaries being 25% lower or higher than the median.
Case Study #2: target-mediated drug disposition model to predict antibody PK and target engagement In addition to characterizing, fitting, and projecting PK as in example 1, antibody drug development often involves projecting target occupancy as a function of dose regimen and antibody affinity, in order to determine the antibody affinity needed or identify the dosing requirements for desired target engagement. In this case study, we consider a two-compartment model with target mediated antibody drug disposition in the central compartment [5]. This model is an extended version of the two-compartment model, in which the non-linear Michaelis-Menten clearance term is replaced with target mediated drug disposition. The model equations are shown below: dC CðnMÞ dt ¼ k on Á X C;free nM ð Þ Á T C;free nM ð Þ À k on Á K D Á C C nM where A C, free , A P, free , and A SC, free represent the amount of free drug in the central, peripheral, and SC compartments; the concentration of free target and the drug:target complex is represented by T C, free(nM) and C c(nM) ; the concentration of free drug in the central and peripheral compartments is denoted by X C, free and X P, free . Unless otherwise noted, the units for antibody amounts, antibody concentrations, and target concentrations are lg/Kg, lg/mL, and ng/mL. Compared to the model presented in Case Study #1, this model has additional parameters including k on , K D , MW Ab , MW T , t 1/2,T , T c,init(nM) , and a CL representing the drug-totarget association rate constant, equilibrium dissociation constant, molecular weight of drug, molecular weight of target, half-life of target, initial concentration of target, and the ratio of the clearance of the bound drug-target complex to free drug clearance. In this model the target is considered to be a soluble target present in circulation. The stoichiometry of binding is taken to be 1:1 for antibody and target to form the complex (see Supplementary File casestudy2_TMDD_equations.pdf for the units of parameters and species in the model).
The instructions in Supplementary Instruction S2 can be followed for working through this study (see Table S5 for initial conditions). The goal is to evaluate the impact of (1) antibody-target binding affinity and (2) dose level on target engagement. The settings for the simulations are described in Table 8. Figure 5a-c show the results for the TMDD model simulations. Comparing Simulations 2.1 and 2.2 (K D = 0.1 and 10 nM), we see that the fraction of the target that is antibody-bound remains high in the simulation with the higher affinity and drops with reducing antibody concentration but it is low over the course of the simulated time course using the low affinity antibody (Fig. 5c ). The exposure of the antibody in both simulations is comparable (Fig. 5a), indicating that while the difference in affinity significantly influences target binding, it does not impact drug PK. The profiles for free and total (bound ? free) target for the two simulations are shown in Fig. 5b. The effect of dose is seen where we compare the target engagement of Simulation 1 (1000 lg/kg dose Q4 W) with Simulation 2.3 (5000 lg/kg dose Q4 W). In this case, the model suggests that with 5 times higher dose and systemic drug exposure, the target remains [ 98% bound over the entire simulation time-window.
The case study demonstrates how scientists can use gPKPDSim to answer some typical questions in drug development such as ''knowing the target levels and turn over rate, at what dose level does the impact of TMDD on the PK profile become negligible?'', ''Is the profile of total antibody different from the free antibody?'', etc. The enduser can also test ''what if'' scenarios for the particular antibody of interest, for example by selecting dose levels, frequency and changing parameter values using sliders. The interactive plots allow visualization of different outputs across multiple simulations with appropriate scales  Table 7 Estimated AUC values for the period of 0-28 days of a 10 mg/kg 4qw IV dosing regimen under different conditions ( Fig. 4c and e) gPKPDSim functionality Scenario AUC d0-28 (lg day/mL) Simulation CL 1 = 6.89 (result of fitting) 4129 Simulation CL 2 = 0.5 9 CL 1 = 3.45 5690 Population simulation CL3 = 6.89 ? 10% CV Case Study #3: physiologic indirect response models to link PK with PD biomarkers Beyond PK and target engagement or inhibition, another important consideration in pharmaceutical sciences is the biological impact as measured by pharmacodynamics (PD) biomarkers. Basic PD models have been classified as direct or indirect models. This section illustrates the use of gPKPDSim with a two-compartment PK model integrated with either one of four basic physiologic indirect response models to characterize pharmacodynamics effects of a drug, whose mechanisms of action can be summarized as inhibition or stimulation of the synthesis or degradation of some endogenous target or mediator [6]. The selection of the desired indirect response model is provided as a variant selection in gPKPDSim. The model is described by the following set of differential equations: ; ð26Þ Response ¼ R 1 for Inhibition on Synthesis R 2 for Inhibition on Degradation R 3 for Stimulation on Synthesis R 4 for Inhibition on Degradation where Eqs. (19)-(22) describe the two-compartment PK model presented in Case Study #1 and Eqs. (23)-(26) describe the four basic physiologic indirect response models and the impact of drug concentration on synthesis (k in ) or degradation (k out ) of factors controlling the response variable (R x ) via an inhibitory (IC 50 ) or stimulatory (EC 50 and E max ) mechanism. Parameters IC 50 and EC 50 have the units of lg/mL, whereas k in , and k out take the units that suit the problem of interest. Four variants are included in the session file, corresponding to the four basic models. Note that only one of these variants should be active at any given time (see Supplementary File cases-tudy3_IDR_TwoCompPK_equations.pdf for the units of parameters and species in the model).
The first example (Simulation 3.1) in this section describes the effects of pyridostigmine on inhibition of cholinesterase, which leads to increased levels of acetylcholine in the neuromuscular junction and subsequently increases muscular activity, which is considered the PD effect in the treatment of myasthenia gravis [10,11]. Here, we aim to replicate the PK profile and the muscle response ( Fig. 6 from [6]) in a patient that received a 5 mg dose of intravenous pyridostigmine. The PK of pyridostigmine has a bi-exponential decay and can be described by a twocompartment model with parameter values listed in Table 9 (the calculations are done for a 70 kg patient). For the PD effect, we use the second indirect response model (R 2 ) with the inhibition on degradation. Profiles generated in Fig. 6a and b characterize the PK profile of pyridostigmine and its PD effect and corroborate the results reported in [6].
The second example (Simulation 3.2) describes the impact of cimetidine on stimulating secretion of prolactin [12]. The plasma concentration of cimetidine after IV administration of 300 mg in patients can be characterized by a two-compartment model. For the PD effect, we use the third indirect response model (R 3 ) with the stimulation on synthesis. Profiles generated in Fig. 6c and d describe cimetidine PK profile and its effect on prolactin concentration and agree with the results reported in Fig. 9 of [6]. The case study presented here showcases how scientists can reproduce and validate models from literature, use them to inform their project and help them design their studies. The instructions for this case study is provided Supplementary Instruction S3 (see Table S5 for initial conditions).
Case Study #4: minimal physiologically based pharmacokinetic (PBPK) model In this example, we use a minimal antibody PBPK model, which incorporates target mediated drug disposition by including target-binding in plasma and two groups of tissues, namely, leaky and tight tissues [7]. As in the original publication, the model assumes that only free antibody flows between compartments and the drug-receptor complex is immobile. The binding equations in the model were revised to fix a typographical error in the original paper [personal communication with the author]. Additionally, the model now enables representation of the target in all compartments, simultaneously. The differential equations for this model are: ð30Þ dA TT;free dt ¼ 1 À r TT ð ÞÁL TT Á X C;free À 1 À r L ð ÞÁL TT Á X TT;free À k on Á X TT;free nM ð Þ Á T TT;free nM dC CðnMÞ dt ¼ k on Á X C;free nM ð Þ Á T C;free nM ð Þ À k on Á K D Á C C nM Similar to Case Study #2, A C, free , A TT, free , A LT, free , and A L, free represent the amount of free drug in the central (plasma), tight tissues, leaky tissues and lymph. The concentration of free target and the drug:target complex in plasma is represented by T C, free(nM) and C c(nM) . The concentration of free drug in plasma is denoted by X C, free . Similar nomenclature was used to define the concentration of drug, target and complex in leaky and tight tissues. Unless otherwise noted, the units for antibody amounts, and concentrations are lg, and lg/L, whereas the units for the concentration of target and complex are nM. The volume of plasma, tight and leaky tissues and lymph are represented by V C , V TT , V LT , and V L . Parameters r TT and r LT represent vascular reflection coefficients for tight and leaky tissues and L TT and L LT are lymph flows in these tissues. Parameters L and r L denote lymphatic flow and reflection coefficient. Parameters k on , K D , MW Ab , and k int represent association rate constant, dissociation constant, molecular weight of drug, and internalization rate of the drug:target complex, whereas k x,deg , and T x,init(nM) denote the degradation rate constant and initial amount of target in compartment x (see Supplementary File cases-tudy4_minPBPK_final.pdf for the units of parameters and species in the model).
Here, we replicate the simulated plasma concentrations in the presence of 10 nM target in plasma, leaky or tight tissues (Fig. 2 from [7]) and compare the results with the case of no target. The parameter values are taken from Table 1 from [7], where CLp should have been listed as 10 -4 L/h/kg. Figure 7a and b show nonlinear profiles when target is present in circulation or tight tissues where nonlinearity is more pronounced for the target in the central compartment. In contrast, in tight tissues, the PK profile appears to be dose proportional (Fig. 7c). When compared to the case with no target, we observed that overall clearance is slightly higher when target is present in tight tissues (Fig. 7d). Note the difference in PK profiles depicted in this paper and Fig. 2 of [7] (this is due to the typographical error in the original binding equations, corrected here). The instructions in Supplementary Box 4 can be followed for generating Fig. 7a-d (the instructions for this case study is provided Supplementary Instruction S4; see Table S5 for initial conditions). Once again, the case study presented here showcases how researchers can use gPKPDSim to corroborate the results from published papers, compare the results with other scenarios that may not have been included in the original publications and use those to informal internal decisions.

Discussion
With the increasing application of PKPD modeling in pharmaceutical research and development, one key challenge is to efficiently deploy models built by expert modelers and enable scientists and collaborators with limited M&S experience to correctly use these models to answer research questions. To do so, we have developed gPKPDSim, a multi-purpose MATLAB Ò application with a user-friendly interface for deploying SimBiology Ò models and allowing scientists to perform various PKPD analyses and model-based tasks, regardless of their modeling expertise. The application has a single screen interface to perform common PKPD tasks and includes additional capabilities that may not currently exist in SimBiology Ò , but ''hides'' the MATLAB Ò code and involved features of SimBiology Ò from end-users. And therefore it removes a significant bottleneck of software training for most end-users. Furthermore, modeling experts need not perform all routine analyses, but can serve instead as advisors and reviewers on these tasks, increasing their bandwidth for more complex modeling and analysis efforts.
In this work, we have used MATLAB Ò and SimBiology Ò as the engine behind gPKPDSim because they provide a powerful and flexible platform to build, simulate, and analyze mechanistic PKPD and systems biology models. Note that several programming and modeling platforms are available to the M&S community and similar applications can be developed in some of these platforms as well. Note that gPKPDSim provides PKPD M&S capabilities to users in an organization with access to SimBiology Ò /MATLAB Ò at no additional cost; however, each organization has to make an informed investment based on their needs in software tools such as SimBiology Ò /MATLAB Ò and others that have costs associated with them.
Along with the gPKPDSim application, we have provided a library of PKPD models including (1) a twocompartment pharmacokinetics model with specific and non-specific clearance, (2) a target-mediated drug disposition (TMDD) model, (3) a physiologic indirect response model and (4) a minimal physiologically based pharmacokinetic (PBPK) model with the drug target represented in the central compartment, leaky and tight tissues. These models have a wide range of applications in drug development such as characterization of PKPD relationships, molecule selection, human PK projection, and preclinical and clinical dose and schedule determination, especially for antibody therapeutics, and the application can be used with additional models (e.g., small molecule PK). The case studies presented in this paper demonstrate how scientists familiar with the concepts of PKPD modeling can efficiently use gPKPDSim to analyze PKPD data and address their drug research and development questions. Furthermore, the application can be used for exploratory research with PKPD models for hypothesis testing and evaluation of what-if scenarios. It also enables easy sharing and communication of the results.
The application is open-source; thus the M&S community can extend software as they see fit for their respective use. Additionally, gPKPDSim works with any SimBiology Ò -based model, providing flexibility to researchers developing their own models. gPKPDSim is also equipped with the NCA functionality, which can handle datasets that have IV and/or extravascular dosing regimens in singledose or multi-dose scenarios with sparse and serial sampling schemes. This is an extra benefit over WinNonlin Ò , in which these different scenarios must be provided in multiple datasets to calculate the NCA results. The population simulation functionality in gPKPDSim is intended for evaluation of parameter variability on PKPD predictions primarily in the preclinical and translational settings; this functionality should not be used for applications that require a formal population PK model. At Genentech, gPKPDSim was initially developed for and is now regularly used by PKPD scientists working on preclinical/translational PKPD of biologics and large molecules. We have found that by empowering pharmaceutical scientists and drug development team members to conduct straightforward model-based tasks, gPKPDSim encourages broader use of PKPD models by collaborators across multiple functions to support drug development, optimize study design, reduce animal usage and help with decision-makings.