ATLAS mPBPK: a MATLAb-based tool for modeling and Simulation of minimal Physiologically Based PharmacoKinetic models.

Minimal physiologically-based pharmacokinetic (mPBPK) models are frequently used to model plasma PK data and utilize and yield physiologically-relevant parameters. Compared to classical compartment and whole-body PBPK modeling approaches, mPBPK models maintain a structure of intermediate physiological complexity that can be adequately informed by plasma PK data. In this tutorial, we present a MATLAb-based tool for modeling and Simulation of mPBPK models (ATLAS mPBPK) of small and large molecules. This tool enables the users to perform: i) PK data visualization, ii) simulation, iii) parameter optimization, and iv) local sensitivity analysis (SA) of mPBPK models in a simple and efficient manner. Along with the theoretical background and implementation of the different tool functionalities, this tutorial includes simulation and SA showcases of small and large molecules with and without target-mediated drug disposition. This article is protected by copyright. All rights reserved. This article is protected by copyright. All rights reserved.

The evaluation of pharmacokinetics (PK) enables the assessment of drug exposure and access to target sites along with an overall understanding of the biological processes determining absorption, distribution, metabolism, and excretion from the body. Based on the question of interest and the amount and quality of the experimental data, PK can be evaluated by noncompartmental, compartmental, and physiologically‐based models.1 Frequently, noncompartmental and simple compartment models are characterized as “top‐down” approaches because their implementation and performance do not require prior knowledge of the system under investigation and are solely based on the observed PK data. On the other hand, physiological models are “bottom‐up” approaches and require a deeper appreciation of the physiology of the species of interest and the mechanisms involved in drug absorption, distribution, metabolism, and excretion. The use of both methods is based on the available PK data and intended purposes. In general, physiological models require numerous equations and a significant amount of data that many times are not easily accessible (e.g., tissue PK), whereas “top‐down” models are easy to implement, but their use is limited to describe the PK of the available data and for limited extrapolations.2

Recently, minimal physiologically‐based pharmacokinetic (mPBPK) models were introduced and shown to provide more useful assessments of PK properties than classical compartment models.3, 4, 5, 6, 7 Compared with whole‐body physiologically‐based pharmacokinetic (PBPK) models, mPBPK models limit physiological information to lumped tissues, and sought parameters can be estimated by the sole use of plasma PK data. The mPBPK models provide a realistic basis for describing plasma PK data and differ from compartment models in ways of initial distribution space, physiological assignments, and restrictions as well as flexibilities in handling different clearance mechanisms. Importantly, because of this general physiological relevance, mPBPK models are attractive frameworks to evaluate drug–drug interactions and provide interspecies scaling.8, 9

Because of the advancements of computer science and in silico methods of calculating physiologically relevant parameters (e.g., partition coefficient), PBPK models have become very popular during the past decade.10 Given this growth of interest, there is increasing availability of commercial platforms that integrate physiologically based methodologies. Examples of such platforms are Simcyp (https://www.certara.com/software/pbpk-modeling-and-simulation/), GastroPlus (https://www.simulations-plus.com/software/gastroplus/), SimBiology (https://www.mathworks.com/products/simbiology.html), and PK‐Sim (http://www.open-systems-pharmacology.org/). Discussion articles on best practices of PBPK model development and their use to address regulatory questions further underline their importance.11, 12 Currently, mPBPK models can be developed through typical software platforms such as WinNonlin (Certara, Princeton, NJ, USA), NONMEM (ICON plc, Gaithersburg, MD, USA), ADAPT (BMSR Biomedical Simulations Resource, University of Southern California, Los Angeles, CA, USA), MATLAB (MathWorks, Natick, MA, USA), and many others as there is no specific computational tool needed for mPBPK modeling and simulation except for numerical integration. However, the critical steps of model development, verification, and validation such as sensitivity analysis (SA) of PK parameters are not part of the software design, and most of the time are cumbersome to perform.13, 14, 15, 16 In this tutorial, we describe a new MATLAB‐based tool for the modeling and simulation of mPBPK models (ATLAS mPBPK). This software offers a number of different predefined mPBPK models based on which the user can perform the following: (i) simulations, (ii) parameter optimization, and (iii) SA. ATLAS mPBPK incorporates models for small and large molecules, with the possibility of incorporating target‐mediated drug disposition (TMDD) in the plasma or interstitial fluid (ISF) of peripheral tissues. Compared with already established PBPK software, ATLAS mPBPK is a tool solely focused on mPBPK modeling that offers an easy and straightforward parameter estimation and SA framework where the user can perform the respective functions by solely using a number of checkboxes/editboxes. Furthermore, it is an open‐source, freely available tool that can be downloaded from SourceForge (https://sourceforge.net/projects/atlas-mpbpk/files/ATLAS_mPBPK_v.1/).17 Details on installation can be found on the README.md file on the ATLAB mPBPK repository.

ATLAS mPBPK User Interface

The ATLAS mPBPK interface is simple and includes a number of checkboxes, drop‐down menus, and push buttons through which the user can perform simulations, parameter estimation, and SA. The interface is shown in Figure  1, and it is divided into five groups of elements. In the Figure  1 a group, the user is able to load the data set of interest, choose the x/y axes, and plot the data (Figure  1 e). By clicking the Load.xlsx data file button, a window opens, and the user can browse computer files to locate the .xlsx file where the PK data of interest are stored. The format of the .xlsx file should include a title row. After loading the data, the user can choose the x and y axes from the drop‐down menus of Figure  1 a. Upon selection, the data are automatically plotted in Figure  1 e of the interface as a typical PK graph. The user can change the x/y axes at will by choosing different values of the Figure  1 a x/y‐axes drop‐down menus as well as changing the axes of the plot from linear–linear to log–linear by clicking the Linear/Log radio button at the upper right corner of the plot.

Figure 1

ATLAS mPBPK user interface. (a) Data visualization elements. (b) Minimal physiologically‐based pharmacokinetic (mPBPK) model selection, simulation time/administration protocol, and target‐binding configuration. (c) mPBPK model parameter set‐up, checkboxes for parameter estimation, configuration of LB, UB for parameter estimation, and checkboxes for sensitivity analysis. (d) Buttons for different ATLAS mPBPK functions execution. (e) Plot for the different ATLAS mPBPK functions illustration. C.I., confidence interval; LB, lower bound; PO, oral; UB, upper bound.

In the Figure  1 b group of elements, the user can choose the mPBPK model of interest and the administration protocol as well as details for target binding and the simulation time. In particular, in the mPBPK model drop‐down menu, the user can choose among two large groups of mPBPK models, small molecule or large molecule, depending on the compound of interest. Figure  2 shows the frameworks of the different mPBPK models that are used in ATLAS mPBPK. Figure  2 a is the mPBPK model for small molecules,7 and Figure  2 b–d is the mPBPK models for large molecules.3, 4 If the user selects Large molecule from the mPBPK model drop‐down menu, an additional Target binding element is shown (Figure  1 b). Through three radio buttons, the user can choose if the compound of interest maintains TMDD in plasma (Central, Figure  2 c), in ISF of peripheral tissues (Peripheral, Figure  2 d), or does not bind to a target (No binding, Figure  2 b). The time of the simulation can be edited in the Simulation time element (Figure  1 b). Finally, in the Administration element, the user can choose among the following four administration protocols: Intravenous (IV), Extravascular (EV), Infusion (INF), and Oral (PO). If user chooses Infusion (INF), a new element appears where the user can set the infusion time (Figure  1 b). Currently, ATLAS mPBPK can handle only single‐dosing protocols. A simultaneous analysis of dose ranges is not possible in this version.

Figure 2

The minimal physiologically‐based pharmacokinetic (mPBPK) models that are used in ATLAS mPBPK. (a) mPBPK model for small molecules. (b) mPBPK model for large molecules with no target binding. (c) mPBPK model for large molecules with target binding in plasma (interstitial fluid). (d) mPBPK model for large molecules with target binding in peripheral tissues. Parameters' names and typical values are shown in Table 1. EV, extravascular; IV, intravenous; PO, oral.

Table 1

Typical parameter values used in ATLAS mPBPK

Parameters	Typical value (human)	Units	Reference
Small‐molecule model
Cardiac output—Q co	5.6	L/minute	7
Portal vein blood flow—Q hep	1.45	L/minute	32
Body weight or extracellular fluid	70	kg	Typical human body weight
Blood volume—Vp	5.2	L/70 kg	32
Liver volume—V hep	1.69	L	7
Highly perfused tissue volume—V1	24.3	L	7
Cardiac output fraction to highly perfused tissues fd1	0.7	—	Drug related or estimated from data
Cardiac output fraction to lower perfused tissues fd2	0.1	—	Drug related or estimated from data
Partition coefficient—K p	0.7	—	Drug related or estimated from data
Hepatic intrinsic clearance—CLintu	0.7	L/minute	Drug related or estimated from data
Nonhepatic clearance—CLnh	0.01	L/minute	Drug related or estimated from data
Large‐molecule model
Lymph flow—L	2.9	L/day	33
Plasma volume—Vp	2.6	L	3
Interstitial fluid—ISF	15.6	L	3
Partition coefficient—Kp	0.8	—	Drug related or estimated from data
Lymph volume—VL	5.2	L	21
Lymph refl. coefficient—sigmaL	0.2	—	Drug related or estimated from data
Vascular reflection coefficient for tight tissues—sigma1	0.95	—	Drug related or estimated from data
Vascular reflection coefficient for leaky tissues—sigma2	0.512	—	Drug related or estimated from data
Plasma clearance—CLp	0.0054–0.03	L/hour	3
Steady state constant—Kss	0.1	nM	Drug related or estimated from data
Target biosynthesis rate—k syn	0.001	nM/hour	Drug related or estimated from data
Free target degradation rate—kdeg	0.1	1/hour	Drug related or estimated from data
Complex internalization rate—k int	0.0117	1/hour	Drug related or estimated from data
For PO/EV
Bioavailability	1	—	Drug related or estimated from data
Absorption rate–ka	0.5	1/hour	Drug related or estimated from data

EV, extravascular; PO, oral.

In the Figure  1 c group of elements, the user can edit the value of the mPBPK model parameters, choose what parameters need to be estimated by also setting their lower and upper bounds, and test certain parameter sensitivities. Based on the mPBPK model that the user chooses (Figure  1 b), the mPBPK parameters shown in Figure  1 c are different and in agreement with the respective mPBPK model (e.g., small‐molecule models vs. large‐molecule models). As initial values, ATLAS mPBPK includes Values typical for human physiology. Typical parameter values along with their bibliographic reference are shown in Table  1. In the Estimate checkboxes, the user can check what parameters to optimize. In the LB, UB editboxes, the user inserts the lower and upper bounds of the parameters with the format: lower bound,upper bound (e.g., 0.1,1). To perform parameter optimization, the lower bound, upper bound (LB, UB) need to be inserted in the tool. After successful optimization, optimal parameters are outputted in the Estimates editboxes with their 95% confidence interval on 95% CI editboxes. Finally, in the Sensitivity checkboxes, the user can select the parameters to be tested for their sensitivity.

The Figure  1 d push‐buttons initiate the respective functions of ATLAS mPBPK. These functions are simulation (Simulate), parameter estimation (Estimate), SA (Sensitivity), and export of the simulation results to an excel file (Export results). For simulations (Simulate), the user must provide the Values of the mPBPK model parameters. For parameter estimation (Estimate), the user must choose the parameters of interest by checking the appropriate checkboxes (Estimate) and provide the lower and upper bounds of the optimization in the LB, UB editboxes. For SA (Sensitivity), the user chooses the parameters of interest in the checkboxes (Sensitivity). Finally, by clicking the Export results button, an excel file named as mPBPK_Results with all of the model numbers is generated in the directory where ATLAS mPBPK is saved.

The user is responsible for the consistency of mPBPK model units. Next to each parameter, square brackets indicate the general group of units required (e.g., mass, volume, time, etc.). Furthermore, in the current ATLAS mPBPK version, there is no specific consideration of active processes (e.g., enzymes, transporters) and formulation parameters. Their effects are indirectly considered through parameters such as absorption rate, nonhepatic clearance, and so on. Lastly, in ATLAS mPBPK, there is no explicit differentiation between subcutaneous and intramuscular administration.

ATLAS mPBPK Simulations

The MATLAB's ode23s function is used to integrate differential equations of the mPBPK models. Supplementary Figure shows part of the ATLAS mPBPK code where the right‐hand side of the differential equations of small (Supplementary Figure a) and large molecules (Supplementary Figure b) are defined. As shown in Supplementary Figure , the user‐supplied parameters (Figure  1 c, Values) are treated as structure arrays (e.g., s.fd1, s.fd2, s.Qhep, etc.). For the mPBPK model of large molecules (Supplementary Figure b), three additional structure values are introduced (s.central, s.peripheral, s.nb) that take values of 0 or 1 depending on the characteristics of TMDD. If there is no TMDD, then s.nb = 1/s.central = 0/s.peripheral = 0; if TMDD is in plasma, then s.nb = 0/s.central = 1/s.peripheral = 0; if TMDD is in the ISF of tissues, then s.nb = 0/s.central = 0/s.peripheral = 1.

Small molecules

The mPBPK model for small molecules as previously published7 incorporates two tissue compartments as well as a liver compartment to describe oral dosing with hepatic first pass (Figure  2 a). The differential equations of the small‐molecule mPBPK model are: where C p is the concentration of the drug in V p (blood or plasma volume), C 1 and C 2 are drug concentrations in tissue compartments 1 (V 1) and 2 (V 2), Q CO is cardiac blood (or plasma) flow, f d1 and f d2 are fractions of Q CO for V 1 and V 2, K p is the partition coefficient for tissues 1 and 2, CLintu and CLnon‐hep are unbound hepatic intrinsic and nonhepatic clearances, Q hep is the portal vein blood flow, C hep is the drug concentration in the liver (V hep), and F G and k a are the prehepatic bioavailability and the absorption rate constant. The constraints of this model are: where BW is the body weight and ECF is the extracellular fluid volume.

Figure  3 shows a representative ATLAS mPBPK simulation for 600 mg intravenous dosing of benzylpenicillin. The parameters of the mPBPK model were taken from ref. 7 and the experimental data from ref. 18. In the original simulations, the physiology‐related parameters (e.g., Q CO, Q HEP, BW, V p) were kept constant to the literature values, whereas drug‐specific parameters (e.g., K p, CLintu/nh, f d) were optimized based on the PK data. Figure  4 further depicts ATLAS mPBPK simulations with respect to the PK data of beta‐lactams taken from ref. 7. Similarly, drug‐related parameters were the only parameters that were estimated.

Figure 3

Simulation example for 600 mg intravenous dosing of benzylpenicillin. The minimal physiologically‐based pharmacokinetic (mPBPK) parameter values were taken from ref. 7 and the PK data from ref. 18. C.I., confidence interval; LB, lower bound; UB, upper bound.

Figure 4

Concordance of simulated ATLAS mPBPK small‐molecule model values based on pharmacokinetic data of beta‐lactams.7 Figure is not an output of ATLAS mPBPK but was constructed by processing ATLAS mPBPK simulation output using MATLAB.

Figure  5 shows further use of the ATLAS mPBPK platform as demonstrated by the fitting and simulation of two corticosteroids, methylprednisolone19 and dexamethasone.20 The intravenous kinetics of a 50‐mg/kg bolus of methylprednisolone in rats were well described using the proposed model structure (Figure  5 a), with clearance from the plasma compartment (CLnh) and K P being estimated with reasonable precision. Both parameters were in good agreement with the compartmental estimates of central clearance and BW‐normalized volume of distribution at steady state (V ss/BW).19 Simulation of the plasma‐concentration time profile of 2.25 mg/kg subcutaneous dexamethasone in rats using parameter values from Song et al.20 yielded good recharacterization of the data (Figure  5 b).

Figure 5

Demonstration of the ATLAS mPBPK tool based on the fitting of methylprednisolone parameters to pharmacokinetic data and simulation of dexamethasone. (a) The intravenous kinetics of a 50 mg/kg bolus of methylprednisolone in rats were well described using the proposed model structure, with clearance from the plasma compartment (CLnh) and the partition coefficient (K P) being estimated with reasonable precision. Both parameters were in good agreement with compartmental estimates of central clearance and BW–normalized volume of distribution at steady state (V ss/BW).19 (b) Simulation of the plasma‐concentration time profile of 2.25 mg/kg subcutaneous dexamethasone in rats using parameter values from Song et al.20 yielded good recharacterization of the data. C.I., confidence interval; LB, lower bound; mPBPK, minimal physiologically‐based pharmacokinetic; UB, upper bound.

Large molecules

No target binding

The schematic of the mPBPK model for large molecules with no target binding is shown in Figure  2 b as previously published.3 The differential equations for the large‐molecule mPBPK model without target binding are: where C p is the concentration of the monoclonal antibody (mAb) in V p (plasma volume), C lymph is the concentration of mAb in V lymph , C tight and C leaky are mAb ISF concentrations in tissues V tight and V leaky. The V tight and V leaky are volumes of ISF in tissues that have continuous and discontinuous or fenestrated capillaries. Based on the report of Sarin,21 muscle, skin, adipose, and brain are tissues assigned to V tight, and all other tissues to V leaky (e.g., liver, kidney, heart, etc.). V lymph is lymph volume that is assumed equal to blood volume, L is total lymph flow that equals the sum of L 1 and L 2 that are lymph flows for V tight and V leaky, σ1 and σ2 are vascular reflection coefficients for V tight and V leaky, σL the lymphatic capillary reflection coefficient (σL = 0.2), and CLp is the plasma clearance. The tissue‐related physiological values are: where K p is the available fraction of ISF for mAb distribution, which is largely determined by antibody size, charge, structure, and other physiochemical properties. Given the similar size and structure of most mAbs, charge will be the primary factor influencing K p, which was designated as 0.8 for native immunoglobulin G subclass 1 (IgG1) and 0.4 for native immunoglobulin subclass 4 (IgG4) according to previous studies.22, 23 The ISF is the total volume of ISF (4.35 mL in mice), and L is the total lymph flow (0.12 mL/hour in mice).

Target binding in plasma (c‐TMDD)

Figure  2 c shows the mPBPK model with TMDD in plasma.4 The equations are the following: where C total is the total concentration of mAb in plasma, C tight and C leaky represent concentrations of mAb in the tissues V tight and V leaky, C lymph is the concentration of mAb in lymph Vlymph, R total refers to the total concentration of target, and AR is the concentration of the drug‐target complex. Rate constants are k syn for target synthesis and k deg for target degradation. The free mAb concentrations in plasma are the following:and the drug‐target complex concentration is:

The V tight, V leaky , L 1, and L 2 related components are shown in Eqs(11), (12), (13), (14). Finally, K ss is the steady‐state binding constant defined by Gibiansky et al.24:where k int is the antibody‐target complex internalization, and k on and k off are antibody‐receptor association and antibody‐receptor dissociation rate constants.

Target binding in peripheral tissues (p‐TMDD)

Figure  2 d shows the mPBPK model with TMDD in the ISF of tissues.4 The equations are the following: where A tight_total, A leaky_total are the total mass of mAb, C tight_free and C leaky_free are the free concentrations of mAb, R tight_total and R leaky_total are the total concentrations of target, and AR tight and AR leaky are the concentrations of drug‐target in V tight and V leaky. Other parameters are similar to those in Eqs. 15–22. Considering that TMDD is mostly associated with antibodies that bind with cell membrane receptors, only free mAb is assumed to be collected in lymph and further recycled back to plasma, and the drug‐receptor complex is immobile in ISF. The free antibody concentrations in V tight and V leaky are the following: and the drug‐target complex concentration in V tight and V leaky are the following:

Figure  6 demonstrates the operation of the ATLAS mPBPK large‐molecule model based on the PK data of Trastuzumab25, which binds in tissues and Mavrillimumab26, which binds in plasma. Parameter values were taken from ref. 4.

Figure 6

Demonstration of ATLAS mPBPK large‐molecule model based on PK data of Trastuzumab25 (target binding in peripheral tissues) and Mavrillimumab26 (target binding in plasma) for a range of doses. Figures are not an output of ATLAS mPBPK but were constructed by processing ATLAS mPBPK simulation output using MATLAB.  c‐TMDD, Target mediated drug disposition in plasma (central); hr, hour; mPBPK, minimal physiologically‐based pharmacokinetic modeling; p‐TMDD, target‐mediated drug disposition in peripheral tissues.

ATLAS mPBPK Parameter Estimation

After the user chooses the parameters for estimation (Figure  1 c, Estimate checkboxes) and inserts their lower and upper bounds (Figure  1 c, LB, UB editboxes), parameter estimation can be executed. Obviously, this requires that the user has previously inserted the PK data of interest for which the parameters will be optimized. Parameter estimation is performed using the nonlinear least‐squares solver of MATLAB “lsqnonlin” and setting the upper and lower bounds equal to those provided by the user (LB, UB). No weighting is added in the parameter estimation. For the 95% confidence interval calculations, the MATLAB function “nlparci” was employed, which  uses the best estimates and residuals and the Jacobian matrix outputs of the “lsqnonlin” function to estimate the Wald (or normal) confidence intervals. The 95% confidence interval of a parameter p is given by the following:where is the optimal parameter value resulting from least squares, tinv (0.975,df) is the Student's t inverse cumulative distribution function for 95% probability, df degrees of freedom (number of data—number of parameters), and diag(v) is the diagonal of the coefficient variance matrix calculated as the following:where J is the Jacobian matrix resulting from least squares, exponent T represents the transpose matrix, and σ 2 is the variance of the residuals. The variance of the residual σ 2 is calculated as the following:where norm () is the Euclidean norm and r the residuals.

ATLAS mPBPK SA

The user can choose the parameters for the evaluation of sensitivity (Figure  7 a, Sensitivity checkboxes) and run a local SA (Figure  7 a). In the SA performed in ATLAS mPBPK, the chosen parameters are varied by 10% of the initial values that are defined in the Value editboxes (Figure  7 a) while the other parameters are kept constant. Next, the sensitivity coefficients (s i,j) for area under the curve (AUC) and maximum concentration (Cmax) are calculated as follows:where ∂fi is the AUC or Cmax difference resulting from simulation with typical parameter values, and simulation with 10% variation to parameter p j and ∂p j is the difference between varied and nominal parameter values. SA produces a figure with two bar graphs (Figure  7 b) indicating the sensitivity coefficients of AUC and Cmax (Eq. (20)).

Figure 7

Sensitivity analysis of Trastuzumab 4 mg/kg.25 (a) ATLAS mPBPK interface showing Trastuzumab simulation and fitting. Red boxes indicate parameters chosen to test their sensitivities. (b) Sensitivity analysis output. Upper panel shows the sensitivity of AUC to the chosen parameters and lower panel the Cmax sensitivity to these parameters. C.I., confidence interval; LB, lower bound; mPBPK, minimal physiologically‐based pharmacokinetic; UB, upper bound.

Summary

One approach to reduce whole‐body PBPK model dimensionality and complexity is to lump some of the physiological compartments of the whole‐body model into a general group of tissues that share certain characteristics.27 Proper lumping enables PK modeling with acceptable loss of the underlying physiological information. Recently, Cao et al.3, 4, 7 and Cao and Jusko3, 4, 7 introduced a number of generic frameworks of mPBPK models that adequately describe the kinetics of a large group of small and large molecules. So far these models have been used to explain different scenarios such as interspecies scaling,8, 28 modeling sex differences in PK, and combined drug effects29, 30, 31 as well as in combination with pharmacodynamic models to investigate drug effects.20

In this tutorial, we present a MATLAB‐based tool for the modeling and simulation of mPBPK models (ATLAS mPBPK). The tool gives users the opportunity to run simulations, parameter estimation, and SA of a number of predefined mPBPK models for small and large molecules in an easy and efficient manner. Currently, the tool does not provide a framework through which the users can implement their own equations; to do so, the user must edit the MATLAB code.

ATLAS mPBPK is released as an open‐source project, and users can download it for free from the SourceForge repository (https://sourceforge.net/projects/atlas-mpbpk/files/ATLAS_mPBPK_v.1/).17 Its use does not require previous MATLAB experience. Collectively, ATLAS mPBPK constitutes a useful addition to the open toolboxes available to quantitative as well as clinical pharmacologists.

Funding

Supported by National Institutes of Health Grant GM‐24211.

Conflict of Interest

The authors declared no competing interests for this work.

Figure S1. Representative ATLAS mPBPK code. (a) Code introducing the right‐hand side of differential equations for mPBPK model of small molecules (i.v. administration), and (b) code introducing the right‐hand side of differential equations for mPBPK model of large molecules (i.v. administration).

Click here for additional data file.