# Clustering Plasma Concentration-Time Curves: Applications of Unsupervised
Learning in Pharmacogenomics

Jackson P. Lautier***Department of Mathematical Sciences, Bentley University,
Waltham, MA, USA Stella Grosser†††Office of Biostatistics, Center for Drug
Evaluation and Research, U.S. Food and Drug Administration, Silver Spring, MD,
USA Jessica Kim22footnotemark: 2 Hyewon Kim‡‡‡Office of Clinical Pharmacology,
Center for Drug Evaluation and Research, U.S. Food and Drug Administration,
Silver Spring, MD, USA Junghi Kim22footnotemark: 2 Corresponding to
<EMAIL_ADDRESS>10903 New Hampshire Avenue, Silver Spring, MD, USA

###### Abstract

Pharmaceutical researchers are continually searching for techniques to improve
both drug development processes and patient outcomes. An area of recent
interest is the potential for machine learning (ML) applications within
pharmacology. One such application not yet given close study is the
unsupervised clustering of plasma concentration-time curves, hereafter,
pharmacokinetic (PK) curves. In this paper, we present our findings on how to
cluster PK curves by their similarity. Specifically, we find clustering to be
effective at identifying similar-shaped PK curves and informative for
understanding patterns within each cluster of PK curves. Because PK curves are
time series data objects, our approach utilizes the extensive body of research
related to the clustering of time series data as a starting point. As such, we
examine many dissimilarity measures between time series data objects to find
those most suitable for PK curves. We identify Euclidean distance as generally
most appropriate for clustering PK curves, and we further show that dynamic
time warping, Fréchet, and structure-based measures of dissimilarity like
correlation may produce unexpected results. As an illustration, we apply these
methods in a case study with 250 PK curves used in a previous pharmacogenomic
study. Our case study finds that an unsupervised ML clustering with Euclidean
distance, without any subject genetic information, is able to independently
validate the same conclusions as the reference pharmacogenomic results. To our
knowledge, this is the first such demonstration. Further, the case study
demonstrates how the clustering of PK curves may generate insights that could
be difficult to perceive solely with population level summary statistics of PK
metrics.

Keywords: CYP2C19, distance metrics, hierarchical clustering, precision
medicine

##  1 Introduction

Plasma concentration time-curves or pharmacokinetic (PK) curves are generated
by plotting drug concentration levels in plasma samples at various time
intervals after the administration of a drug product (Shargel and Yu, 2016).
As such, they are an essential component of characterizing drug disposition,
which is an important prerequisite to determine or modify dosing regimens for
individuals and groups of patients (Shargel and Yu, 2016).  Because of recent
machine learning (ML) and artifical intelligence (AI) work related generally
to drug development (Zhang et al., 2022) and patient outcomes
(Thirunavukkarasu and Karuppasamy, 2022), it is natural to suppose that any ML
and AI applications with the potential to enhance the understanding or
interpretation of PK curves would be of significant interest within the
broader field of pharmaceutical research. Indeed, this is the case (e.g., Koch
et al., 2020; Zame et al., 2020; McComb et al., 2021). Absent from these
studies is the use of clustering techniques on data sets of PK curves,
however, and we thus present to our knowledge the first applied study of
clustering techniques for use on PK data.

We may use the observation that PK curves are time series data to utilize
existing literature related to the ML clustering of time series data as a
starting point to cluster similar PK curves. Broadly, time series clustering
is a technique for grouping time series data based on their similarity.
Clustering techniques for time series data have grown rapidly and have been
applied successfully in a wide range of domains including medicine (e.g.,
personalized drug design), environmental science, and many more (Javed et al.,
2020). In particular, there exists a large set of well-established time series
dissimilarity (or distance) measures (Montero and Vilar, 2014). It is not
straightforward to define dissimilarity between PK curves, and we illustrate a
potential shape-based distance interpretation for two PK curves with five
concentration sampling points in Figure 1. The novelty of this paper is
twofold. First, we narrow down the broad field of potential time series
dissimilarity measures to select a suitable one for use on PK curves.
Specifically, we find the following five dissimilarity measures most
applicable: correlation, Fréchet, dynamic time warping (DTW), temporal
correlation coefficient (CORT), and Euclidean. Of these five, we identify
Euclidean distance as generally most appropriate for clustering PK
concentration curves, and we summarize the pros and cons of all five of these
dissimilarity measures in Table 1.

Figure 1:  Illustrative Example of Distance Between Two PK Curves. It is not
straightforward to measure the dissimilarity between two PK curves for use
within a ML clustering exercise. The above is a geometric interpretation of
“distance” between two PK curves with five plasma concentration sampling time
points.

Second, we present a novel case study to illustrate the merits of clustering
PK curves. The case study objective is to use ML clustering to independently
analyze data from a concluded pharmacogenomic study that attempted to tailor
treatment strategies to the individual patient level. The case study data set
consists of PK curves, along with genetic and demographic information from 250
observations, and it spans nine Phase 1 studies. As such, the concentration
sampling time points vary between PK curve data observations. This is a well-
known issue in time series clustering (Montero and Vilar, 2014). While there
are dissimilarity measures in Section 2 that are able to handle PK curves with
a differing number of sampling times, we find such methods may not yield
desirable results (see Section 2 and the Appendix for details). Hence, we find
Euclidean distance using only the shared concentration sampling time points to
be most effective. Lastly, the case study demonstrates how clustering PK
curves with Euclidean distance can provide additional insights in comparison
to PK analysis based solely on PK parameters (e.g., area-under-the-time-
concentration curve (AUC), maximum concentration ($C_{\max}$), time-until-
maximum concentration ($T_{\max}$), etc.).

The paper proceeds as follows. The methodological review occurs in Section 2.
The case study then follows in Section 3, and Section 4 concludes. For
reference, the Appendix provides extended details on the dissimilarity
measures of Section 2, and the Supplementary Material provides additional
details on clustering methods, cluster selection criteria, and the case study
results.

##  2 Methods

Clustering algorithms are usually categorized by families, such as
hierarchical clustering, partition-based, model-based, and density-based
clustering (e.g., Hastie et al., 2009; James et al., 2013). We situate our
paper within hierarchical clustering with a particular focus on similarity
measures for PK curves and defer discussion of other algorithms to the
Supplemental Material. (Hierarchical clustering will also be used exclusively
within the case study, see Section 3.2.)

An important step in hierarchical clustering is to find an appropriate
distance or dissimilarity measure between data objects to be clustered, which
is done in Section 2.3. Section 2.4 then briefly reviews the well-known
problem of deciding on the final number of clusters (e.g., James et al.,
2013).

###  2.1 Hierarchical Clustering

Hierarchical clustering is a “bottom-up” clustering technique with an
attractive feature of producing a tree-based representation of observations
called a dendrogram. Colloquially, data objects that are relatively similar
are to be grouped into the same cluster, while data objects that are
relatively dissimilar are to be grouped into separate clusters. The
dendrogram, therefore, represents the relationships of similarity among all
objects in a data set (see bottom of Figure 2).

At the beginning of the algorithm, each data object is treated as a single
cluster. The two most similar clusters are then merged into a new single
cluster, and this new cluster becomes an updated data object with a value that
is determined by averaging its now two members. There are other methods to
determine the new cluster value besides averaging, but we defer this
discussion at present for ease of exposition (the precise vernacular is
linkage, see James et al. (2013) for details). The merging process continues
until all original data objects eventually merge into a single cluster, as
represented by the very top of the dendrogram.

The dendrogram itself does not report an optimal number of clusters; it is
better thought of as a visualization of similarity (or dissimilarity) within a
data set given a particular measure of dissimilarity. (We discuss
dissimilarity measures more thoroughly in Section 2.3.). To interpret the
amount of similarity between two PK curves on a dendrogram, it is necessary to
find the vertical point where the two curves first fuse. It is an error to
associate horizontal proximity of two curves on the $x$-axis of the dendrogram
with similarity. For example, consider the labeled subjects X98, X100, and Y8
on the bottom of Figure 2. Despite the fact subjects X100 and Y8 are quite
close in terms of horizontal labels, they do not fuse until the very top of
the dendrogram. Thus, X100 and Y8 should be considered quite dissimilar, on a
relative basis, among all PK curves within the complete data set. On the other
hand, X98 and X100 are much further apart in horizontal labels than X100 and
Y8, but they fuse much sooner vertically. Hence, X98 and X100 should be
interpreted as relatively more similar than X100 and Y8. Finally, horizontal
ordering of labeled subjects has no bearing on the interpretation of the
clustering outcome; it is akin to the horizontal ordering of bars on a bar
chart.

Figure 2:  An Illustration of Hierarchical Clustering for PK Curves. (Top) An
example of two resulting groups of 100 simulated PK curves with the group mean
curves denoted by the dashed and solid black curves. The simulated PK curves
were generated using R Core Team (2022) with the package linpk (Rich, 2022).
(Bottom) The red squares indicate the two resulting clusters using Euclidean
distance (1) for a given scenario. Simulated subjects X98, X100, and Y8 are
labeled to demonstrate that it is a mistake to associate a proximity on the
horizontal axis (i.e., Y8 and X100) with similarity. That is, subjects X98 and
X100 fuse on the vertical axis much earlier than X100 and Y8 and are thus
relatively more similar than X100 and Y8. Indeed, even with two clusters
selected, X100 and Y8 are in separate clusters. The horizontal dashed line
indicates where a “cut” to the dendrogram tree would occur to create two
clusters (see Section 2.4 for more details).

###  2.2 An Illustrative Example

We demonstrate the potential effectiveness of hierarchical clustering for PK
curves with an illustrative example. Consider two PK curves, $C_{1}$ and
$C_{2}$, from a one compartment linear PK model, assuming first order
absorption and first order elimination after oral administration.
Specifically, we assume the first dosage is administered at $t=0$. That is,

$C(t)=\frac{Dk_{a}}{V_{c}(k_{a}-k_{el})}\bigg{[}\exp(-k_{el}t)-\exp(-k_{a}t)\bigg{]},$

where $D$ is the dosage, $V_{c}$ is the central volume, $k_{a}$ is the first
order absorption, and $k_{el}=\mathcal{K}/V_{c}$ (where $\mathcal{K}$ is the
clearance). The curve $C_{1}$ has parameters $\mathcal{K}_{1}=66.65$,
$V_{c,1}=189.84$, and $k_{a,1}=30.36$. The curve $C_{2}$ has parameters
$\mathcal{K}_{2}=66.67$, $V_{c,2}=93.42$, and $k_{a,2}=43.90$. Both curves
assume $D=100$. This implies $C_{\max,1}=0.5$ and $C_{\max,2}=1$ with
$\text{AUC}_{1}=\text{AUC}_{2}=1.5$.

We then introduce random errors to simulate 100 realizations for each curve.
The resulting curves are displayed in Figure 2 (top), with $C_{1}$ represented
by the solid line and $C_{2}$ represented by the dashed line. The shaded
regions around PK curves $C_{1}$ and $C_{2}$ represent the range of 100
simulated curves. From the plot of the simulated PK curves in Figure 2 (top),
we have two distinct PK curve groups. We assume each PK curve has 15
measurements to consider a time series data vector of length 15. We then
perform hierarchical clustering combined with Euclidean distance to
demonstrate this unsupervised ML clustering technique is capable of
identifying the two distinct groups. Specifically, in our hypothetical
simulation setting, we find that hierarchical clustering with Euclidean
distance obtains a mean accuracy rate of 99.9% with a standard deviation of
0.26% in 1,000 replicates. The mean Rand Index (Rand, 1971) was 0.998 (i.e., 1
is perfect, 0.50 indicates no contribution) with a standard deviation of
0.0052. An example of a resulting dendrogram appears in Figure 2 (bottom). We
will further discuss how to select the number of clusters in Section 2.4.

###  2.3 Measures of Dissimilarity

The concept of dissimilarity draws from the mathematical concept of distance
and includes some intuitive properties. For example, any measure of distance
should be nonnegative, and two identical inputs should result in a distance of
zero. It is also preferable that a distance function be both symmetric and
satisfy the triangle inequality. These ideas may be formalized (e.g., Rudin,
1976, Definition 2.15).

Montero and Vilar (2014) review a large set of well-established time series
dissimilarity measures. Different distance metrics are designed to measure
different attributes, and it is important to determine a suitable
dissimilarity measure for a given data set to avoid reaching misleading
conclusions. There are four broad categories of dissimilarity measures: model-
free, model-based, complexity-based, and prediction-based (Montero and Vilar,
2014). Model-based approaches assume the time series data has been generated
from a time-dependent parametric model, typically an AutoRegressive Integrated
Moving Average (ARIMA). ARIMA models assume a sequential time dependence
whereby the current realization of a time series can be explained by some
number of its past realizations. These models often rely on assumptions, such
as stationarity and independent, identically distributed Gaussian random
errors, that are generally not satisfied for PK curves (though the latter may
be appropriate for accounting for concentration measurement or sampling error,
an application outside the scope of our analysis). See Shumway and Stoffer
(2006) for additional details. Complexity-based dissimilarity can suffer from
interpretation challenges. Finally, clustering completed PK curves will
generally not require making predictions about future levels. We can thus
consider model-free distance metrics as the best potential candidates to
support the hierarchical clustering of PK curves.

The model-free dissimilarity measures category is itself extensive (Montero
and Vilar, 2014). One helpful step is the further grouping of model-free
dissimilarity measures into shape-based or structure-based (e.g., Lin and Li,
2009; Corduas, 2010). Shape-based dissimilarities are based on geometric
profile comparisons at a localized level. Alternatively, structure-based
dissimilarities attempt to compare underlying dependence structures. Figure 3
presents an illustration of the meaningful difference in results using shape-
or structure-based dissimilarity measures for PK curves. Of the five choices
we find most suitable for PK curves, Euclidean distance, Fréchet, and DTW are
shape-based model-free dissimilarity measures. Correlation is a structure-
based distance metric, and CORT attempts to straddle both shape and structure.
A summary of each of these five distance metrics, including a list of pros and
cons, may be found in Table 1. For ease of exposition, we will present
Euclidean distance formally within this section and defer details of Fréchet,
DTW, correlation, and CORT to Appendix. The many other model free
dissimilarity measures not suitable for PK curves fall outside our scope, and
thus further discussion has been omitted for brevity. Montero and Vilar (2014)
provide extended details for interested readers.

Figure 3:  Dissimilarity Measure Behavior Illustration: Shape vs. Structure. A
pharmacologist comparing PK curves by AUC and $C_{\max}$ would likely prefer
Subject II and Subject III to measure as relatively more similar than either
Subject II or Subject III with Subject I. This is the case with Euclidean,
Fréchet, DTW, and CORT. Conversely, the correlation-based distance evaluates
Subject I and Subject II as identical because Subject II’s concentration level
at each time point is four times that of Subject I and hence linearly
dependent. Further, Euclidean, Fréchet, and DTW consider Subjects II and III
as most similar, followed by Subjects I and III, with Subjects I and II as the
most dissimilar. However, CORT considers Subjects I and III as the most
dissimilar, because it considers the linear dependence between Subjects I and
II. Hence, it may not be appropriate to use a purely correlation-based measure
of dissimilarity to cluster PK curves by geometric shape. In another
application, such as attempting to identify PK curves from diluted samples,
however, using a structure-based distance (correlation, CORT) may have a more
desirable performance. All distance metric calculations were performed using
the R package TSclust (Montero and Vilar, 2014). For reference, we have
$\\{\text{AUC},C_{\max},T_{\max}\\}$ of $\\{84,16,4\\}$, $\\{334,64,4\\}$, and
$\\{246,62,4\\}$ for Subjects I, II, and III, respectively. The PK parameters
were calculated using the R package PKNCA (Denney et al., 2015).

Euclidean distance is a point-by-point comparative measurement summarized into
a single number. Despite its canonical familiarity, we present a formal
definition for completeness and practioner-orientated readers. Suppose
$\bm{x}=\\{x_{1},x_{2},\ldots,x_{t}\\}^{\top}$ and
$\bm{y}=\\{y_{1},y_{2},\ldots,y_{t}\\}^{\top}$ are two $t$-dimensional vectors
of PK measurements, where $x_{j}$ and $y_{j}$ are measured at the same time
for $j\in\\{1,\ldots,t\\}$. Euclidean Distance, also sometimes called
$L_{2}$-distance, is defined as

$d_{L_{2}}(\bm{x},\bm{y})=\sqrt{\sum_{i=1}^{t}(x_{i}-y_{i})^{2}}.$ (1)

Because of the simplicity of (1), Euclidean distance is relatively
straightforward to visualize and interpret as a measure of the difference in
geometric shape between two PK curves. For example, the Euclidean distance
between $\text{Curve}_{1}$ and $\text{Curve}_{2}$ in Figure 1 is simply
$d_{L_{2}}(\text{Curve}_{1},\text{Curve}_{2})=\sqrt{\sum_{i=1}^{t}d_{i}^{2}}.$
Additionally, we find Euclidean distance computationally inexpensive, stable,
and capable of clustering PK curves effectively. The major shortcoming of
Euclidean distance, however, is the data restrictions required to ensure the
point-by-point comparison is valid. Specifically, the samples of each PK curve
must be made at the same time points or, alternatively, only measurement times
shared by all PK curves in a data set should be used. It is not unreasonable
to suggest interpolation or extrapolation in such instances of varied sampling
times (e.g., Keogh and Ratanamahatana, 2005), but its effectiveness for PK
curves will require further study. At present, we avoid pursuing interpolation
techniques for fear of introducing a difficult to measure confounding factor
into our analysis, though we acknowledge its exploration may be a viable path
of future research.

Both DTW and Fréchet are distance metrics that work as searching algorithms to
minimize distance over a set of comparisons between two PK curves. Notably,
the algorithms of each allow for temporal distortion between PK curves, which
means concentration samples at different times may end up being directly
compared (see Figure 10 for an illustration). Because “time” is commonly one
of the important factors characterizing drug response (e.g., $C_{\max}$,
$T_{\max}$), we find this inherent property of DTW and Fréchet to be
potentially undesirable. Indeed, such flexibility can return unexpected
results in the context of PK curves, as we demonstrate in Figure 4 with a
hypothetical example. A major distinction between Euclidean distance and the
DTW and Fréchet distance metrics is that the latter two will compute a valid
distance between two PK curves with either differing concentration sampling
times or even a different total number of sampling times. Contrast this to
Euclidean distance, which is restricted equivalent to point-by-point
comparisons only. This apparent attractive property has received some
criticism (Ratanamahatana and Keogh, 2004), and we further find it is not
enough to overlook the potential concerns identified in Figure 4. Finally,
because of the exhaustive searching required in the Fréchet and DTW
dissimilarity measures, they are relatively more computationally expensive
than Euclidean distance, especially for lengthy time series vectors within a
large observation data set. More details about DTW and Fréchet are illustrated
in Appendix.

Figure 4:  Dissimilarity Measure Behavior Illustration: Fréchet, DTW. If a
pharmacologist is interested in $T_{\max}$ in addition to measures such as AUC
and $C_{\max}$, then it may be desirable to classify Subjects A and B as
relatively more dissimilar than Subjects A and C. Fréchet and DTW are
minimization algorithms that consider more than single point-by-point
comparisons, however, and thus they each identify Subjects A and B as most
similar. Euclidean, CORT, and the correlation-based distance metric each
consider Subjects A and C as most similar. Because Subject C is a twice
multiple of Subject A, the correlation-based distance considers the two curves
identical, while CORT also considers Subject A and C as relatively more
similar than Subjects A and B or Subjects B and C. All distance metric
calculations were performed using the R package TSclust (Montero and Vilar,
2014). For reference, we have $\\{\text{AUC},C_{\max},T_{\max}\\}$ of
$\\{37,15,3\\}$, $\\{38,15,5\\}$, and $\\{74,30,3\\}$ for Subjects A, B, and
C, respectively. The PK parameters were calculated using the R package PKNCA
(Denney et al., 2015).

The lone structure-based dissimilarity measure we find applicable to cluster
PK curves is based on Pearson’s correlation coefficient (Golay et al., 1998),
$d_{COR}$, which we define formally in (5). The distance $d_{COR}$ considers a
linear relationship between two PK curves in evaluating dissimilarity. For
example, in Figure 3, $d_{COR}$ will recognize that Subjects A and C are
perfectly linearly correlated and treat them as “similar”, despite the notable
geometric differences. Some positive attributes of $d_{COR}$ are that it
offers a straightforward interpretation and is computationally inexpensive. As
with Euclidean distance, it fails for two PK curves with a different number of
sampling times.

Finally, the dissimilarity measure CORT, $d_{CORT}$ is both shape-based and
structure-based. It’s formal definition may be found in (7). Of the five
dissimilarity measures we discuss, $d_{CORT}$ is the only one that
simultaneously considers both geometric shape and structure (Chouakria and
Nagabhushan, 2007). Computational time and the ability to interpret $d_{CORT}$
can vary depending on the choice of its embedded parameters. See the brief
discussion near (7) or Chouakria and Nagabhushan (2007) for details.

To close this subsection, it is worth remarking that the four additional
model-free distance metrics we identify as suitable for clustering PK curves –
DTW, Fréchet, correlation, CORT – also perform with a similar level of
accuracy as that of Euclidean distance for the 200 simulated curves in the
illustrative example of Section 2.2 and Figure 2. For more detail on other
model-free dissimilarity measures and time series dissimilarity measures
overall (including the five we discuss), see Montero and Vilar (2014). We
emphasize that of all the additional model-free dissimilarity measures in
Montero and Vilar (2014), we believe that none outside of the five we
summarize in Table 1 are suitable for clustering PK curves. For reference, the
Appendix provides additional dissimilarity measure details.

Table 1:  Summary of dissimilarity measures suitable for PK curves. The
specific formula references are Euclidean (1), Fréchet (4), DTW (3),
correlation (5), and temporal (7). Type refers to shape- versus structure-
based (see Section 2.3). The pros and cons are based on observations of the
authors in utilizing the dissimilarity measures within a ML clustering
application of PK curves.

Distance | Type | Formula |  Pros |  Cons
---|---|---|---|---
Euclidean | shape | $d_{L_{2}}(\bm{x},\bm{y})$ |  • Interpretability • Computation time • Stability • Performance |  • Only valid if two PK curves share the same measurement times
Fréchet | shape | $d_{F}(\bm{x},\bm{y})$ |  • Valid if two PK curves have different measurement times • Valid if two PK curves have a different number of measurement times |  • Interpretability • Computation time • Potentially unexpected results for PK curves
Dynamic Time Warping (DTW) | shape | $d_{DTW}(\bm{x},\bm{y})$ |  • Valid if two PK curves have different measurement times • Valid if two PK curves have a different number of measurement times • Many open-source tools • Many computational options |  • Interpretability • Computation time • Potentially unexpected results for PK curves
Correlation | structure | $d_{COR}(\bm{x},\bm{y})$ |  • Interpretability • Computation time • Identifies dependence (linear) between PK curves |  • Not suitable to cluster based on the geometric shape of PK curves • Only valid if two PK curves have the same number of sampling times
Temporal | both | $d_{CORT}(\bm{x},\bm{y})$ |  • Only metric to consider both shape and structure of PK curves |  • Only valid if two PK curves have the same number of sampling times • Requires more user choice, such as $\phi_{k}$ or $d(\bm{x},\bm{y})$

###  2.4 Selecting the Number of Clusters

The last decision for cluster analysis is to select the appropriate number of
clusters. There is no current statistical consensus criteria to arrive at a
definitive validation of the resulting clusters (James et al., 2013). As such,
deciding on the optimal number of clusters often relies on interpretations
that may vary with the nature of the underlying data or even the desired
clustering resolution of the application.

For hierarchical clustering, a visual inspection of the dendrogram is
generally used (e.g., Figure 2) as a starting point. More specifically, we
decide where to “cut” the tree; i.e., draw a horizontal line that intersects
with the dendrogram branches. Each intersection between the horizontal line
and the dendrogram branches indicates a cluster (i.e., all observations below
the intersection point belong to that cluster). This is illustrated in Figure
2 and discussed more thoroughly in James et al. (2013, Figure 10.9, pg. 392).

Metrics, referred to as cluster validation indices (CVIs), are also available
for estimating the quality of partitions produced by clustering algorithms and
for determining the number of clusters in data. One such commonly used CVI is
the Calinski-Harabasz index (Arbelaitz et al., 2013), which is done by
comparing the average distance within a cluster against the average distance
between clusters. The Calinski-Harabasz is calculated by comparing the average
distance within a cluster against the average distance between clusters. Let
$N$ be the number of PK curves with $t$ measurement times, denoted by a vector
$\bm{x}_{i}=(x_{1},\ldots,x_{t})^{\top}$, $i=1,\ldots,N$. Further assume there
are $K$ clusters, each of which is denoted $c_{j}$, $j=1,\ldots K$. If we let
$n_{j}$ denote the number of PK curves in cluster $c_{j}$, $j=1,\ldots,K$,
then we can define the centroid of a cluster as its mean vector,

$\bar{c}_{j}=\frac{1}{n_{j}}\sum_{\bm{x}_{i}\in c_{j}}\bm{x}_{i}.$

Similarly, the centroid of the entire PK data set is

$\bar{\bm{x}}=\frac{1}{N}\sum_{i=1}^{N}\bm{x}_{i}.$

The Calinski-Harabasz index for this particular set of clusters,
$C=\\{c_{1},\ldots,c_{K}\\}$, denoted $\text{CH}(C)$, is then

$\text{CH}(C)=\frac{(N-K)\sum_{c_{k}\in
C}n_{j}d_{L_{2}}(\bar{c}_{k},\bar{\bm{x}})}{(K-1)\sum_{c_{k}\in
C}\sum_{\bm{x}_{i}\in c_{k}}d_{L_{2}}(\bm{x}_{i},\bar{c}_{k})}.$ (2)

In words, the Calinski-Harabasz CVI is a ratio of the separation between
clusters divided by the cohesion or similarity within a cluster. A user would
then select the number of clusters that generates the highest value of
$\text{CH}(C)$. We present an example application with PK curve data in
Section 3.2. Other CVIs and techniques are described in Supplementary Material
S.2 and the references (e.g., Rousseeuw, 1987; Kim and Ramakrishna, 2005;
Suzuki and Shimodaira, 2006; Saitta et al., 2007; Arbelaitz et al., 2013).

##  3 Case Study

Formal evaluations by the U.S. Food and Drug Administration (FDA) have
concluded there is sufficient evidence to suggest that subgroups of patients
with certain genetic variants, or genetic variant-inferred phenotypes are
likely to have an altered drug metabolism or differential therapeutic effects.
For example, the FDA recommends a reduced dosage for patients classified as a
poor metabolizer (PM) in many therapeutic areas (U.S. Food and Drug
Administration, 2021). Similar results may also be found in the literature
(e.g., Hicks et al., 2013, 2015; Lee et al., 2018; Lima et al., 2021).

In the original pharmacogenomic analysis of our case study data, subject
genetic information (and corresponding phenotypes) were used in conjunction
with subject PK curve samples to identify phenotype groups that had greater
exposure to the active ingredient after drug administration (as measured by
AUC). Hence, we may perform a validating case study analysis for the use of
hierarchical clustering with Euclidean distance on PK curve data by performing
an unsupervised ML clustering on the raw PK curve data alone (i.e., ignoring
the subject’s genetic information). In other words, if the methods of Section
2, using only the PK curve data, can independently recover the same
pharmacogenomic-based drug therapy conclusion as our reference study, then we
will have demonstrable evidence of the effectiveness of hierarchical
clustering with Euclidean distance on PK curve data. Indeed, we find not only
can we validate the pharmacogenomic-based drug therapy results with
unsupervised clustering on just PK curve data, but we also glean additional
insights that would not have been easily obtained without the use of our
proposed clustering approach. We elaborate as follows.

###  3.1 Data

Our data is compiled from nine Phase 1 studies in which we have each subject’s
genetic information of the cytochrome P450 (CYP) enzyme and individual PK
curve data. At the conclusion of the nine pharmacogenomic studies, it was
suggested to use a one-half dosage reduction for patients who are known to be
poor metabolizers (PMs) of the CYP enzyme based on an observed 2.3-fold
increase in drug exposure in PMs compared to that in non-PMs. The sample of
250 total observations consists of six PMs, 52 intermediate metabolizers
(IMs), 113 extensive metabolizers (EMs), 63 rapid metabolizers (RMs), and 16
ultra-rapid metabolizers (UMs). Table 2 summarizes basic PK curve metrics for
each phenotype and genotype. We can see the average dose-normalized
$\text{AUC}_{\text{last}}$ is noticeably higher in the PMs in comparison to
subjects with other metabolizer classifications, and conversely, the UMs have
the lowest dose-normalized $\text{AUC}_{\text{last}}$. The pattern of dose-
normalized $C_{\max}$ is similar. Figure 5 visualizes the average
concentration at each timepoint by metabolizer classification.

A complication with this data is that the sampling times of PK concentration
are not consistent for each of the nine Phase 1 studies. Therefore, without
any adjustments (or interpolation) to the data, only the distance metrics
Fréchet and DTW from Table 1 may be applied. The Fréchet and DTW distance
metrics are likely to distort time between concentration sampling points
between PK curves (e.g., Figures 4, 10 and Appendix), however, and may thus
produce results that are difficult to reconcile with PK principles. Given this
and because we desire to avoid the potential confounding effects of
interpolation, we proceed by using only the common sampling time points across
all 250 PK curve observations. Specifically, we find that all 250 subjects
share PK measurements at 0, 1, 2, 4, 6, and 8 hours, and we proceed using only
these shared concentration sampling time points. For completeness, we note
that DTW and Fréchet results on all time measurements were indeed
unsatisfactory in our testing, though the detailed results have been omitted
for brevity.

Table 2:  Summary of PK Metrics by Metabolizer Status. The maximum concentration ($C_{\max}$) and area-under-the-time-concentration Curve (AUC) are normalized by dosage. AUC is calculated until the last observation ($\text{AUC}_{\text{last}}$). A complete PK curve is equivalent to one observation and some subjects may have more than one PK curve. For reference, the phenotype abbreviations are Poor Metabolizer (PM), Intermediate Metabolizer (IM), Extensive Metabolizer (EM), Rapid Metabolizer (RM), and Ultra-Rapid Metabolizer (UM).  |  |  | $C_{\max}$ | $\text{AUC}_{\text{last}}$
---|---|---|---|---
CYP Enzyme |  | (ng/ML)/Dose(mg) | (ng*hr/ML)/Dose(mg)
Phenotype | Genotype | # Obs. | Mean | St.Dev. | Mean | St.Dev.
PM | *2/*2 | 6 | 9.03 | 1.64 | $42.09$ | $8.64$
IM | *1/*2 | 35 | 6.21 | 2.69 | $25.64$ | $11.18$
IM | *2/*17 | 17 | 6.03 | 2.17 | $26.49$ | $9.39$
EM | *1/*1 | 113 | 5.44 | 2.45 | $20.06$ | $9.63$
RM | *1/*17 | 63 | 4.88 | 2.42 | $18.49$ | $12.30$
UM | *17/*17 | 16 | 4.29 | 2.68 | $15.29$ | $8.77$
Total Sample |  | 250 | 5.46 | 2.55 | 21.11 | 11.32
Figure 5:  Average PK Curve by Metabolizer Classification. PK curves are
normalized by dosage. The CYP metabolizer categories are Ultra-rapid
Metabolizer (UM), Rapid Metabolizer (RM), Extensive Metabolizer (EM),
Intermediate Metabolizer (IM), and Poor Metabolizer (PM). A complete PK curve
is considered as one observation and some subjects may have more than one PK
curve. The Whisker plots show one standard deviation above and below the mean
of the data.

###  3.2 Results

We find unsupervised clustering on PK curve data alone, without using any
subject genetic information, independently validated the result of the
reference study analysis. Specifically, we find the hierarchical clustering
results with Euclidean distance indicate that a PM subject is likely to have
greater drug exposure (as measured by cluster-level AUC) than a subject
classified in the other metabolizer categories (IM, EM, RM, UM). This result
is consistent with the proposal of the original pharmacogenomic study which
suggests to use a one-half dosage reduction for patients of PMs of the CYP
enzyme.

To obtain these results, we employ hierarchical clustering to the data
reviewed in Section 3.1 by first ignoring each subject’s genetic information.
After the clusters have been established from just the PK curve data, we then
align the PK curve cluster index with each subject’s genetic information (or
corresponding metabolizer status). The optimal number of clusters is selected
based on a visual review of the dendrogram together with the Calinski-Harabasz
index ($\text{CH}(C)$). Figure 6 presents a dendrogram for the sample of 250
PK curves produced by Euclidean distance over the shared measurement times. We
also present a summary of the $\text{CH}(C)$ by number of clusters in Figure
7. Recall that $\text{CH}(C)$ should be maximized. We can see that
$\text{CH}(C)$ obtains a clear peak at four clusters. (We briefly note that it
is not uncommon to observe non-monotonic behavior of $\text{CH}(C)$ in other
studies (Litos et al., 2022).) Visually, the dendrogram of Figure 6 also
suggests that four clusters would be a reasonable cutting point. Therefore,
given that $\text{CH}(C)$ is maximized at four clusters and confirmation via a
visual inspection of the vertical fuse points of the dendrogram, we elect to
use four as the optimal number of clusters. The four clusters are also
indicated in Figure 6, and the six PMs have been labeled on the dendrogram for
easy identification. We also remark here that hierarchical clustering provides
a complete view of the structural heterogeneity of the data via the
dendrogram, which would not otherwise be easy to obtain from an analysis based
solely on summary PK metrics.

Figure 6:  Case Study Clustering Results: Dendrogram. The dendrogram for the 250 PK curve observations in Section 3.1 that results from clustering the distance matrix calculated by Euclidean distance over the shared measurement times. Our recommended four clusters are highlighted by the red squares. The cluster locations of the poor metabolizer (PM) subjects are indicated by the labels. All other metabolizer categories are not labeled. Table 3 presents summary PK metrics by cluster.  Figure 7:  Case Study Clustering Results: Calinski-Harabasz Index. The value of $\text{CH}(C)$ (2) for a potential number of clusters, $K$, for $K\in\\{1,\ldots,8\\}$ based on the dendrogram of Figure 6. We desire to maximize the $\text{CH}(C)$, which implies four clusters, but subjective factors should also be considered (i.e., the dendrogram, see Figure 6.) These results were calculated with the R package fpc (Hennig, 2023). Table 3:  Detailed Case Study Clustering Results. Summary statistics for $\text{AUC}_{\text{last}}$ and $C_{\max}$ represent the arithmetic average by cluster with the standard deviation (St.Dev.) by cluster provided in parenthesis. AUC is calculated until the last observation ($\text{AUC}_{\text{last}}$). Subject metabolizer status was not used to perform the clustering; only the PK curves were used to generate the clusters. For reference, the abbreviations are Poor Metabolizer (PM), Intermediate Metabolizer (IM), Extensive Metabolizer (EM), Rapid Metabolizer (RM), and Ultra-Rapid Metabolizer (UM).  Cluster | 1 | 2 | 3 | 4 | Total
---|---|---|---|---|---
# Obs. (%) | 204 (81.6%) | 39 (15.6%) | 3 (1.2%) | 4 (1.6%) | 250
| Distribution of Metabolizer Status by Cluster |
PM | 1 (16.67%) | 5 (83.33%) | 0 (0%) | 0 (0%) | 6 (100%)
IM | 35 (67.31%) | 16 (30.77%) | 0 (0%) | 1 (1.92%) | 52 (100%)
EM | 99 (87.61%) | 11 (9.73%) | 2 (1.77%) | 1 (0.88%) | 113 (100%)
RM | 55 (87.30%) | 5 (7.94%) | 1 (1.59%) | 2 (3.17%) | 63 (100%)
UM | 14 (87.50%) | 2 (12.50%) | 0 (0%) | 0 (0%) | 16 (100%)
| Avg. (St.Dev.) PK Metrics by Cluster |
$\text{AUC}_{\text{last}}$ (ng*hr/mL)/Dose(mg) | 17.12 (6.76) | 36.27 (8.28) | 60.01 (10.75) | 47.10 (10.69) | 21.11 (11.32)
$C_{\max}$ (ng/mL)/Dose(mg) | 4.66 (1.89) | 9.04 (2.17) | 10.22 (0.45) | 7.48 (1.16) | 5.46 (2.55)

In our proposed four clusters summarized in Table 3, a majority (204 out of
250) of observed PK curves are grouped into Cluster 1. Cluster 2 is also
relatively large with a membership of 39 observed PK curves. Both Clusters 3
and 4 are quite small, with 3 and 4 members, respectively. Figure 8 highlights
the PK curves that comprise each cluster. As shown in Figure 8, we have
performed the clustering based on the shape of the full PK curve without using
any other information about the observations or subjects. Recall our goal is
to see if we can independently recover the same potential of increased drug
exposure for PMs by using unsupervised ML on only subject PK curve data. To
this end, we connect the clustering labels with metabolizer status and PK
parameters.

Figure 8:  Case Study Clustering Results: Visualizing Each Cluster. A visual
representation of the four suggested clusters summarized in Figure 6 and Table
3.

The details are as follows. Of the 204 observations grouped into Cluster 1,
only 17% (1 of 6) are PM subjects. Compare this to 67% (35 of 52) of IM
observations, and over 87% for EM (99 of 113), RM (55 of 63), and UM (14 of
16) observations. From a statistical perspective, this suggests that it is
likely that IM, EM, RM, and UM subjects would follow the PK profile of Cluster
1, which has an average dose-normalized AUC of 17.12 (ng*hr/mL)/Dose(mg). Most
PM subjects are in Cluster 2, however (5 of 6; 83%). This suggests that it is
more likely that a PM subject PK curve would follow that of Cluster 2, which
has an average dose-normalized AUC of 36.26 (ng*hr/mL)/Dose(mg). In other
words, the hierarchical clustering analysis we performed using only PK curves
suggests that a PM subject is likely to have a PK curve with a dose-normalized
AUC of approximately double that of observations with other phenotypes. This
corresponds very closely to the suggestion of a dosage reduction for patients
with PM in the reference studies. We remark that a binomial-exact hypothesis
test (e.g., Conover, 1971, pg. 97-104) also supports that PM subjects likely
belong to Cluster 2 and the metabolizer categories (IM, EM, RM, UM) belong to
Cluster 1, the details of which have been omitted for brevity.

The unsupervised ML analysis also provides additional insights not reported in
the reference studies. For example, Figure 6 and Table 3 indicate seven
observations in Clusters 3 and 4 that have both a higher AUC and $C_{\max}$
than Cluster 2 and are not PM subjects: IM (1), EM (3), and RM (3). These
observations may represent an important subset of the population with a
stronger reaction to the active ingredient than the general population, or
they may simply represent measurement errors or outliers. In either case, the
granularity of hierarchical clustering suggests seven subjects worthy of
further investigation that were not identified in the reference studies.

Furthermore, while the percentage of PM subjects in Cluster 2 is markedly
higher than the percentage of PM subjects in Clusters 1, 3, or 4, there are
phenotypes of each other type of metabolizer category (IM, EM, RM, UM) that
are also clustered into Cluster 2, albeit in much smaller percentages versus
Cluster 1. A further study may reveal additional factors potentially affecting
absorption of the active ingredient, which may help direct the research of
precision medicine into alternative directions (e.g., Do IM subjects also
warrant caution or a reduced dosage in some cases?). Although not performed
here, the clustered PK curve data could also be used to search for patterns
beyond metabolizer status, such as subject weight, body mass index, age, or
sex. We find this case study reveals the clustering technique we employ can
identify additional heterogeneity among subject’s drug exposure. For the
benefit of future research, our case study process is summarized in Figure 9.
As a brief comment to close the section, we also considered the other distance
metrics described in Section 2 (DTW, Fréchet, correlation, CORT). Only DTW was
able to produce results in agreement with Euclidean distance. For additional
detail, see the Supplemental Material.

Figure 9:  Suggested PK Curve Clustering Workflow.

##  4 Discussion

While ML and AI has attracted attention in the clinical research area in
recent years, clustering techniques have not been studied or applied in
understanding or interpretation of PK curves. The purpose of this paper was to
examine the applicability and potential of clustering algorithms in the
context of pharmacokinetics. We believe our paper contributes to the
literature on the use of ML in drug development. Specifically, we showed that
the use of hierarchical clustering on PK curve data with Euclidean distance
over the shared concentration measurement times is a potentially valuable
enhancement to pharmacological studies.

Important outcomes of our study were twofold. First was to recognize that PK
curve data is a time series data object and subsequently determine an optimal
measure of PK curve similarity. This required narrowing down a long list of
well-established time series dissimilarity measures (e.g., Montero and Vilar,
2014) down to five candidates: correlation, Euclidean, Fréchet, DTW, and CORT.
Of these five, we found that structure-based measures (i.e., correlation and
CORT) may produce undesirable results if we are not interested in any
potential linear relationship between PK curves (e.g., Figure 3). Further, we
found the potential temporal distortion of DTW and Fréchet may also produce
potentially misleading or difficult to interpret results within the context of
PK curves (e.g., Figure 4, Figure 10). We therefore suggest Euclidean
distance, which performs well, has low computational cost, and is easy to
interpret. For convenience, Table 1 provides a complete summary of our
observations related to each of these five potential distance metrics.

Second, we present the first known case study using unsupervised learning on
PK curve data. Specifically, we utilized PK curve data from a pharmacogenomic
study that suggested a reduced dosage reduction for patients known to be PMs
of the CYP enzyme. We found that clustering PK curve data with hierarchical
clustering based on Euclidean distance, while ignoring subject genetic
information, was able to independently arrive at the same conclusion: PMs
likely need a reduced dosage based on an elevated drug exposure. Additionally,
our clustering results found areas of potential further study, which were not
identified in the standard PK metric analysis reference study.

We conclude with brief comments on suggested future work. Despite our
encouraging results, a shortcoming of our approach is that Euclidean distance
requires a one-to-one correspondence between concentration sampling time
points. As we found, this may require removing some concentration measurements
from PK curve data, especially if they are compiled from different clinical
trial phases. Keogh and Ratanamahatana (2005) suggest interpolation, but this
will require further study for PK curves. DTW and Fréchet allow for
comparisons of time series of different lengths, but the possible temporal
distortion was not desirable in our application. Therefore, we suggest future
research related to finding dissimilarity measures that are suitable for PK
curves with irregular or a different number of concentration sampling time
points. At present, we leave this question open to future study, and suggest
the methods we studied herein as a current best practice.

##  5 Acknowledgements

This project was supported in part by an appointment to the Research
Participation Program at the U.S. Food and Drug Administration by the Oak
Ridge Institute for Science and Education (ORISE) through an interagency
agreement between the U.S. Department of Energy and the U.S. Food and Drug
Administration. Jackson Lautier’s contribution was supported in part by the
Center for Drug Evaluation and Research Intramural Funding Program and a
National Science Foundation Graduate Research Fellowship under Grant No. DHE
1747453. The simulation studies, data analyses, and manuscript drafts were
completed during Jackson Lautier’s ORISE training program at the FDA. We
further thank seminar participants at the 2022 Joint Statistical Meetings in
Washington, DC for providing helpful comments.

##  6 Declaration of Interest

The present study reflects the views of the authors and should not be
construed to represent the views or recommendations of the U.S. Food and Drug
Administration. The authors have no other potential conflicts of interest to
report.

## References

  * Arbelaitz et al. (2013) O. Arbelaitz, I. Gurrutxaga, J. Muguerza, J. M. Pérez and I. Perona (2013).  “An extensive comparative study of cluster validity indices.”  Pattern Recognition 46, 243–256.
  * Berndt and Clifford (1994) D. J. Berndt and J. Clifford (1994).  “Using dynamic time warping to find patterns in time series.”  In Proceedings of the 3rd International Conference on Knowledge Discovery and Data Mining, AAAIWS’94, p. 359–370. AAAI Press.
  * Chouakria and Nagabhushan (2007) A. D. Chouakria and P. N. Nagabhushan (2007).  “Adaptive dissimilarity index for measuring time series proximity.”  Advances in Data Analysis and Classification 1, 5–21.
  * Conover (1971) W. Conover (1971).  Practical Nonparametric Statistics.  New York: John Wiley & Sons.
  * Corduas (2010) M. Corduas (2010).  “Mining time series data: A selective survey.”  In F. Palumbo, C. N. Lauro and M. J. Greenacre (eds.), Data Analysis and Classification, pp. 355–362. Berlin, Heidelberg: Springer Berlin Heidelberg.
  * Denney et al. (2015) W. S. Denney, S. Duvvuri and C. Buckeridge (2015).  “Simple, automatic noncompartmental analysis: The PKNCA R package.”  Journal of Pharmacokinetics and Pharmacodynamics 42, 11–107.  R package version 0.10.1.
  * Eiter and Mannila (1994) T. Eiter and H. Mannila (1994).  “Computing discrete Fréchet distance.”  Technische Universität.
  * Giorgino (2009) T. Giorgino (2009).  “Computing and visualizing dynamic time warping alignments in R: The dtw package.”  Journal of Statistical Software 31, 1–24.
  * Golay et al. (1998) X. Golay, S. Kollias, G. Stoll, D. Meier, A. Valavanis and P. Boesiger (1998).  “A new correlation-based fuzzy logic clustering algorithm for fMRI.”  Magnetic Resonance in Medicine 40, 249–260.
  * Hahsler et al. (2019) M. Hahsler, M. Piekenbrock and D. Doran (2019).  “dbscan: Fast density-based clustering with R.”  Journal of Statistical Software 91, 1–30.
  * Hastie et al. (2009) T. Hastie, R. Tibshirani and J. Friedman (2009).  The Elements of Statistical Learning: Data Mining, Inference, and Prediction, Second Edition.  Springer New York, NY.
  * Hennig (2023) C. Hennig (2023).  fpc: Flexible procedures for clustering.  R package version 2.2-10.
  * Hicks et al. (2015) J. K. Hicks, J. R. Bishop, K. Sangkuhl, D. J. Müller, Y. Ji, S. G. Leckband, J. S. Leeder, R. L. Graham, D. L. Chiulli, A. LLerena, T. C. Skaar, S. A. Scott, J. C. Stingl, T. E. Klein, K. E. Caudle and A. Gaedigk (2015).  “Clinical Pharmacogenetics Implementation Consortium (CPIC) guideline for CYP2D6 and CYP2C19 genotypes and dosing of selective serotonin reuptake inhibitors.”  Clinical Pharmacology and Therapeutics 98, 127 – 134.
  * Hicks et al. (2013) J. K. Hicks, J. J. Swen, C. F. Thorn, K. Sangkuhl, E. D. Kharasch, V. L. Ellingrod, T. C. Skaar, D. J. Müller, A. Gaedigk and J. C. Stingl (2013).  “Clinical Pharmacogenetics Implementation Consortium (CPIC) guideline for CYP2D6 and CYP2C19 genotypes and dosing of tricyclic antidepressants.”  Clinical Pharmacology and Therapeutics 93, 402 – 408.
  * James et al. (2013) G. James, D. Witten, T. Hastie and R. Tibshirani (2013).  An Introduction to Statistical Learning with Applications in R.  Springer New York, NY.
  * Javed et al. (2020) A. Javed, B. S. Lee and D. M. Rizzo (2020).  “A benchmark study on time series clustering.”  Machine Learning with Applications 1, 1–13.
  * Keogh and Ratanamahatana (2005) E. Keogh and C. Ratanamahatana (2005).  “Exact indexing of dynamic time warping.”  Knowledge and Information Systems 7, 358–386.
  * Kim and Ramakrishna (2005) M. Kim and R. Ramakrishna (2005).  “New indices for cluster validity assessment.”  Pattern Recognition Letters 26, 2353–2363.
  * Koch et al. (2020) G. Koch, M. Pfister, I. Daunhawer, M. Wilbaux, S. Wellmann and J. E. Vogt (2020).  “Pharmacometrics and machine learning partner to advance clinical data analysis.”  Clinical Pharmacology & Therapeutics 107, 926–933.
  * Lee et al. (2018) C. R. Lee, V. B. Sriramoju, A. Cervantes, L. A. Howell, N. Varunok, S. Madan, K. Hamrick, M. J. Polasek, J. A. Lee, M. Clarke, J. D. Cicci, K. E. Weck and G. A. Stouffer (2018).  “Clinical outcomes and sustainability of using CYP2C19 genotype-guided antiplatelet therapy after percutaneous coronary intervention.”  Circulation: Genomic and Precision Medicine 11, e002069.
  * Lima et al. (2021) J. Lima, C. Thomas, J. Barbarino, Z. Desta, S. V. Driest, N. E. Rouby, J. Johnson, L. Cavallari, V. Shakhnovich, D. Thacker, S. Scott, M. Schwab, C. Uppugunduri, C. Formea, J. Franciosi, K. Sangkuhl, A. Gaedigk, T. Klein, S. Roseann and F. Takahisa (2021).  “Clinical Pharmacogenetics Implementation Consortium (CPIC) guideline for CYP2C19 and proton pump inhibitor dosing.”  Clinical Pharmacology and Therapeutics 109, 1417 – 1423\.
  * Lin and Li (2009) J. Lin and Y. Li (2009).  “Finding structural similarity in time series data using bag-of-patterns representation.”  In M. Winslett (ed.), Scientific and Statistical Database Management, pp. 461–477. Berlin, Heidelberg: Springer Berlin Heidelberg.
  * Litos et al. (2022) A. Litos, E. Intze, P. Pavlidis and I. Lagkouvardos (2022).  “Cronos: A machine learning pipeline for description and predictive modeling of microbial communities over time.”  Frontiers in Bioinformatics 2, 1–11.
  * McComb et al. (2021) M. McComb, R. R. Bies and M. Ramanathan (2021).  “Machine learning in pharmacometrics: Opportunities and challenges.”  British Journal of Clinical Pharmacology 88, 1482 – 1499\.
  * Montero and Vilar (2014) P. Montero and J. A. Vilar (2014).  “TSclust: An R package for time series clustering.”  Journal of Statistical Software 62, 1–43.
  * Mueen and Keogh (2016) A. A. Mueen and E. J. Keogh (2016).  “Extracting optimal performance from dynamic time warping.”  Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining pp. 2129–2130.
  * R Core Team (2022) R Core Team (2022).  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.
  * Rand (1971) W. M. Rand (1971).  “Objective criteria for the evaluation of clustering methods.”  Journal of the American Statistical Association 66, 846–850.
  * Ratanamahatana and Keogh (2004) C. Ratanamahatana and E. Keogh (2004).  “Everything you know about dynamic time warping is wrong.”  Proceedings of the Third Workshop on Mining Temporal and Sequential Data pp. 1–11.
  * Rich (2022) B. Rich (2022).  linpk: Generate Concentration-Time Profiles from Linear PK Systems.  R package version 1.1.2.
  * Rousseeuw (1987) P. J. Rousseeuw (1987).  “Silhouettes: A graphical aid to the interpretation and validation of cluster analysis.”  Journal of Computational and Applied Mathematics 20, 53–65.
  * Rudin (1976) W. Rudin (1976).  Principles of Mathematical Analysis.  McGraw-Hill, Inc.
  * Saitta et al. (2007) S. Saitta, B. Raphael and I. F. C. Smith (2007).  “A bounded index for cluster validity.”  In IAPR International Conference on Machine Learning and Data Mining in Pattern Recognition.
  * Sardá-Espinosa (2019) A. Sardá-Espinosa (2019).  “Time-series clustering in R using the dtwclust package.”  The R Journal 11, 22–43.
  * Schwarz (1978) G. Schwarz (1978).  “Estimating the dimension of a model.”  The Annals of Statistics 6, 461 – 464.
  * Scrucca et al. (2016) L. Scrucca, M. Fop, T. B. Murphy and A. E. Raftery (2016).  “mclust 5: Clustering, classification and density estimation using Gaussian finite mixture models.”  The R Journal 8, 289–317.
  * Shargel and Yu (2016) L. Shargel and A. B. Yu (2016).  Applied Biopharmaceutics & Pharmacokinetics, Seventh Edition.  McGraw Hill Education.
  * Shumway and Stoffer (2006) R. H. Shumway and D. S. Stoffer (2006).  Time Series Analysis and Its Applications with R Examples.  Springer Science+Business Media, LLC.
  * Suzuki and Shimodaira (2006) R. Suzuki and H. Shimodaira (2006).  “pvclust: An R package for assessing the uncertainty in hierarchical clustering.”  Bioinformatics 22, 1540–1542.
  * Thirunavukkarasu and Karuppasamy (2022) M. K. Thirunavukkarasu and R. Karuppasamy (2022).  “Forecasting determinants of recurrence in lung cancer patients exploiting various machine learning models.”  Journal of Biopharmaceutical Statistics pp. 1–15.
  * U.S. Food and Drug Administration (2021) U.S. Food and Drug Administration (2021).  “Table of pharmacogenetic associations.”  https://www.fda.gov/medical-devices/precision-medicine/table-pharmacogenetic-associations (Accessed 04-18-2023).
  * Zame et al. (2020) W. R. Zame, I. Bica, C. Shen, A. Curth, H.-S. Lee, S. Bailey, J. Weatherall, D. Wright, F. Bretz and M. van der Schaar (2020).  “Machine learning for clinical trials in the era of COVID-19.”  Statistics in Biopharmaceutical Research 12, 506–517.
  * Zhang et al. (2022) D. Zhang, J. Song, S. Dharmarajan, T. H. Jung, H. Lee, Y. Ma, R. Zhang and M. Levenson (2022).  “The use of machine learning in regulatory drug safety evaluation.”  Statistics in Biopharmaceutical Research pp. 1–5.

## Appendix: Measures of Dissimilarity

We provide additional details on the dissimilarity measures DTW, Fréchet,
correlation, and CORT discussed in Section 2.3. Let us first consider DTW
(Berndt and Clifford, 1994) in some depth, as it is a distance metric of great
popularity. For this exposition, we will use similar notation from Giorgino
(2009), which is a particularly helpful resource. Suppose we desire to compare
two time series, $\bm{x}=(x_{1},\ldots,x_{N})$ and
$\bm{y}=(y_{1},\ldots,y_{M})$. The counting symbol $i=1,\ldots,N$ will only be
used for indexing $\bm{x}$ and $j=1,\ldots,M$ for $\bm{y}$. Let
$d(x_{i},y_{j})$ be a local dissimilarity function between any pair of
elements $(x_{i},y_{j})$, $i=1,\ldots,N$ and $j=1,\ldots,M$. The DTW algorithm
rests on the warping curve, $\pi(k)$, $k=1,\ldots,T$. Let
$\pi_{x}(k)\in\\{1,\ldots,N\\}$ and $\pi_{y}(k)\in\\{1,\ldots,M\\}$. The
warping curve is then

$\pi(k)=\big{(}\pi_{x}(k),\pi_{y}(k)\big{)}.$

The warping functions, $\pi_{x}(k)$, $\pi_{y}(k)$, remap the indexes of
$\bm{x}$ and $\bm{y}$, respectively. Given $\pi$, we then calculate the
average accumulated distortion between the warped time series $\bm{x}$ and
$\bm{y}$,

$d_{\pi}(\bm{x},\bm{y})=\sum_{k=1}^{T}d(x_{\pi_{x}(k)},y_{\pi_{y}(k)})m_{\pi}(k)/M_{\pi},$

where $m_{\pi}(k)$ is a per-step weighting coefficient and $M_{\pi}$ is a
normalizing constant (the latter ensures the accumulated distortions are
comparable along different possible paths). A number of constraints are
imposed on $\pi$, see Giorgino (2009) for details. Finally, we have

$d_{DTW}(\bm{x},\bm{y})=\min_{\pi}d_{\pi}(\bm{x},\bm{y}).$ (3)

Informally, the goal of DTW is to pick the deformation of the time axes of
$\bm{x}$ and $\bm{y}$ which brings the two time series as close as possible to
each other. It is interesting that the elements $x_{i}$, $i=1,\ldots,N$ and
$y_{j}$, $j=1,\ldots,M$ only enter into the DTW computation through
$d(x_{i},y_{j})$; otherwise, DTW is just an indexing algorithm. It is typical
to choose Euclidean distance for $d(x_{i},y_{j})$. Further, we do not require
that $M=N$, and so time series of different lengths and measurement times may
be used. The results may not always align with expectations, however,
especially if a user drops the right endpoint constraint:
$\\{\pi_{x}(T)=N,\pi_{y}(T)=M\\}$.

To see this, consider Figure 10, which presents two hypothetical dosage-
normalized PK curves. The larger curve (circle points) is generated by
tenfolding the smaller curve (triangle points). In this case, statistical
normalization of each PK curve would not be appropriate (i.e., $z$-normalizing
would result in two identical time series, whereas we are interested in the
relative magnitude). This is counter to the standard advice in the literature
(Mueen and Keogh, 2016) and is a product of our specific application to
dosage-normalized PK curves.

Figure 10 shows an example mapping from the DTW algorithm. It’s clear how the
index remapping works to get the shortest distance. An analyst would need to
consider if such a remapping is theoretically appropriate for their PK curve
clustering application. Because $T_{\max}$ is an important factor in analyzing
drug response, it is possible such temporal distortion may be undesirable.
Though not pictured, we also analyzed the DTW algorithm with a subsample of
the smaller PK curve. The subsample has a final measurement at hour 48, and
the full sample has a final measurement at hour 72. As such, we relaxed the
end point constraint. This resulted in the algorithm assigning all mappings to
the time zero measurement in the larger PK curve. Once again, an analyst would
need to consider if such a remapping is theoretically appropriate for their PK
curve clustering application.

Figure 10:  Illustration of DTW Temporal Distortion. The larger PK curve
(circle points) has concentration measurements,
$\bm{x}=\\{0,2,4,6,9,10,7.5,5,2,0.5,0.1,0\\}$, at hours
$\\{0,0.5,1,2,3,4,8,12,24,36,48,72\\}$. The smaller PK curve (triangle points)
has concentration measurements $\bm{y}=\bm{x}/10$ at the same measurement
times. DTW exploits temporal distortion to minimize the distance between PK
curves (i.e., the time index is allowed to be flexible in the distance
calculation). For example, it matches the maximum point on the smaller curve
with many points on the larger curve (dashed lines). If $T_{\max}$ is an
important factor in characterizing drug response for a given application, this
distortion may be undesirable. The optimized index matching between the two
curves using DTW was found with the default settings of the R package dtw
(Giorgino, 2009).

We now briefly consider the Fréchet, correlation, and CORT dissimilarity
measures. The Fréchet (Eiter and Mannila, 1994) dissimilarity measure shares
some general properties with DTW in that it is an index-based searching
algorithm. Formally, let $\mathcal{Q}$ be the set of all possible sequences of
$q$ pairs that preserve the order of observations of the form
$\bm{r}=\big{(}(x_{a_{1}},y_{b_{1}}),\ldots,(x_{a_{q}},y_{b_{q}})\big{)}$,
with $a_{i}\in\\{1,\ldots,N\\}$ and $b_{j}\in\\{1,\ldots,M\\}$ such that
$a_{1}=b_{1}=1$, $a_{q}=N$, $b_{q}=M$, and $a_{i+1}=a_{i}$ or
$a_{i+1}=a_{i}+1$ and $b_{j+1}=b_{j}$ or $b_{j+1}=b_{j}+1$, for
$i\in\\{1,\ldots,N-1\\}$, $j\in\\{1,\ldots,M-1\\}$. Then Fréchet distance is
defined as

$d_{F}(\bm{x},\bm{y})=\min_{\bm{r}\in\mathcal{Q}}\bigg{(}\max_{i=1,\ldots,q}\lvert
x_{a_{i}}-y_{b_{i}}\rvert\bigg{)}.$ (4)

As with DTW, Fréchet distance may be calculated for PK curve time series of
different lengths.

For the remainder, assuming $\bm{x}$ and $\bm{y}$ have the same length, $t$.
That is, $\bm{x}=(x_{1},\ldots,x_{t})$ and $\bm{y}=(y_{1},\ldots,y_{t})$. The
correlation dissimilarity measure is based on Pearson’s correlation
coefficient (Golay et al., 1998):

$d_{COR}(\bm{x},\bm{y})=\sqrt{2\big{(}1-\text{COR}(\bm{x},\bm{y})\big{)}},$
(5)

where

$\text{COR}(\bm{x},\bm{y})=\frac{\sum_{i=1}^{t}(x_{i}-\bar{x})(y_{i}-\bar{y})}{\sqrt{\sum_{i=1}^{t}(x_{i}-\bar{x})^{2}}\sqrt{\sum_{i=1}^{t}(y_{i}-\bar{y})^{2}}}.$

To formally define CORT, we must first define a temporal correlation
coefficient

$\text{CORT}(\bm{x},\bm{y})=\frac{\sum_{i=1}^{t-1}(x_{i+1}-x_{i})(y_{i+1}-y_{i})}{\sqrt{\sum_{i=1}^{t-1}(x_{i+1}-x_{i})^{2}}\sqrt{\sum_{i=1}^{t-1}(y_{i+1}-y_{i})^{2}}}.$
(6)

The ratio in (6) ensures $\text{CORT}(\bm{x},\bm{y})\in[-1,1]$. The adaptive
dissimilarity measure, $d_{CORT}$, modulates between the raw values of each
vector $\bm{x}$ and $\bm{y}$ using the coefficient
$\text{CORT}(\bm{x},\bm{y})$. Formally,

$d_{CORT}(\bm{x},\bm{y})=\phi_{k}\big{(}\text{CORT}(\bm{x},\bm{y})\big{)}d(\bm{x},\bm{y}),$
(7)

where $\phi_{k}$ is an adaptive tuning function dependent on parameter $k$
that adjusts a conventional raw distance metric, $d(\bm{x},\bm{y})$, according
to the value of $\text{CORT}(\bm{x},\bm{y})$. Typically, $d(\bm{x},\bm{y})$ is
one of (1), (3), or (4). Both $d_{COR}$ and $d_{CORT}$ fail for two time
series of different lengths. For additional details, we suggest Montero and
Vilar (2014).

## Supplementary Materials

The following is intended as an online companion supplement to the manuscript,
Clustering plasma concentration-time curves: Applications of unsupervised
learning in pharmacogenomics. Please attribute any citations to the original
manuscript. This companion first presents additional details related to
Section 2. Specifically, Section S.1 briefly reviews other clustering
techniques beyond hierarchical clustering and Section S.2 expands on
techniques to select the number of clusters. This companion concludes by
presenting additional analysis related to the case study of Section 3.

###  S.1 Clustering Methods

For completeness, we briefly describe partition-based, model-based, and
density-based clustering techniques.

Partition-based clustering requires us to select the number of clusters, $K$,
beforehand, and each observation is then explicitly assigned to only one
cluster. Effectively, it is a combinatorial optimization problem to minimize
intra (within)-cluster distance while maximizing inter (between)-cluster
distance. The algorithm first assigns $K$ centroids, typically by randomly
selecting $K$ objects in the data set. These become the initial clusters, and
then a distance is calculated for all objects in the data and each centroid,
after which each data object is then assigned to its nearest centroid. The
centroids are then updated, and the process repeats iteratively until objects
no longer change clusters. A common example is the k-medoid algorithm along
with partition-around-medoids (PAM) centroid prototyping. Please see Sardá-
Espinosa (2019) for additional details.

For model-based clustering, a Gaussian mixture model is common, which assumes
each observation is specified through a multivariate Gaussian mixture
distribution of $G$ components. Each component of the mixture is typically
associated with a cluster. The task is then to fit each unknown parameter of
the underlying Gaussian distributions and $G$ through maximum likelihood
estimation. Additional details may be found in Scrucca et al. (2016).

Density-based clustering methods, conversely, do not assume parametric
distributions or consider variance. Therefore, density-based clusters may be
arbitrary in shape and can handle varying amounts of noise or outliers. Like
with Gaussian mixture model-based clustering, we do not need prior knowledge
of how to determine the number of clusters. A common technique is the DBSCAN
density-based clustering algorithm, which assigns observations to the same
cluster if they are density-reachable. For additional details, please consult
Hahsler et al. (2019).

The selection process for the number of clusters for the partition-based
clustering, Gaussian Mixture Model and DBSCAN algorithms is generally like
hierarchical clustering (see Section 2.4 and Supplementary Material S.2). In
some cases for the Gaussian Mixture Model, however, it is common to make the
selection based on the Bayesian-Information-Criterion (BIC) (Schwarz, 1978) in
addition to CVIs.

###  S.2 Cluster Selection Criteria

When selecting the number of clusters, there is often no clear choice. The
decision can rely on interpretations that may vary with the nature of the
underlying data or even the desired clustering resolution of the practitioner
(i.e., data management versus formal research). As mentioned in Section 2 and
illustrated in Section 3, many practitioners rely on a visual inspection of
the dendrogram in hierarchical clustering.

In addition, there is an automated bootstrapping method (Suzuki and
Shimodaira, 2006), in which resampling procedures are employed to assess the
potential sampling uncertainty of hierarchical clusters. The main idea of the
bootstrapping approach is to obtain replicates of the dendrogram by repeatedly
applying the cluster analysis to the resampled data. This allows a user to
identify the number of clusters that appear in a large percentage of
iterations.

Internal cluster validity indices (CVIs) attempt to evaluate a proposed number
of clusters based only on the data used to define the clusters. This is done
by comparing the average distance within a cluster against the average
distance between clusters. Ideally, the average distance within a cluster will
be small, indicating similar observations have been grouped, while the average
distance between clusters will be large indicating different observations have
been treated as so. Generally, the varying definitions of within and between
cluster distances are usually ratioed for each possible number of clusters.
This allows for a user to compare the CVI value for each possible number of
clusters. Arbelaitz et al. (2013) compares performance of 30 CVIs in many
different environments with different characteristics. There is a selection of
internal CVIs to be maximized: Silhouette (Rousseeuw, 1987), Dunn (Arbelaitz
et al., 2013), COP (Arbelaitz et al., 2013), Calinski-Harabasz (Arbelaitz et
al., 2013), Score Function (Saitta et al., 2007); or minimized: Davies-Bouldin
(Arbelaitz et al., 2013), Modified Davies-Bouldin (Kim and Ramakrishna, 2005).

###  S.3 Case Study: Additional Details

The following figures are a summary of results related to hierarchical
clustering of a distance matrix for the PK curve data and case study of
Section 3. In Section 3.2, a thorough discussion of hierarchical clustering
with Euclidean distance on the shared measurement times of the 250 PK curve
observations was presented. In this supplement, we present numerical results
for the other four measures of dissimilarity summarized in Table 1:
correlation, Fréchet, Dynamic time warping (DTW), and temporal (CORT). As with
the data of Section 3.2, each of Figures S.1, S.2, S.3, and S.4 is a summary
of clustering analysis performed on the shared measurement times of the 250 PK
curve observations. We can see that the choice of a distance metric can have a
large influence on the results. Of the four dissimilarity measures we
consider, only DTW performs in a way consistent with Euclidean distance.
Unless otherwise stated, the calculations were performed using the same
statistical packages referenced in the main manuscript.

Figure S.1:  Detailed Case Study Clustering Results: Correlation Distance. Summary statistics for $\text{AUC}_{\text{last}}$ and $C_{\max}$ represent the arithmetic average by cluster. AUC is calculated until the last observation ($\text{AUC}_{\text{last}}$). The number of clusters was selected by a visual inspection of the dendrogram in combination with the Calinski-Harabasz Index. For reference, the abbreviations are Poor Metabolizer (PM), Intermediate Metabolizer (IM), Extensive Metabolizer (EM), Rapid Metabolizer (RM), and Ultra-Rapid Metabolizer (UM). The cluster locations of the PMs are indicated in the dendrogram for reference.  | Dendrogram | Calinski-Harabasz (CH) Index by Cluster
---|---
|
|  Cluster | 1 | 2 | Total
---|---|---|---
# Obs. | 238 | 12 | 250
| Distribution of Metabolizer Status by Cluster |
PM | 6 (100%) | 0 (0%) | 6 (100%)
IM | 48 (92.31%) | 4 (7.69%) | 52 (100%)
EM | 107 (94.69%) | 6 (5.31%) | 113 (100%)
RM | 61 (96.83%) | 2 (3.17%) | 63 (100%)
UM | 16 (100%) | 0 (0%) | 16 (100%)
| Avg. PK Metrics by Cluster |
$\text{AUC}_{\text{last}}$ (ng*hr/mL)/Dose(mg) | 20.72 | 28.74 | 21.11
$C_{\max}$ (ng/mL)/Dose(mg) | 5.50 | 4.59 | 5.46
Figure S.2:  Detailed Case Study Clustering Results: Fréchet Distance. Summary statistics for $\text{AUC}_{\text{last}}$ and $C_{\max}$ represent the arithmetic average by cluster. AUC is calculated until the last observation ($\text{AUC}_{\text{last}}$). The number of clusters was selected by a visual inspection of the dendrogram in combination with the Calinski-Harabasz Index. For reference, the abbreviations are Poor Metabolizer (PM), Intermediate Metabolizer (IM), Extensive Metabolizer (EM), Rapid Metabolizer (RM), and Ultra-Rapid Metabolizer (UM). The cluster locations of the PMs are indicated in the dendrogram for reference.  | Dendrogram | Calinski-Harabasz (CH) Index by Cluster
---|---
|
|  Cluster | 1 | 2 | Total
---|---|---|---
# Obs. | 33 | 217 | 250
| Distribution of Metabolizer Status by Cluster |
PM | 4 (66.67%) | 2 (33.33%) | 6 (100%)
IM | 12 (23.08%) | 40 (76.92%) | 52 (100%)
EM | 12 (10.62%) | 101 (89.38%) | 113 (100%)
RM | 3 (4.76%) | 60 (95.24%) | 63 (100%)
UM | 2 (12.50%) | 14 (87.50%) | 16 (100%)
| Avg. PK Metrics by Cluster |
$\text{AUC}_{\text{last}}$ (ng*hr/mL)/Dose(mg) | 36.77 | 18.72 | 21.11
$C_{\max}$ (ng/mL)/Dose(mg) | 9.39 | 4.86 | 5.46
Figure S.3:  Detailed Case Study Clustering Results: DTW Distance. Summary statistics for $\text{AUC}_{\text{last}}$ and $C_{\max}$ represent the arithmetic average by cluster. AUC is calculated until the last observation ($\text{AUC}_{\text{last}}$). The number of clusters was selected by a visual inspection of the dendrogram in combination with the Calinski-Harabasz Index. For reference, the abbreviations are Poor Metabolizer (PM), Intermediate Metabolizer (IM), Extensive Metabolizer (EM), Rapid Metabolizer (RM), and Ultra-Rapid Metabolizer (UM). The cluster locations of the PMs are indicated in the dendrogram for reference.  | Dendrogram | Calinski-Harabasz (CH) Index by Cluster
---|---
|
|  Cluster | 1 | 2 | 3 | 4 | Total
---|---|---|---|---|---
# Obs. | 95 | 149 | 3 | 3 | 250
| Distribution of Metabolizer Status by Cluster |
PM | 6 (100%) | 0 (0%) | 0 (0%) | 0 (0%) | 6 (100%)
IM | 32 (61.54%) | 19 (36.54%) | 0 (0%) | 1 (1.92%) | 52 (100%)
EM | 36 (31.86%) | 75 (66.37%) | 2 (1.77%) | 0 (0%) | 113 (100%)
RM | 14 (22.22%) | 46 (73.02%) | 1 (1.59%) | 2 (3.17%) | 63 (100%)
UM | 7 (43.75%) | 9 (56.25%) | 0 (0%) | 0 (0%) | 16 (100%)
| Avg. PK Metrics by Cluster |
$\text{AUC}_{\text{last}}$ (ng*hr/mL)/Dose(mg) | 28.89 | 14.97 | 55.31 | 45.24 | 21.11
$C_{\max}$ (ng/mL)/Dose(mg) | 7.68 | 3.94 | 8.80 | 6.93 | 5.46

Figure S.4:  Detailed Case Study Clustering Results: CORT Distance. Summary statistics for $\text{AUC}_{\text{last}}$ and $C_{\max}$ represent the arithmetic average by cluster. AUC is calculated until the last observation ($\text{AUC}_{\text{last}}$). The number of clusters was selected by a visual inspection of the dendrogram in combination with the Calinski-Harabasz Index. For reference, the abbreviations are Poor Metabolizer (PM), Intermediate Metabolizer (IM), Extensive Metabolizer (EM), Rapid Metabolizer (RM), and Ultra-Rapid Metabolizer (UM). The cluster locations of the PMs are indicated in the dendrogram for reference.  | Dendrogram | Calinski-Harabasz (CH) Index by Cluster
---|---
|
|  Cluster | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | Total
---|---|---|---|---|---|---|---|---|---
# Obs. | 169 | 34 | 18 | 12 | 3 | 7 | 4 | 3 | 250
| Distribution of Metabolizer Status by Cluster |
PM | 1 (16.67%) | 0 (0%) | 4 (66.67%) | 1 (16.67%) | 0 (0%) | 0 (0%) | 0 (0%) | 0 (0%) | 6 (100%)
IM | 30 (57.69%) | 4 (7.69%) | 5 (9.62%) | 7 (13.46%) | 0 (0%) | 3 (5.77%) | 1 (1.92%) | 2 (3.85%) | 52 (100%)
EM | 85 (75.22%) | 13 (11.50%) | 6 (5.31%) | 1 (0.88%) | 2 (1.77%) | 4 (3.54%) | 1 (0.88%) | 1 (0.88%) | 113 (100%)
RM | 45 (71.43%) | 11 (17.46%) | 1 (1.59%) | 3 (4.76%) | 1 (1.59%) | 0 (0%) | 2 (3.17%) | 0 (0%) | 63 (100%)
UM | 8 (50.00%) | 6 (37.50%) | 2 (12.50%) | 0 (0%) | 0 (0%) | 0 (0%) | 0 (0%) | 0 (0%) | 16 (100%)
| Avg. PK Metrics by Cluster |
$\text{AUC}_{\text{last}}$ | 17.55 | 16.07 | 38.12 | 39.90 | 60.01 | 18.41 | 47.10 | 33.67 | 21.11
(ng*hr/mL)/Dose(mg) |  |  |  |  |  |  |  |  |
$C_{\max}$ | 5.10 | 3.54 | 10.22 | 8.22 | 10.22 | 2.79 | 7.48 | 6.48 | 5.46
(ng/mL)/Dose(mg) |  |  |  |  |  |  |  |  |