### Population Pharmacokinetics

Population pharmacokinetics is the study of the sources and correlates of variability in drug concentrations among individuals who represent the target patient population receiving clinically relevant doses of a drug of interest. The goal of population pharmacokinetic analysis to identify pathophysiologic factors that cause changes in the dose-concentration relationship and the extent of these changes so that, if such changes are associated with clinically significant shifts in the therapeutic index, dosage can be appropriately modified. Population pharmacokinetic models adhere into a hierarchical structure. At the initial stage, the relationship between concentration and time (pharmacokinetics) is modeled for an individual patient. At the second stage, pharmacokinetic parameters that define each of the individual patients' drug concentration profiles are assigned some distributional form, after accounting for relevant covariate information. A primary aim of a population analysis is to determine covariates that are important predictors of pharmacokinetic parameters. A Bayesian model then requires a third stage in which prior distributions are specified for the parameters defining the second-stage distributional form and the intra-individual variance parameters. Such a Bayesian model defines and estimates the variability that is observed both in individual concentrations and between different individuals' pharmacokinetic parameters. This model framework makes it is possible to determine appropriate patient-specific dosage regimens that ensure the attainment of desirable drug concentrations.

The estimation method most commonly used in population pharmacokinetics and nonlinear mixed effect modeling in general is based on a maximum likelihood approach. Maximum Likelihood (ML) is an alternative to the least squares objective function; it seeks to maximize the *likelihood* or *log*-likelihood function (or to minimize the negative log-likelihood function). In general terms, the likelihood function is defined as:

\text{L}=\text{F(Y,Model)}={\Pi}_{\text{i=1}}^{\text{n}}{\text{{p[y}}_{\text{i}}{\text{,ModelParameters(x}}_{\text{i}}\text{)]}}

The probability (now called *L*, the *likelihood*) is predicted in the sample data, given the respective regression model. Provided that all observations are independent of each other, this likelihood is the geometric sum (Π, across *i* = 1 to *n* cases) of probabilities for each individual observation (*i*) to occur, given the respective model and parameters (*θ*'s) for the *x* values. As it is customary to express this function as a natural logarithm, the geometric sum becomes a regular arithmetic sum (Σ, across *i* = 1 to *n* cases). The larger the likelihood of the model, the larger is the probability of the dependent variable values to occur in the sample and the better is the fit of the model to the data. If all assumptions for standard multiple regression are met, then the standard least squares estimation method will yield results identical to the maximum likelihood method. If the assumption of equal error variances across the range of the *x* variable(s) is violated, then the weighted least squares method will yield maximum likelihood estimates.

The typical structural model is chosen from one of several compartmental models which incorporate the route of administration as a fixed input into the model with certain assumptions (i.e., linear or zero order input). Compartmental models are, for the most part, empirical even though they may incorporate some mechanistic assumptions so they appear more realistic. Numerically, they are generally easier to handle as opposed to mechanistic models. Complex mechanistic and/or highly parameterized structural models can be accommodated as well of course. The prediction engine discussed herein is not limited by the nature of the model definition.

The framework for the mixed effect modeling approach to population pharmacokinetic analysis can be defined as follows: for i = 1, ... n individuals in a population of interest, let x_{ij}, j = 1, ... n_{j} represent the design points on which the y_{ij} responses are observed. In the pharmacokinetic (PK) setting, x_{ij} are typically the sampling time points and y_{ij} are the observed concentrations in the biologic matrix of interest (usually plasma or blood). Hence, the PK response can be described by

y_{ij} = *f*(*θ*
_{i}, x_{ij}) + *ε*
_{ij}

where the function *f* denotes the structural model; *θ*
_{i} is the p × 1 parameter vector for the i^{th} individual and *ε*
_{ij} are the independently and identically distributed (*i.i.d*.) error terms assumed to be normal random variables with a zero mean and a variance (*σ*
_{e}
^{2}) which may depend on the mean concentration. The *ε*
_{ij}'s account for the intraindividual variability and may incorporate model misspecification or other unresolved (or incorrect) error partitioning. In most population pharmacokinetic software the structural model is chosen from a library of compartmental models, expressed as a closed form system of equations or defined via differential equations. The probability density function which accounts for the within-individual variability only as

p(y_{ij} | *θ*
_{i}, x_{ij})

The intra-individual variation about the ith individual is defined when the distribution *ε*
_{i} of is specified. The second stage model defines the between-individual variability in the parameters as follows

*θ*
_{i} = *θ* + *η*
_{i}

where *θ* is the mean parameter vector for the population and *η*
_{i} are the individual deviations assumed to be *i.i.d*. and normal with zero mean vector and covariance matrix Ω. The expression of *θ*
_{i} shown (additive) is one of numerous ways that individual *θ*'s can be defined. In addition, the population *θ* can be expressed as a function of covariates (β_{i}). The covariate matrix Ω captures both the variance and covariance among the *η*'s. The density of the second stage model can then be defined as

p(*θ*
_{i} | *θ*, Ω, β_{i})

where β_{i} represents the individual patient covariate data (i.e., age, sex, race, etc.). The third stage of the mixed effect model approach would represent a Bayesian representation in which the model would contain the prior distributions of the population parameters as mentioned previously.

### Prediction Models for Toxicity and Adverse Events

PK/PD models can be developed to explore the relationship between drug exposure and observed toxicity as well. Endpoints can be expressed as a dichotomous categorical variable representing the occurrence of an adverse effect (AE, e.g., nausea, vomiting) or drug reaction (1 = yes, 0 = no). The population PD data are viewed as a probabilistic outcome and analyzed using a logistic regression model [27]. The probability of an event for individual *i* at sampling time *j* is given by *p*
_{
ij
}. The ratio of the probability of that event occurring vs. the probability that it does not occur is given by the Odds Ratio; the log of the Odds Ratio is known as the logit function (*λ*
_{
ij
}):

\begin{array}{ll}{p}_{ij}=ex{p}^{{\lambda}_{ij}}/(1+ex{p}^{{\lambda}_{ij}})\hfill & Probability\hfill \\ \frac{{p}_{ij}}{(1-{p}_{ij})}=ex{p}^{{\lambda}_{ij}}\hfill & OddsRatio\hfill \\ log\left[\frac{{p}_{ij}}{(1-{p}_{ij})}\right]={\lambda}_{ij}\hfill \end{array}

Independent variables and covariates (*x*
_{
ij
}) will be incorporated into the model via the logit function, with population typical population (*θ*) and individual random effect (*η*
_{
i
}) parameters to be estimated:

*λ*
_{
ij
}= *f*(*x*
_{
ij
}, *θ*, *η*
_{
i
})

Covariate effects, and random effects can modulate the predicted probability in a positive or negative direction, with the probability constrained between the values of 0 and 1.

In these analyses it is expected that each individual will contribute only one observation for each outcome endpoint. Thus the individual random effects (*η*
_{
i
}) will be fixed at a value of zero and a naïve-pooled data analysis will be conducted. This approach has also been suggested as a check for nonlinear mixed-effects models of dichotomous outcome data[28]. The predicted likelihood (*l*
_{
ij
}for individual *i* and sampling time *j*) of the data (*y*
_{
ij
}), given the model and parameters will be described by a binomial probability density function:

{l}_{ij}={[{\mathrm{exp}}^{{\lambda}_{ij}}/(1+{\mathrm{exp}}^{{\lambda}_{ij}})]}^{{y}_{ij}}\cdot {[1/(1+{\mathrm{exp}}^{{\lambda}_{ij}})]}^{1-{y}_{ij}}

The likelihood for the entire population PD data set is simply the product of likelihoods across all individuals and data points. Diagnostic plots and the minimum value of the objective function are used to guide model building and assess goodness-of-fit.

### Methotrexate (MTX) Model

The administration of methotrexate (MTX) to children with cancer was chosen as our initial setting to develop the first drug dashboard prototype. The difficulty in effectively administering high-dose MTX to oncology patients lie in balancing efficacy and safety. Increased MTX exposure has been shown to be predictive of greater efficacy [29–31], while increased MTX concentrations and prolonged exposure time have also been linked to toxicity [32]. Due to the high inter- and intra-patient variability in methotrexate pharmacokinetics, monitoring of methotrexate plasma concentrations in individual patients has become a standard procedure in order to identify patients at risk of toxicity. Typically, patient plasma concentrations are monitored starting at 24 hours post infusion until MTX plasma concentrations fall below 0.1 to 1 μM [33–37], with adjunct rescue therapy implemented as needed.

The occurrence of methotrexate-induced renal toxicity further complicates chemotherapy administration. Although methotrexate-induced nephrotoxicity is a relatively rare occurrence, it is none-the-less a life threatening complication of methotrexate therapy [38]. Since methotrexate is mainly cleared from systemic circulation via glomerular filtration and renal secretion, delayed drug elimination is a product of this nephrotoxicity. This results in prolonged drug exposure and elevated plasma concentrations. As a result of this increased exposure, severe adverse events such as myelosuppression, mucositis, and hepatitis become more prevalent and severe.

Numerous studies have been conducted to examine the feasibility and reliability of applying Bayesian forecasting approaches to predicting MTX pharmacokinetics. The goals of these studies have been to predict MTX concentrations at later times or the time that MTX concentrations fall below a threshold value [39–41], MTX dose adjustment [42, 43], or providing guidance for rescue administration in the case of elevated MTX concentrations for prolonged time periods [44]. The Bayesian prediction models developed thus far have concentrated on those patients with normal renal function, and are not applicable in the case of severe renal dysfunction secondary to high-dose MTX administration.

We have developed a population pharmacokinetic model to implement as a Bayesian predictor of MTX concentrations in patients with normal renal function and MTX-induced renal dysfunction. Plasma concentrations from patients with normal renal function and patients with MTX-induced renal dysfunction were obtained from standard clinical monitoring. The model was constructed from methotrexate dosing histories and monitored drug concentrations in 240 patients. The original dataset contained 2176 observations covering a range of one to 56 observations per patient (an average of 9 observations per patient). The age range was from 1 to 80 years with a weight range of 6.6 to 157 kg. The gender distribution was approximately 48% male (52% female). Hence, our underlying patient diversity allowed us to include and consider relevant size and demographic dependencies. The model was developed using NONMEM version VI[45].

Methotrexate disposition is described by a two-compartment model with first-order elimination. Although MTX clearance changes over time in patients with renal dysfunction, clearance is approximated with a simple model defined by two different clearance distributions for the two populations. Inter-subject variability in PK parameters was expressed using an exponential error model:

{P}_{i}=\widehat{P}\mathrm{exp}({\eta}^{Pi})

where:

*P*
_{
i
}is the estimated parameter value for individual *i*

\widehat{P} is the typical population value (geometric mean) of the parameter

*η*
^{Pi}are individual-specific interindividual random effects for individual *i* and parameter *P* and are assumed to be distributed: *η* ~ *N(0, ω*
^{2}
*)* with covariance defined by the inter-individual covariance matrix Ω.

The residual error was described by an additive expression on a log-transformed scale (i.e., proportional error model):

\mathrm{ln}({C}_{ij})=\mathrm{ln}({\widehat{C}}_{ij})+{\epsilon}_{ij}

where:

*C*
_{
ij
}is the *j* th measured observation in individual *i*

{\widehat{C}}_{ij} is the *j* th model predicted value in individual *i*

*ε*
_{
ij
}is the additive residual random error for individual *i* and measurement *j* and is assumed to be independently and identically distributed

The Bayesian forecasting model utilizes the NONMEM PRIOR subroutine to incorporate population priors into the model. Fixed effects parameters obtained from the final pop PK model were implemented for the initial Bayesian model. Prior distributions of the fixed effects parameters were obtained from the variance-covariance matrix from the final pop PK model as well. Prior distributions for random effects parameters were specified as an inverse Wishart distribution. Clearance was implemented as a mixture model, where a patient is assigned to a population (normal or impaired clearance) based on the probability of that patient belonging to either population given their MTX plasma concentrations. The Bayesian forecasting model was evaluated using MTX plasma concentrations that were not used during model construction. The model reliably predicts future MTX plasma concentrations from two prior concentrations in all patients except a small number who develop renal toxicity at delayed times (> 48 hours). In these patients, the addition of a third concentration after 48 hours increases the precision of the prediction of concentrations at later times.

### Dashboard Design and System Architecture

The MTX dashboard was developed based on a three-tier architecture comprising a back end database tier, a business logic middle tier and a data presentation/user interface tier at the front (see Figure 2). The database tier consists of patient records from our electronic medical records system (Sunrise Clinical Manager, SCM) merged with data from patient registration system (IDX), lab data management system (Clarity) and adverse event management system (proprietary). Required data from these systems are extracted and loaded into relational tables within the staging area of the PKB. The data fields are then processed systematically for gaps and manually filled (from patient charts) when any of the missing data are critical for functionality. Such gaps are then noted for future improvements in the data collection process. The level of data validation is minimal and sufficient for testing; a comprehensive data validation approach has been outlined for implementation prior to production release. Various views and summary tables are created from the relational tables for quick retrieval by the application. The current version of the application retrieves data from staging area tables via views and summary tables. Eventually, a multidimensional data mart will permit ad hoc querying and drill down analysis. Upon its completion, data from the staging area will be transformed and loaded into the PKB datamart.

The middle tier consists of rules and processing logic required to collect and prepare data for user presentation (alerts, filters, aggregations, derived values and predictions). Predictions are conducted in an external computational platform – our modeling and simulation (M&S) workbench. This platform can execute code in a variety of languages provided they can run in a batch mode. Of note, the M&S workbench can currently accommodates many of the standard prediction engines used to forecast PK and PK/PD relationships (NONMEM, SAS, SPLUS and R). Details of analytical run processing using NONMEM with the workbench are described below. While the workbench can perform various data processing functions and analytics including generation of plots and figures, it is important to note that all PKB related analytics are gated in the middle tier through logic to ensure that minimally required data sets are available for each patient or sets of patients for meaningful analysis (e.g., appropriate data density to make predictions, etc).

The user interface is a currently web-based and utilizes a combination of HTML, JavaScript and XML content. We are in the process of migrating to an AJAX paradigm and implementing Web 2.0 standards (the new standard for live HTML content). Upon successful transition to Web 2.0, it is possible to readily support more user interface types such as java applet, java application or stay completely within the browser. The Workbench component can co-reside on the same server or installed on another server within the intranet or even hosted anywhere on the internet. The system is designed to transfer data and messages in an authenticated and completely secure manner. When a physician invokes a given forecasting function within PKB, the data required of that patient for the analysis of interest is extracted and written to a flat file in the required application format. The system then edits a previously prepared control file (with job settings and analytical inputs) with necessary changes and makes a call to the Workbench requesting an execution. The communication and data transfer mechanisms are implemented via message queues, thereby ensuring guaranteed delivery, fault tolerance and scalability. The Workbench execution module listening to the queue receives the request and submits the jobs to NONMEM or whatever application has been specified for that analysis by the experts and calls the process manager module to monitor for completion. Upon completion, the process manager calls an application specific parser module to parse and return the results back to the calling application (middle-tier). The middle-tier receives the results and formats it into an XML object and returns the same to the browser for display. If the output of the analysis calls for graphical display of results, the middle-tier sends the unformatted output to the graphing engine and forwards the generated graphs back to the browser for display. The graphs can be setup to support drill downs, mouseovers or invoke further analysis. This ensures that analysis of any kind can be performed via the M&S Workbench to support drug/disease dash boards of varying complexities.