# Probabilistic Learning of Treatment Trees in Cancer

Tsung-Hung Yao
Department of Biostatistics
University of Michigan at Ann Arbor
Ann Arbor, MI 48109
<EMAIL_ADDRESS>
&Zhenke Wu
Department of Biostatistics
University of Michigan at Ann Arbor
Ann Arbor, MI 48109
<EMAIL_ADDRESS>
&Karthik Bharath
School of Mathematical Sciences
University of Nottingham
Nottingham, UK
<EMAIL_ADDRESS>
&Jinju Li
Department of Biostatistics
University of Michigan at Ann Arbor
Ann Arbor, MI 48109
<EMAIL_ADDRESS>
&Veerabhadran Baladandayuthapan
Department of Biostatistics
University of Michigan at Ann Arbor
Ann Arbor, MI 48109
<EMAIL_ADDRESS>
Corresponding author.

###### Abstract

Accurate identification of synergistic treatment combinations and their
underlying biological mechanisms is critical across many disease domains,
especially cancer. In translational oncology research, preclinical systems
such as patient-derived xenografts (PDX) have emerged as a unique study design
evaluating multiple treatments administered to samples from the same human
tumor implanted into genetically identical mice. In this paper, we propose a
novel Bayesian probabilistic tree-based framework for PDX data to investigate
the hierarchical relationships between treatments by inferring treatment
cluster trees, referred to as treatment trees (R${}_{\textrm{x}}$-tree). The
framework motivates a new metric of mechanistic similarity between two or more
treatments accounting for inherent uncertainty in tree estimation; treatments
with a high estimated similarity have potentially high mechanistic synergy.
Building upon Dirichlet Diffusion Trees, we derive a closed-form marginal
likelihood encoding the tree structure, which facilitates computationally
efficient posterior inference via a new two-stage algorithm. Simulation
studies demonstrate superior performance of the proposed method in recovering
the tree structure and treatment similarities. Our analyses of a recently
collated PDX dataset produce treatment similarity estimates that show a high
degree of concordance with known biological mechanisms across treatments in
five different cancers. More importantly, we uncover new and potentially
effective combination therapies that confer synergistic regulation of specific
downstream biological pathways for future clinical investigations. Our
accompanying code, data, and shiny application for visualization of results
are available at: https://github.com/bayesrx/RxTree.

_Keywords_ Approximate Bayesian Computation $\cdot$ Dirichlet Diffusion Trees
$\cdot$ Patient Derived Xenograft $\cdot$ Precision Medicine, Tree-Based
Clustering

##  1 Introduction

According to the World Health Organization, cancer is one of the leading
causes of death globally, with $\sim$10 million deaths in 2020 (Ferlay et al.,
2020). Despite multiple advances over the years, systematic efforts to predict
efficacy of cancer treatments have been stymied due to multiple factors,
including patient-specific heterogeneity and treatment resistance (Dagogo-Jack
and Shaw, 2018; Groisberg and Subbiah, 2021). Given that the evolution of
tumors relies on a limited number of biological mechanisms, there has been a
recent push towards combining multiple therapeutic agents, referred to as
“combination therapy” (Sawyers, 2013; Groisberg and Subbiah, 2021). This is
driven by the core hypothesis that combinations of drugs act in synergistic
manner, with each drug compensating for the drawbacks of other drugs. However,
despite higher response rates and efficacy in certain instances (Bayat
Mokhtari et al., 2017), combination therapy can lead to undesired drug-drug
interactions, lower efficacy, or severe side effects (Sun et al., 2016).
Consequently, it is highly desirable to advance the understanding of
underlying mechanisms that confer synergistic drug effects and identify
potential favorable drug-drug interaction mechanisms for further
investigations.

Given that not all possible drug combinations can be tested on patients in
actual clinical trials, cancer researchers rely on preclinical “model” systems
to guide the discovery of the most effective combination therapies (note,
models have a different contextual meaning here). In translational oncology,
preclinical models assess promising treatments and compounds, before they are
phased into human clinical trials. The traditional mainstay of such
preclinical models has been cell-lines, wherein cell cultures derived from
human tumors are grown in an in vitro controlled environment. However, it has
been argued that they do not accurately reflect the true behavior of the host
tumor and, in the process of adapting to in vitro growth, lose the original
properties of the host tumor, thus leading to limited clinical relevance and
successes (Tentler et al., 2012; Bhimani et al., 2020). To overcome these
challenges, there has been a push towards more clinically relevant model
systems that maintain a high degree of fidelity to human tumors. One such
preclinical model system is Patient-Derived Xenograft (PDX) wherein tumor
fragments obtained from cancer patients are directly transplanted into
genetically identical mice (Hidalgo et al., 2014; Lai et al., 2017). Compared
to traditional oncology models such as cell-lines (Yoshida, 2020), PDX models
maintain key cellular and molecular characteristics, and are thus more likely
to mimic human tumors and facilitate precision medicine. More importantly,
accumulating evidence suggests responses (e.g. drug sensitivity) to standard
therapeutic regimens in PDXs closely correlate with patient clinical data,
making PDX an effective and predictive experimental model across multiple
cancers (Topp et al., 2014; Nunes et al., 2015).

Figure 1: PDX experimental design and tree-based representation. Panel A: an
illustrative PDX dataset with five treatments (row) and eight patients
(column). Mice in a given column are implanted with tumors from the same
patient and receive different treatments (across rows). The level of tumor
responses are shown along a color gradient. Panel B: a tree structure that
clusters the treatments and quantifies the similarity among mechanisms. Two
treatments (1 and 4) are assumed to have different but known biological
mechanisms (in different colors); the rest three treatments (2,3, and 5) have
unknown mechanisms (in gray). The tree suggests two treatment groups are
present ($\\{1,2\\}$ and $\\{3,4,5\\}$) that may correspond to two different
known mechanisms. The horizontal position of “$\triangle$” represents the
divergence time (defined in Section 2.1) and the mechanism similarity for
treatments $\\{3,4,5\\}$. In a real data analysis, the tree (topology and
divergence times) is unknown and is to be inferred from PDX data.

##### PDX experimental design and key scientific questions.

Briefly, in a standard PDX preclinical experiment, a set of common treatments
are tested, and each treatment is given to multiple mice with tumors implanted
from the same (matched) patient (see conceptual schema in Figure 1(A)).
Treatment responses (e.g. tumor size) are then evaluated, resulting in a data
matrix (treatments $\times$ patients) as depicted in the heatmap in Figure
1(A). The PDX experimental protocol, often referred to as “co-clinical
trials,” mirrors a real human clinical trial using mouse “avatars” (Clohessy
and Pandolfi, 2015). Thus the protocol serves as a scalable platform to: (a)
identify underlying plausible biological mechanisms responsible for tumor
growth and resistance, and (b) evaluate the effectiveness of drug combinations
based on mechanistic understanding (Lunardi and Pandolfi, 2015). In this
context, the (biological) mechanism refers to the specific mechanism of action
of a treatment, which usually represents a specific target, such as an enzyme
or a receptor (Grant et al., 2010). From the perspective of treatment
responses as data, responses are the consequences of the downstream biological
pathways from the corresponding interaction between a treatment and the
target/mechanism.

Ideally, treatments with the same target/mechanism should induce similar
responses and engender mechanism-related clustering among treatments.
Evidently then, a sensible clustering of treatments would not only partition
treatments into clusters but also explicate how the clusters relate to one
another; in other words, a hierarchy among treatment clusters is more likely
to uncover plausible mechanisms for combinations of treatments with “similar”
responses when compared to “flat” clusters (e.g., $k$-means clustering). Such
response-based identification of potential synergistic effects from
combinations of treatments will augment understanding from known mechanistic
synergy. In our application, using tree-based clustering, we assume known
entities at the leaves, i.e., the different treatments. The treatments are
assumed to act upon potentially distinct biological pathways, resulting in
different levels of responses across the treated mice. In this paper, we use
PDX response data on the leaves to infer a hierarchy over treatments that may
empirically characterize the similarity in the targeted mechanistic pathways.
The primary statistical goals are to (i) define and estimate a general metric
measuring the similarity within any subset comprising two or more treatments,
and (ii) facilitate (i) by conceptualizing and inferring an unknown hierarchy
among treatments.

##### Tree-based representations for PDX data.

To this end, we consider a tree-based construct to explore the hierarchical
relationships between treatments, referred to as treatment tree
(R${}_{\textrm{x}}$-tree, in short). We view such a tree structure as a
representation of clustering of treatments based on mechanisms that confer
synergistic effects, wherein similarities between mechanisms are captured
through branch lengths. Hierarchy among treatments can be interpreted through
branch lengths (from the root) that are potentially reflective of different
cancer processes; this would then help identify common mechanisms and point
towards treatment combinations disrupting oncological processes if
administered simultaneously.

We will focus on rooted trees. The principal ingredients of a rooted tree
comprise a root node, terminal nodes (or, leaves), internal nodes and branch
lengths. In the context of the R${}_{\textrm{x}}$-tree for PDX data, the
leaves are observed treatment responses, whereas internal nodes and branch
lengths are unobserved. Internal nodes are clusters of treatments, and lengths
of branches between nodes are indicative of strengths of mechanism
similarities. The root is a single cluster consisting of all treatments. This
leads to the following interpretation: at the root all treatments share a
common target or mechanism; length of path from the root to the internal node
(sum of branch lengths) at which two treatments split into different clusters
measures mechanism similarity between the two treatments. Thus treatments that
stay clustered “longer” have higher mechanism similarities.

Throughout, we will use ‘tree’ when describing methodology for an abstract
tree (acyclic graphs with distinguished root node) and ‘treatment tree’ or
‘R${}_{\textrm{x}}$-tree’ when referring to the latent tree within the
application context.

##### An illustrative example.

A conceptual R${}_{\textrm{x}}$-tree and its interpretation is illustrated in
Figure 1 where five treatments (1 to 5) are applied on eight patients’ PDXs
(Figure 1(A)) with the corresponding (unknown true) R${}_{\textrm{x}}$-tree
(Figure 1(B)) based on the PDX data. Assume two treatment groups based on
different mechanisms – treatments $\\{1,2\\}$ and treatments $\\{3,4,5\\}$;
further, suppose treatment $4$ is approved by the Food and Drug Administration
(FDA). The heatmap in Panel (A) visualizes the distinct levels of response
profiles to the five treatments so that treatments closer in the tree are more
likely to have similar levels of responses. The R${}_{\textrm{x}}$-tree
captures the mechanism similarity by arranging treatments $\\{1,2\\}$ and
$\\{3,4,5\\}$ to stay in their respective subtrees longer and to separate the
two sets of treatment early in the tree. Based on the R${}_{\textrm{x}}$-tree,
treatments $\\{3,5\\}$ share high mechanism similarity values with treatment
$4$; treatment $5$ is the closest to the treatment $4$, suggesting the most
similar synergistic mechanism among all the evaluated treatments $1$ to $5$.

##### Existing methods and modeling background.

The Pearson correlation is a popular choice to assess mechanism similarity
between treatments (Krumbach et al., 2011), but is inappropriate to examine
multi-way similarity. A tree-structured approach based on a (binary)
dendrogram obtained from hierarchical clustering of cell-line data using the
cophenetic distance (Sokal and Rohlf, 1962) was adopted in Narayan et al.
(2020); their approach, however, failed to account for uncertainty in the
dendrogram, which is highly sensitive to measurement error in the response
variables as well distance metrics (we show this via simulations and in real
data analyses). In this paper, we consider a model for PDX data parameterized
by a tree-structured object representing the R${}_{\textrm{x}}$-tree. The
model is derived from the Dirichlet diffusion tree (DDT) (Neal, 2003)
generative model for (hierarchically) clustered data. The DDT engenders a data
likelihood and a prior distribution on the tree parameter with support in the
space of rooted binary trees. We can then use the posterior distribution to
quantify uncertainty about the latent R${}_{\textrm{x}}$-tree.

##### Summary of novel contributions and organization of the article.

Our approach based on the DDT model for PDX data results in three main novel
contributions:

  1. (a)

Derivation of a closed-form likelihood that encodes the tree structure. The
DDT specification results in a joint distribution on PDX data, treatment tree
parameters and other model parameters. By marginalizing over unobserved data
that correspond to internal nodes of the tree, we obtain a new multivariate
Gaussian likelihood with a special tree-structured covariance matrix, which
completely characterizes the treatment tree (Proposition 1 and Lemma 1).

  2. (b)

Efficient two-stage algorithm for posterior sampling. Motivated by the form of
marginal data likelihood in (a), we decouple the Euclidean and tree parameters
and propose a two-stage algorithm that combines an approximate Bayesian
computation (ABC) procedure (for Euclidean parameters) with a Metropolis-
Hasting (MH) step (for tree parameters). We demonstrate via multiple
simulation studies the superiority of our hybrid approach over approaches
based on classical single-stage MH algorithms (Sections 4.2 and 4.1).

  3. (c)

Corroborating existing, and uncovering new, synergistic combination therapies.
We define and infer a new similarity measure that accounts for inherent
uncertainty in estimating a latent hierarchy among treatments. As a result,
the maximum a posteriori R${}_{\textrm{x}}$-tree and the related mechanism
similarity show high concordance with known existing biological mechanisms for
monotherapies and uncover new and potentially useful combination therapies
(Sections 5.3 and 5.4).

Of particular note is contribution (c), where we leverage a recently collated
PDX dataset from the Novartis Institutes for BioMedical Research - PDX
Encyclopedia [NIBR-PDXE, (Gao et al., 2015)] that interrogated multiple
targeted therapies across five different cancers. Our pan-cancer analyses of
the NIBR-PDXE dataset show a high degree of concordance with known existing
biological mechanisms across different cancers; for example, a high
mechanistic similarity is suggested between two agents currently in clinical
trials: CGM097 and HDM201 in breast cancer and colorectal cancer, known to
target the same gene MDM2 (Konopleva et al., 2020). In addition, our model
uncovers new and potentially effective combination therapies. For example,
exploiting knowledge of the combination therapy of a class of agents targeting
the PI3K-MAPK-CDK pathway axes – PI3K-CDK for breast cancer, PI3K-ERBB3 for
colorectal cancer and BRAF-PI3K for melanoma – confers possible synergistic
regulation for prioritization in future clinical studies.

The rest of the paper is organized as follows: we first review our
probabilistic formulation for PDX data based on the DDT model and present the
marginal data likelihood and computational implications in Section 2. In
Section 3, we derive the posterior inference algorithm based on a two-stage
algorithm. In Section 4, we conduct two sets of simulations to evaluate the
operating characteristics of the model and algorithm. A detailed analysis of
the NIBR-PDXE dataset, results, biological interpretations and implications
are summarized in Section 5. The paper concludes by discussing implications of
the findings, limitations, and future directions.

##  2 Modeling R${}_{\textrm{x}}$-tree via Dirichlet Diffusion Trees

Given a PDX experiment with $I$ correlated treatments and $J$ independent
patients, we focus on the setting with $1\times 1\times 1$ design (one animal
per PDX model per treatment) with no replicate response for each treatment and
patient. A PDX experiment produces an observed data matrix
$\mathbf{X}_{I\times J}=[\bm{X}_{1},\ldots,\bm{X}_{I}]^{\sf T}$ where
$\bm{X}_{i}=[X_{i1},\ldots,X_{iJ}]^{\sf T}$ is data under treatment $i$ across
$J$ patients; let the observed response column for each patient be
$\bm{X}_{\cdot,j}=[x_{1j},\ldots,x_{Ij}]^{\sf
T}\in\mathbb{R}^{I},j=1,\ldots,J$.

In this paper, the observed treatment responses are continuous and we model
the responses through a generative model that results in a Gaussian likelihood
with a structured covariance:

$\displaystyle\bm{X}_{\cdot,j}|\bm{\Sigma}^{\mathcal{T}},\sigma^{2}\stackrel{{\scriptstyle
iid}}{{\sim}}{\sf N}_{I}(0,{\bm{\Sigma}^{\mathcal{T}}}),\quad j=1,\ldots,J,$
(1)

where the $\bm{\Sigma}^{\mathcal{T}}$ is a tree-structured covariance matrix
that encodes the tree ${\mathcal{T}}$. In particular,
$\bm{\Sigma}^{{\mathcal{T}}}=\left\\{\bm{\Sigma}^{\mathcal{T}}_{i,i^{\prime}},i,i^{\prime}=1,\ldots,I\right\\}$
encodes the tree ${\mathcal{T}}$ through two constraints (Lapointe and
Legendre, 1991; McCullagh, 2006):

$\displaystyle\bm{\Sigma}^{{\mathcal{T}}}_{i^{\prime},i}=\bm{\Sigma}^{{\mathcal{T}}}_{i,i^{\prime}}\geq
0;\thickspace$
$\displaystyle\bm{\Sigma}^{{\mathcal{T}}}_{i,i}\geq\bm{\Sigma}^{{\mathcal{T}}}_{i,i^{\prime}},$
(2)
$\displaystyle\bm{\Sigma}^{{\mathcal{T}}}_{i,i^{\prime}}\geq\min\\{\bm{\Sigma}^{\mathcal{T}}_{i,i{{}^{\prime\prime}}}\
,\ \bm{\Sigma}^{\mathcal{T}}_{i^{\prime},i{{}^{\prime\prime}}}\\}$
$\displaystyle\text{ for all }i\neq i^{\prime}\neq i{{}^{\prime\prime}}.$ (3)

Each element $\bm{\Sigma}^{{\mathcal{T}}}_{i,i}$ is the covariance between
treatments $i$ and $i^{\prime}$ and measures their similarity. The inequality
(2) imposes the symmetry of covariance matrix and ensures the divergence of
all leaves. The tree structure is characterized by the ultrametric inequality
(3) that ensures $\bm{\Sigma}^{{\mathcal{T}}}$ bijectively maps to a tree
${\mathcal{T}}$; for more details on the relationship between the covariance
$\bm{\Sigma}^{{\mathcal{T}}}$ and the tree ${\mathcal{T}}$ see McCullagh
(2006) and Bravo et al. (2009). Of note, mean parameterized models (e.g. mixed
effects models) are inappropriate for uncovering the tree parameter under the
given data structure since the latent tree is completely encoded in covariance
matrix $\bm{\Sigma}^{\mathcal{T}}$.

A Bayesian formulation requires an explicit prior distribution on
$\bm{\Sigma}^{{\mathcal{T}}}$ which satisfies constrains (2) and (3); this
requirement is far from straightforward since the set of tree-structured
matrices is complicated (e.g., it is not a manifold (McCullagh, 2006)). We
instead consider the Dirichlet Diffusion tree (DDT) model (Neal, 2003) for
hierarchically clustered data which provides two useful ingredients:

  1. 1.

a prior is implicitly specified on the latent treatment tree, comprising the
root, internal nodes, leaves, and branch lengths;

  2. 2.

upon integrating out the internal nodes, a tractable Gaussian likelihood on
PDX data with tree-structured covariance is specified.

We first provide a brief description of the DDT model proposed by Neal (2003)
and its joint density on data and tree (Section 2.1). Subsequently, we derive
an expression for the likelihood and demonstrate how it can be profitably
employed to develop a generative model for PDX data and carry out
R${}_{\textrm{x}}$-tree estimation (Section 2.2 and 2.3).

###  2.1 The Generative Process of DDT

The DDT prescribes a fragmentary, top-down mechanism to generate a binary tree
(acyclic graph with a preferred node or vertex referred to as the root),
starting from a root, containing $J$-dimensional observed responses
$\bm{X}_{i}$ at $I$ leaves/terminal nodes; each node in the tree has either 0
or 2 children excepting the root which has a solitary child. This prescription
manipulates dynamics of a system of $I$ independent Brownian motions
$B_{1},\ldots,B_{I}$ on $\mathbb{R}^{J}$ in a common time interval
$t\in[0,1]$. As shown in Figure 2(A), all Brownian motions $B_{i}(t)$ start at
the same point at time $t=0$, location of which is the root
$\bm{0}\in\mathbb{R}^{J}$, and diverge at time points in $[0,1]$ and locations
in $\mathbb{R}^{J}$ before stopping at the time $t=1$ at locations
$\bm{X}_{i}$. The Brownian trajectories and their divergences engender the
tree structure as shown in Figure 2(A).

Specifics on when and how the Brownian motions diverge are as follows: the
first Brownian motion $B_{1}(t)$ starts at $t=0$ and generates $\bm{X}_{1}$ at
$t=1$; a second independent Brownian motion $B_{2}(t)$ starts at the same
point at $t=0$, branches out from the first Brownian motion at some time $t$,
after which it generates $\bm{X}_{2}$ at time $1$. The probability of
divergence in a small interval $[t,t+dt]$ is given by a _divergence function_
$t\mapsto a(t)$, assumed as in Neal (2003) to be of the form
$a(t)=c(1-t)^{-1}$ for some divergence parameter $c>0$. Inductively then, the
vector of observed responses to treatment $i$, $\bm{X}_{i}$, is generated by
$B_{i}(t)$, which follows the path of previous ones. If at time $t$,
$B_{i}(t)$ has not diverged and meets the previous divergent point, it will
follow one of the existing path with the probability proportional to the
number of data points that have previously traversed along each path.
Eventually, given $B_{i}(t)$ has not diverged at time $t$, it will do so in
$[t,t+dt]$ with probability $a(t)dt/m$, where $m$ is the number of data points
that have previously traversed the current path.

From the illustration in panel (A) of Figure 2, we note that $B_{3}$ diverges
from the $B_{1}$ and $B_{2}$ at time $t_{1}$ at location $\bm{X}_{1}^{\prime}$
and at $t=1$ is at location $\bm{X}_{3}$, which is the $J$-dimensional
response vector for treatment 3; this creates a solitary branch of length
$t_{1}$ from the root and an unobserved internal node at location
$\bm{X}_{1}^{\prime}$. Continuing, given three Brownian motions $B_{1}$,
$B_{2}$ and $B_{3}$, $B_{4}$ does not diverge before $t_{1}$ and meet the
previous divergent point $t_{1}$. $B_{4}$ chooses to follow the path of
$B_{3}$ with probability $1/3$ at $t_{1}$ and finally diverges from $B_{3}$ at
time $t_{2}>t_{1}$ at location $\bm{X}_{2}^{\prime}$; this results in
observation $\bm{X}_{4}$ for treatment 4 and an unobserved internal node at
$\bm{X}_{2}^{\prime}$, and so on. As a consequence, the binary tree that
arises from the DDT comprises of:

  1. (i)

an unobserved root at the origin in $\mathbb{R}^{J}$ at time $t=0$;

  2. (ii)

observed data $\mathbf{X}=[\bm{X}_{1},\ldots,\bm{X}_{I}]^{\sf
T}\in\mathbb{R}^{I\times J}$ situated at the leaves of the tree;

  3. (iii)

unobserved internal nodes
$\mathbf{X}^{I}=[\bm{X}_{1}^{\prime},\ldots,\bm{X}^{\prime}_{I-1}]^{\sf
T}\in\mathbb{R}^{(I-1)\times J}$;

  4. (iv)

unobserved times $\bm{t}=(t_{1},\ldots,t_{I-1})^{\sf T}\in[0,1]^{I-1}$ that
characterize lengths of branches;

  5. (v)

unobserved topology $\mathcal{T}$ that links (i)-(iv) into a tree structure,
determined by the number of data points $\bm{X}_{i}$ that have traversed
through each segment or branch.

Conceptually, observed data at the leaves $\bm{X}_{1},\ldots,\bm{X}_{I}$
collectively form the observed PDX responses generated through a process
involving a few parameters: tree-related parameters $(\mathcal{T},\bm{t})$ and
the locations of internal nodes $\bm{X}^{\prime}_{i}$. The tree $\mathcal{T}$
clusters $I$ treatments as a hierarchy of $(I-1)$ levels (excluding the last
level containing leaves). At level $0<d\leq I-1$ of the hierarchy,
characterized by the pair $(\bm{X}^{\prime}_{d},t_{d})$, the $I$ treatments
are clustered into $d+1$ groups; a measure of similarity (or dissimilarity)
between treatment clusters at levels $d$ and $d+1$ is given by the branch
length $t_{d+1}-t_{d}$.

Figure 2: (A) A binary tree with $I=5$ leaves underlying the diffusion
dynamics. The observed response vector $\bm{X}_{i},i=1,\ldots,I$ is generated
by the Brownian motion up to $t=1$. The unobserved response vector
$\bm{X}_{d}^{\prime},d=1,\ldots,(I-1)$ at the divergence is generated by the
Brownian motion at time $t_{d}$. (B) A tree-structured matrix
$\bm{\Sigma}^{\mathcal{T}}$ that encapsulates the tree ${\mathcal{T}}$. See
the Proposition 1 for the definition of $\bm{\Sigma}^{\mathcal{T}}$.

We now give a brief description of how the joint density of
$(\mathbf{X},\mathbf{X}^{I},\bm{t},\mathcal{T})$ can be derived; for more
details we direct the reader to Neal (2003) and Knowles and Ghahramani (2015).
For a fixed $c>0$ that governs the divergence function $a(t)=c(1-t)^{-1}$,
probabilities associated with the independent Brownian motions
$B_{1},\ldots,B_{I}$ induce a joint (Lebesgue) density on the generated tree.
Note that the binary tree arising from the DDT is encoded by the triples
$\\{(t_{d},\bm{X}^{\prime}_{d},\bm{X}_{i}),d=1,\ldots,I-1;i=1,\ldots,I\\}$. An
internal node at $\bm{X}^{\prime}_{d}$ contains $l_{d}$ and $r_{d}$ leaves
below to its left and right with $m_{d}=l_{d}+r_{d}$. If each of the Brownian
motions is scaled by $\sigma^{2}>0$, then given $\mathcal{T}$ and a branch
with endpoints $(t_{u},\bm{X}^{\prime}_{u})$ and $(t_{v},\bm{X}^{\prime}_{v})$
with $0<t_{u}<t_{v}<1$, from properties of a Brownian motion we see that
$\bm{X}^{\prime}_{v}\sim
N_{J}(\bm{X}^{\prime}_{u},\sigma^{2}(t_{v}-t_{u})\bm{I}_{J})$, and the
(Lebesgue) density of ${\mathcal{T}}$ can be expressed as the product of
contributions from its branches. Then the joint density of all nodes, times
and the tree topology is given by

$\displaystyle
P(\mathbf{X},\mathbf{X}^{I},\bm{t},\mathcal{T}|c,\sigma^{2})=\Pi_{[u,v]\in\mathcal{S}(\mathcal{T})}\frac{(l_{v}-1)!(r_{v}-1)!}{(l_{v}+r_{v}-1)!}c(1-t_{v})^{cJ_{l_{v},r_{v}}-1}N_{J}(\bm{X}^{\prime}_{u},\sigma^{2}(t_{v}-t_{u})\bm{I}_{J})$
(4)

where $\mathcal{S}(\mathcal{T})$ is the collection of branches and
$\mathbf{X}^{I}_{(I-1)\times
J}=[\bm{X}^{\prime}_{1},\ldots,\bm{X}^{\prime}_{(I-1)}]^{\sf T}$ are
unobserved locations of the internal nodes. On each branch $[u,v]$, the first
term $\frac{(l_{v}-1)!(r_{v}-1)!}{(l_{v}+r_{v}-1)!}$ represents the chance the
branch containing $l_{v}$ and $r_{v}$ leaves to its left and right
respectively; $c(1-t_{v})^{cJ_{l_{v},r_{v}}-1}$ represents the probability of
diverging at $t_{v}$ with $l_{v}$ and $r_{v}$ leaves, where
$J_{l+r}=H_{l_{v}+r_{v}-1}-H_{l_{v}-1}-H_{r_{v}-1}$ with
$H_{n}=\sum_{i=1}^{n}1/i$ is the $n$th harmonic number.

The joint density is hence parameterized by $(c,\sigma^{2})$, where $c$ plays
a crucial role in determining the topology $\mathcal{T}$: through the
divergence function $a(t)$, it determines the propensity of the Brownian
motion to diverge from its predecessors; consequently, a small $c$ engenders
later divergence and a higher degree of similarity among treatments in PDX.
The latent tree has two components: (i) topology $\mathcal{T}$ and (ii) vector
of divergence times $\bm{t}$ determining branch lengths. We refer to
$(c,\sigma^{2})$ as the _Euclidean parameters_ and $(\mathcal{T},\bm{t})$ as
_tree parameters_.

###  2.2 Prior on tree and closed-form likelihood

The joint density in (4) factors into a prior
$P(\bm{t},\mathcal{T}|c,\sigma^{2})$ on the tree parameter through
$(\mathcal{T},\bm{t})$ and a density
$P(\mathbf{X},\mathbf{X}^{I}|\bm{t},\mathcal{T},c,\sigma^{2})$ that is a
product of $J$-dimensional Gaussians on the internal nodes and leaves. The
prior distribution on the latent tree is thus implicitly defined through the
Brownian dynamics and is parameterized by $(\mathcal{T},\bm{t})$ with
hyperparameters $(c,\sigma^{2})$. In (4) the product is over the set of
branches $\mathcal{S}(\mathcal{T})$, and the contribution to the prior
$P(\mathcal{T},\bm{t}|c,\sigma^{2})$ from each branch $[u,v]$ is
$\frac{(l_{v}-1)!(r_{v}-1)!}{(l_{v}+r_{v}-1)!}c(1-t_{v})^{cJ_{l_{v},r_{v}}-1}$,
which is free of $\sigma^{2}$; on the other hand, the contribution to
$P(\mathbf{X},\mathbf{X}^{I}|\bm{t},\mathcal{T},c,\sigma^{2})$ from $[u,v]$ is
the $J$-dimensional
$N_{J}(\bm{X}^{\prime}_{u},\sigma^{2}(t_{v}-t_{u})\bm{I}_{J})$, which is
independent of $c$. The likelihood function based on the observed $\mathbf{X}$
is thus obtained by integrating out the unobserved internal nodes
$\mathbf{X}^{I}$ from
$P(\mathbf{X},\mathbf{X}^{I}|\bm{t},\mathcal{T},\sigma^{2})$. Accordingly, our
first contribution is to derive a closed-form likelihood function for
efficient posterior computations; to our knowledge, this task is currently
achieved only through sampling-based or variational methods (Neal, 2003;
Knowles and Ghahramani, 2015).

Denote as ${\sf MN}_{I\times J}(M,U,V)$ the matrix normal distribution of an
$I\times J$ random matrix with mean matrix $M$, row covariance $U$, and column
covariance $V$, and let $\bm{I}_{k}$ denote the $k\times k$ identity matrix.
Evidently, $\mathbf{X}$ follows a matrix normal distribution since Gaussian
laws of the Brownian motions imply that
$[\mathbf{X},\mathbf{X}^{I}]=[\bm{X}_{1},\ldots,\bm{X}_{I},\bm{X}^{\prime}_{1},\ldots,\bm{X}^{\prime}_{(I-1)}]^{\sf
T}$ follow a matrix normal distribution.

######  Proposition 1.

Under the assumption that the root is located at the origin in
$\mathbb{R}^{J}$, the data likelihood
$\mathbf{X}|\sigma^{2},\mathcal{T},\bm{t}\sim{\sf MN}_{I\times
J}(\bm{0},\sigma^{2}{\bm{\Sigma}}^{{\mathcal{T}}},\bm{I}_{J})$, where
${\bm{\Sigma}^{\mathcal{T}}}=\Big{(}\bm{\Sigma}^{\mathcal{T}}_{i,i^{\prime}}\Big{)}$
is an $I\times I$ tree-structured covariance matrix satisfying (2) and (3)
with $\bm{\Sigma}^{\mathcal{T}}_{i,i}=1$ and
$\bm{\Sigma}^{\mathcal{T}}_{i,i^{\prime}}=t_{d}$, for $i\neq i^{\prime}$ where
$i,i^{\prime}=1,\ldots,I$ and $d=1,\ldots,I-1$.

Proposition 1 asserts that use of the DDT model leads to a centered Gaussian
likelihood on PDX data $\bm{X}$ with a tree-structured covariance matrix.
Proposition 1 also implies that each patient independently follows the normal
distribution of (1) with an additional scale parameter $(\sigma^{2})$ from the
Brownian motion:

$\displaystyle\bm{X}_{\cdot,j}|\bm{\Sigma}^{\mathcal{T}},\sigma^{2}\stackrel{{\scriptstyle
iid}}{{\sim}}{\sf N}_{I}(0,{\sigma^{2}\bm{\Sigma}^{\mathcal{T}}}),\quad
j=1,\ldots,J.$ (5)

By setting $\bm{\Sigma}^{\mathcal{T}}_{i,i^{\prime}}=t_{i,i^{\prime}}$ as the
divergence time of $i$ and $i^{\prime}$, $\bm{\Sigma}^{\mathcal{T}}$ satisfies
(2) and (3) and encodes the tree ${\mathcal{T}}$. For example, consider a
three-leaf tree with
$\bm{\Sigma}^{\mathcal{T}}_{i,i^{\prime}}=t_{i,i^{\prime}}$, inequality (3)
implies that for the three leaves, say, $i,i^{\prime}$ and $i^{\prime\prime}$,
one of the following conditions must hold: (i)
$t_{i^{\prime},i^{\prime\prime}}\geq t_{i,i^{\prime}}=t_{i,i^{\prime\prime}}$;
(ii) $t_{i,i^{\prime\prime}}\geq
t_{i,i^{\prime}}=t_{i^{\prime},i^{\prime\prime}}$; (iii) $t_{i,i^{\prime}}\geq
t_{i,i^{\prime\prime}}=t_{i^{\prime},i^{\prime\prime}}$. We then obtain a tree
containing 1) a subtree of two leaves with a higher similarity and 2) a
singleton clade with a lower similarity between the singleton leaf and the two
leaves in the first subtree. In particular, if
$t_{i^{\prime},i^{\prime\prime}}\geq t_{i,i^{\prime}}=t_{i,i^{\prime\prime}}$
holds, the three-leaf tree has leaf $i$ diverging earlier before the subtree
of $(i^{\prime},i^{\prime\prime})$.

###  2.3 Decoupling Tree and Euclidean Parameters for Efficient Sampling.

In the full joint density in (4) the Euclidean and tree parameters are
confounded across row and column dimensions of $\mathbf{X}$, and this may
result in slow mixing of chains using traditional MCMC algorithms (Turner et
al., 2013). State-of-the-art posterior inference on
$(c,\sigma^{2},\mathcal{T},\bm{t})$ can be broadly classified into sampling-
based approaches (e.g., Knowles and Ghahramani, 2015) and deterministic
approaches based on variational message passing (e.g., Knowles et al., 2011,
VMP). Variational algorithms can introduce approximation errors to the joint
posterior via factorization assumptions (e.g.,mean-field) and choice of
algorithm is typically determined by the speed-accuracy trade-off tailored for
particular applications. On the other hand, in classical MCMC-based algorithms
for DDT we observed slow convergence in the sampling chains for $c$ and
$\sigma^{2}$ with high autocorrelations for the corresponding chains, owing to
possibly the high mutual dependence between $c$ in the divergence function and
the tree topology ${\mathcal{T}}$ , resulting in slow local movements in the
joint parameter space of model and tree parameters (Simulation II in Section
4.2).

Notwithstanding absence of the parameter $c$ in the Gaussian likelihood, the
dependence, and information about, $c$ is implicit: the distribution of
divergence times $\bm{t}$ that populate $\bm{\Sigma}^{\mathcal{T}}$ are
completely determined by the divergence function $t\mapsto c(1-t)^{-1}$. In
other words, $c$ can indeed be estimated from treatment responses
$\\{\bm{X}_{\cdot,j}\\}$ using the likelihood. From a sampling perspective,
however, form of the likelihood obtained by integrating out the internal nodes
$\mathbf{X}^{I}$, suggests an efficient two-stage sampling strategy that
resembles the classical collapsed sampling (Liu, 1994) strategy in MCMC
literature: first draw posterior samples of $(c,\sigma^{2})$ and then proceed
to draw posterior samples of $(\mathcal{T},\bm{t})$ conditioned on each sample
of $(c,\sigma^{2})$.

##  3 R${}_{\textrm{x}}$-tree Estimation and Posterior Inference

In line with the preceding discussion, we consider a two-stage sampler for
Euclidean and tree parameters. While in principle MCMC techniques could be
used in both stages, we propose to use a hybrid ABC-MH algorithm.
Specifically, we use an approximate Bayesian computation (ABC) scheme to draw
weighted samples of $(c,\sigma^{2})$ in the first stage followed by a
Metropolis-Hastings (MH) step that samples $({\mathcal{T}},\bm{t})$ given ABC
samples of $(c,\sigma^{2})$ in the second stage. Motivation for using ABC in
the first stage stems from: (i) availability of informative statistics; (ii)
generation of better quality samples of the tree (compared to a single-stage
MH); and (iii) better computational efficiency. We refer to Section 4.2 for
more details.

###  3.1 Hybrid ABC-MH Algorithm

ABC is a family of inference techniques that are designed to estimate the
posterior density $\textrm{pr}(\theta|\mathcal{D})$ of parameters $\theta$
given data $\mathcal{D}$ when the corresponding likelihood
$\textrm{pr}(\mathcal{D}|\theta)$ is intractable but fairly simple to sample
from. Summarily, ABC approximates $\textrm{pr}(\theta|\mathcal{D})$ by
$\textrm{pr}(\theta|\bm{S}_{obs})$ where $\bm{S}_{obs}$ is a $d$-dimensional
summary statistic that ideally captures most information about $\theta$. In
the special case where $\bm{S}_{obs}$ is a vector of sufficient statistics, it
is well known that
$\textrm{pr}(\theta\mid\mathcal{D})=\textrm{pr}(\theta\mid\bm{S}_{obs})$. To
generate a sample from the partial posterior distribution
$\textrm{pr}(\theta\mid\bm{S}_{obs})$, ABC with rejection sampling proceeds
by: (i) simulating $N^{\textrm{syn}}$ values
$\theta_{l},l=1,\ldots,N^{\textrm{syn}}$ from the prior distribution
$\textrm{pr}(\theta)$; (ii) simulating datasets $D_{l}$ from
$\textrm{pr}(\mathcal{D}|\theta_{l})$; (iii) computing summary statistics
$\bm{S}_{l},l=1,\ldots,N^{\textrm{syn}}$ from $\mathcal{D}_{l}$; (iv)
retaining a subset of $\\{\theta_{l_{s}},s=1\ldots,k\\}$ of size
$k<N^{\textrm{syn}}$ that corresponds to ‘small’
$\|\bm{S}_{l_{s}}-\bm{S}_{obs}\|$ values based on some threshold. Given pairs
$\\{(\theta_{l_{s}},\bm{S}_{l_{s}})\\}$, the task of estimating the partial
posterior translates to a problem of conditional density estimation, e.g.,
based on Nadaraya-Waston type estimators and local regression adjustment
variants to correct for the fact that $\bm{S}_{l_{s}}$ may not be exactly
$\bm{S}_{obs}$; see Sisson et al. (2018) for a comprehensive review. To
implement ABC, the choice of summary statistics is central.

We detail the specialization of ABC to the marginal posterior distributions of
$c$ and $\sigma^{2}$ in Section 3.1.1. Given any pair of $(c,\sigma^{2})$, we
can sample trees from a density function up to an unknown normalizing constant
based on an existing MH algorithm (Knowles and Ghahramani, 2015). Our proposal
is to condition on the posterior median of $(c,\sigma^{2})$ of ABC-weighted
samples from the first stage, when sampling the trees in the second stage;
clearly, other choices are also available. This strategy produced comparable
MAP trees and inference of other tree-derived results relative to tree samples
based on full ABC samples of $c$ and $\sigma^{2}$.

Pseudo code for the two-stage algorithm is presented in the Supplementary
Material Algorithm S1. We briefly describe below its key components.

####  3.1.1 Stage 1: Sampling Euclidean Parameters $(c,\sigma^{2})$ using ABC

Accuracy and efficiency of the ABC procedure is linked to two competing
desiderata on the summary statistics: (i) informative, or ideally sufficient;
(ii) low-dimensional.

Summary statistic for $\sigma^{2}$. From the closed-form likelihood in
Equation (5), a sufficient statistic of $\sigma^{2}\bm{\Sigma}^{\mathcal{T}}$
is easily available, using which we construct a summary statistics for
$\sigma^{2}$.

######  Lemma 1.

With $\mathbf{X}$ as the observed data, the statistic
$\bm{T}:=\sum_{j}\bm{X}_{\cdot,j}\bm{X}_{\cdot,j}^{\sf T}$ is sufficient for
$\sigma^{2}{\bm{\Sigma}^{\mathcal{T}}}$ and follows a Wishart distribution
$W_{I}(J,\sigma^{2}\bm{\Sigma}^{{\mathcal{T}}})$, where
$\bm{X}_{\cdot,j}=[x_{1j},\ldots,x_{Ij}]\in\mathbb{R}^{I}$. Then with
$S^{(\sigma^{2})}:=\frac{tr(\bm{T})}{IJ}$ we have
$E[S^{(\sigma^{2})}]=\sigma^{2}$ and
$Var[S^{(\sigma^{2})}]=\frac{2\sigma^{4}tr(({\bm{\Sigma}^{\mathcal{T}}})^{2})}{I^{2}J}$.

Due to the normality of $\mathbf{X}$ in (5), and the Factorization theorem
(Casella and Berger, 2001), we see that $\bm{T}$ is complete and sufficient
for $\sigma^{2}\bm{\Sigma}^{\mathcal{T}}$ and $\bm{T}\sim
W_{I}(J,\sigma^{2}{\bm{\Sigma}^{\mathcal{T}}})$. Well-known results about the
trace and determinant of $\mathbf{X}$ (see for e.g. Mathai (1980)) provide the
stated results on the mean and variance of $tr(\bm{T})$. Owing to its
unbiasedness, we choose $S^{(\sigma^{2})}={tr(\bm{T})}/{IJ}$ as the summary
statistic for $\sigma^{2}$ and examine its performance through simulations in
Section 4; other choices are assessed in the Supplementary Material Section
S4.1.

Summary statistic for $c$. Based on the matrix normal distribution of
Proposition 1, the divergence parameter $c$ does not appear in the observed
data likelihood. Any statistic based on the entire observed data set
$\mathbf{X}$ is sufficient, but not necessarily informative about $c$. In DDT,
the prior distribution of the vector of branching times $\bm{t}$ is governed
by divergence parameter $c$ via the divergence function $a(t;c)$. Thus an
informative summary statistic for $c$ can be chosen by assessing its
information about $\bm{t}$. For example, tighter observed clusters indicate
small $c$ (e.g., $c<1$), where the level of tightness is indicated by the
branch lengths from leaves to their respective parents. We construct summary
statistics for $c$ based on a dendrogram estimated via hierarchical clustering
of $\mathbf{X}$ based on pairwise distances
$\delta_{i,i^{\prime}}:=\|\bm{X}_{i}-\bm{X}_{i^{\prime}}\|,i\neq i^{\prime}$.
The summary statistics $\bm{S}^{(c)}$ we choose is a ten-dimensional
concatenated vector comprising the 10th, 25th, 50th, 75th and 90th percentiles
of empirical distribution of: (i) $\delta_{i,i^{\prime}}$; (ii) branch lengths
associated with leaves of the dendrogram. Other candidate summary statistics
for $c$ are examined in Supplementary Material Section S4.1.

####  3.1.2 Stage 2: Sampling Tree Parameters $({\mathcal{T}},\bm{t})$ using
Metropolis-Hastings

For the second stage, we proceed by choosing a representative value
$(c_{0},\sigma^{2}_{0})$ chosen from the posterior sample of $(c,\sigma^{2})$,
which in our case is the posterior median. Then a Metropolis-Hastings (MH)
algorithm to sample from
$\textrm{pr}((\mathcal{T},\bm{t})|c_{0},\sigma^{2}_{0},\mathbf{X})$; recall
that the R${}_{\textrm{x}}$tree is characterized by both the topology
${\mathcal{T}}$ and divergence times $\bm{t}$. In particular, after
initialization (e.g., the dendrogram obtained from hierarchical clustering),
we first generate a candidate tree $(\mathcal{T}^{\prime},\bm{t}^{\prime})$
from the current tree $(\mathcal{T},\bm{t})$ in two steps: (i) detaching a
subtree from the original tree; (ii) reattaching the subtree back to the
remaining tree. Acceptance probabilities for a candidate tree can be computed
exactly and directly using the explicit likelihood in (5), without which they
would have to be calculated iteratively (Neal, 2003; Knowles and Ghahramani,
2015). See Supplementary Material Section S2.2 for details of the proposal
function and the acceptance probabilities.

######  Remark 1.

In order to use the explicit likelihood in (5) from Proposition 1 to generate
observed data $\mathbf{X}$, a tree-structured covariance
$\bm{\Sigma}^{\mathcal{T}}$ needs to be specified, whose entries in-turn
depend on the parameter $c$ through the divergence function. It is not
straightforward to fix or sample a $\bm{\Sigma}^{\mathcal{T}}$ since its
entries need to satisfy the inequalities (3). It is easier to generate data
$\mathbf{X}$ directly using the DDT generative mechanism in the ABC stage, and
this is the approach we follow and is described in Supplementary Section S2.

Summarily, there are three main advantages to using the explicit likelihood
from Proposition 1: (i) decoupling of Euclidean and tree parameters to enable
an efficient two-stage sampling algorithm; (ii) direct and exact computation
of tree acceptance probabilities in MH stage; (iii) determination of
informative sufficient statistic for $\sigma^{2}$ (Lemma 1).

###  3.2 Posterior Summary of R${}_{\textrm{x}}$-Tree,
$({\mathcal{T}},\bm{t})$

While quantifying uncertainty concerning the tree parameters
$({\mathcal{T}},\bm{t})$ is of main interest, we note that, from definition of
the DDT, this is influenced by uncertainty in the model parameters. In
particular, the first stage of ABC-MH produces weighted samples and we
calculate the posterior median by fitting an intercept-only quantile
regression with weights (see details in the Supplementary Material Section
S2.1). For the R${}_{\textrm{x}}$-tree, we consider global and local tree
posterior summaries that capture uncertainty in the latent hierarchy among all
and subsets of treatments.

Flexible posterior inference is readily available based on $L$ posterior
samples of $({\mathcal{T}},\bm{t})$ from the MH step. It is possible to
construct correspond tree-structured covariance matrices
$\bm{\Sigma}^{\mathcal{T}}$ from sample $({\mathcal{T}},\bm{t})$. Instead, we
compute:

  1. (a)

a global maximum a posteriori (MAP) estimate of the R${}_{\textrm{x}}$-tree
that represents the overall hierarchy underlying the treatment responses;

  2. (b)

local uncertainty estimates of co-clustering probabilities among a subset
$\mathcal{A}\subset\\{1,\ldots,I\\}$ of treatments based on posterior samples
of the corresponding subset of divergence times.

##### Posterior co-clustering probability functions.

We elaborate on the local summary (b). Suppose
$\mathcal{A}=\\{i,i^{\prime},i^{\prime\prime}\\}$ consists of three
treatments. Given a tree topology ${\mathcal{T}}$, note that at every
$t\in[0,1]$ a clustering of all $I$ treatments is available and the clustering
changes only at times $0<t_{1}<\cdots<t_{I-1}$. Consequently, for a given tree
topology ${\mathcal{T}}$ drawn from its posterior, we can compute for every
level $t\in[0,1]$ a posterior probability that $i,i^{\prime}$ and
$i^{\prime\prime}$ belong to the same cluster. Such a posterior probability
can be approximated using Monte Carlo on the $L$ posterior samples.
Accordingly, we define the estimated posterior co-clustering probability (PCP)
function associated with $\mathcal{A}$ as,

$\displaystyle\text{PCP}_{\mathcal{A}}(t)$
$\displaystyle=\frac{\sum_{l=1}^{L}\mathbb{I}_{[0,t^{(l)}_{i,i^{\prime},i^{\prime\prime}})}(t)}{L},$

where $\mathbb{I}_{B}$ is the indicator function on the set $B$ and
$t_{i,i^{\prime},i^{\prime\prime}}^{(l)}$ is the divergence time of
$\mathcal{A}=\\{i,i^{\prime},i^{\prime\prime}\\}$ in the $l$-th tree sample.
Essentially, the $\text{PCP}_{\mathcal{A}}(t)$ can be viewed as the proportion
of tree samples with $\\{i,i^{\prime},i^{\prime\prime}\\}$ having the most
recent common ancestor later than $t$.

For every subset $\mathcal{A}$, the function $[0,1]\ni
t\mapsto\text{PCP}_{\mathcal{A}}(t)\in[0,1]$ is non-increasing starting at 1
and ending at 0, and reveals propensity among treatments in $\mathcal{A}$ to
cluster as one traverses down an (estimate of) R${}_{\textrm{x}}$-tree
starting at the root: a curve that remains flat and drops quickly near 1
indicates higher relative similarity among the treatments in $\mathcal{A}$
relative to the rest of the treatments. A scalar summary of
$\text{PCP}_{\mathcal{A}}(t)$ is the area under its curve known as integrated
PCP $\text{iPCP}_{\mathcal{A}}$, which owing to the definition of
$\text{PCP}_{\mathcal{A}}(t)$, can be interpreted as the expected (or average)
chance of co-clustering for treatments in $\mathcal{A}$.

Figure 3 illustrates an example of a three-way $\text{iPCP}_{\mathcal{A}}$
with $\mathcal{A}=\\{i,i^{\prime},i{{}^{\prime\prime}}\\}$ for a PDX data with
$I$ treatments and $J$ patients (Figure 3(A)). Given $L=3$ posterior trees
samples (Figure 3(B)) drawn from the PDX data, we first calculate the whole
$\text{PCP}_{\mathcal{A}}(t)$ function by moving the time $t$ from $0$ to $1$.
Starting from time $t=0$, no treatment diverges at time $t=0$ and the
$\text{PCP}_{\mathcal{A}}(t)$ is $1$. At time $t^{\prime}$, treatments diverge
in one out of the three posterior trees and $\text{PCP}_{\mathcal{A}}(t)$
therefore drops from $1$ to $2/3$. Moving the time toward $t=1$, treatments
diverge in all trees and the $\text{PCP}_{\mathcal{A}}(t)$ drops to $0$. The
$\text{iPCP}_{\mathcal{A}}$ then can be obtained by the area under the
$\text{PCP}_{\mathcal{A}}(t)$.

Figure 3: Posterior tree summaries. (A) The input PDX data with $I$ treatments
and $J$ patients, and treatments
$\mathcal{A}=\\{i,i^{\prime},i{{}^{\prime\prime}}\\}$ are of interest. (B)
$\text{PCP}_{\mathcal{A}}(t)$ and $\text{iPCP}_{\mathcal{A}}$ for treatments
$\mathcal{A}$ based on $L=3$ posterior trees. The relevant divergence times
are represented by a “$\triangle$” in each posterior tree sample. For example,
at time $t^{\prime}$, the treatments in $\mathcal{A}$ diverge in one out of
the three trees. Because $\text{PCP}_{\mathcal{A}}(t^{\prime})$ is defined by
the proportion of posterior tree samples in which $\mathcal{A}$ has not
diverged up to and including $t^{\prime}$, it drops from $1$ to $2/3$.

######  Remark 2.

In the special case of $\mathcal{A}=\\{i,i^{\prime}\\}$ for two treatments,
the definition of iPCPA can be related to the cophenetic distance (Sokal and
Rohlf, 1962; Cardona et al., 2013) and, moreover, extends definition of the
cophenetic distance to multiple trees. Given two treatments $i$ and
$i^{\prime}$ in a single tree, let $t_{d}$ be the time at which their
corresponding Brownian paths diverge. Then
$\text{PCP}_{\mathcal{A}}(t)=\mathbb{I}_{[0,t_{d})}(t)$ and
$\text{iPCP}_{\mathcal{A}}=t_{d}$; this implies that the cophenetic distance
is $2(1-{t_{d}})$ and thus $\text{iPCP}_{\mathcal{A}}$ and the cophenetic
distances uniquely determines the same tree structure. For $L>1$ trees, a
Carlo average of divergence times of $L$ trees leads to the corresponding
$\text{iPCP}_{\mathcal{A}}$.

######  Remark 3.

Given $I$ treatments, since pairwise cophenetic distances from one tree
determines a tree (Lapointe and Legendre, 1991; McCullagh, 2006), one might
consider summarizing and represent posterior trees in terms of an $I\times I$
matrix $\bm{\Sigma}$ consisting of entries $\text{iPCP}_{\\{i,i^{\prime}\\}}$
for every pair of treatments of $(i,i^{\prime})$, estimated from the posterior
sample of trees. However, $\bm{\Sigma}$ need not to be a tree-structured
matrix that uniquely encodes a tree. It is possible to project $\bm{\Sigma}$
on to the space of tree-structured matrices (see for e.g., Bravo et al.
(2009)) but the projection might result in a non-binary tree structure. We
discuss this issue and its resolution in Supplementary Material Section S3.

##  4 Simulations

Accurate characterization of similarities among any subset of treatments is
central to our scientific interest in identifying the promising treatment
subsets for further investigation. In addition, we have introduced a two-stage
algorithm to improve our ability to efficiently draw tree samples from the
posterior distribution (similarly for the Euclidean parameters). To
demonstrate the modeling and computational advantages, we conduct two sets of
simulations. The first simulation shows that the proposed model estimates the
similarity (via iPCP) better than alternatives, even when the true data
generating mechanisms deviate from DDT assumptions in terms of the form of
divergence function, prior distribution for the unknown tree, and normality of
the responses. The second simulation illustrates the computational efficiency
of the proposed two-stage algorithm in producing higher quality posterior
samples of Euclidean parameters, resulting in more accurate subsequent
estimation of an unknown tree and iPCPs, two key quantities to our
interpretation of real data results.

###  4.1 Simulation I: Estimating Treatment Similarities

We first show that iPCPs estimated by DDT are closer to the true similarities
(operationalized by functions of elements in the true divergence times in
$\bm{\Sigma}^{\mathcal{T}}$) under different true data generating mechanisms
that may follow or deviate from the DDT model assumptions in three distinct
aspects (the form of divergence function, the prior distribution over the
unknown tree, and normality).

##### Simulation setup.

We simulate data by mimicking the PDX breast cancer data (see Section 5) with
$I=20$ treatments and $J=38$ patients. We set the true scale parameter as the
posterior median $\sigma_{0}^{2}$ and the true tree ${\mathcal{T}}_{0}$ as the
MAP tree that are estimated from the breast cancer data; We consider four
scenarios to represent different levels of deviation from the DDT model
assumptions:

  * (i)

No deviation of the true data generating mechanism from the fitted DDT models:
given $\sigma^{2}_{0}$ and ${\mathcal{T}}_{0}$, simulate data based on the DDT
marginal data distribution (Equation (5));

The true data generating mechanism deviates from the fitted DDT in terms of:

  * (ii)

divergence function: same as in (i), but the true tree is a random tree from
DDT with misspecified divergence function, $a(t;r)=\frac{r}{(1-t)^{2}},r=0.5$;

  * (iii)

prior for tree topology: same as in (i), but the true tree is a random tree
from the coalescence model (generated by function `rcoal` in R package `ape`),
and,

  * (iv)

marginal data distribution: same as in (i), but the marginal likelihood is a
centered multivariate $t$ distribution with degree-of-freedom four and scaled
by $\sigma^{2}_{0}\bm{\Sigma}^{\mathcal{T}_{0}}$.

For each of four true data generating mechanisms above, we simulate $B=50$
replicate data sets. In the following, we use the DDT model and the two-stage
algorithm for all estimation regardless of the true data generating
mechanisms. For DDT, we ran the two-stage algorithm where the second stage is
implemented with five parallel chains. For each chain, we ran $10,000$
iterations, discarded first $9,000$ trees and combined five chains with a
total of $5,000$ posterior tree samples.

First, we compute the iPCPs for all pairs of treatment combinations following
the definition of $\textrm{iPCP}_{\mathcal{A}}$ where
$\mathcal{A}=\\{i,i^{\prime}\\},1\leq i<i^{\prime}\leq I$. Two alternative
approaches to defining and estimating similarities between treatments are
considered: (i) similarity derived from agglomerative hierarchical clustering,
and (ii) empirical Pearson correlation of the two vectors of responses
$\bm{X}_{i}$ and $\bm{X}_{i^{\prime}}$, for $i\neq i^{\prime}$. In particular,
for (i), we considered five different linkage methods (Ward, Ward’s D2,
single, complete and Mcquitty) with Euclidean distances. Given an estimated
dendrogram from hierarchical clustering, the similarity for a pair of
treatments is defined by first normalizing the sum of branch lengths from the
root to leaf as $1$, and then calculating the area under of the co-clustering
curve (AUC) obtained by cutting the dendrogram at various levels from 0 to 1.
For three- or higher-way comparisons, (i) can still produce an AUC based on a
dendrogram obtained from hierarchical clustering, while the empirical Pearson
correlation in (ii) is undefined hence not viable as a comparator beyond
assessing pairwise treatment similarities.

##### Performance metrics.

For treatment pairs $\mathcal{A}=\\{i,i^{\prime}\\}$, to assess the quality of
estimated treatment similarities for each of the methods above (DDT-based,
hierarchical-clustering-based, and empirical Pearson correlation), we compare
the estimated values against the true branching time
$\Sigma_{i,i^{\prime}}^{{\mathcal{T}}_{0}}$; similarly when assessing recovery
of three-way treatment similarities, e.g.,
$\mathcal{A}=\\{i,i^{\prime},i^{\prime\prime}\\}$,
$\Sigma^{{\mathcal{T}}_{0}}_{i,i^{\prime},i^{\prime\prime}}$ is defined as the
time when $\\{i,i^{\prime},i^{\prime\prime}\\}$ first branches in the true
tree ${\mathcal{T}}_{0}$. In particular, for replication data set
$b=1,\ldots,B$, let $\widehat{\Sigma}^{(b)}_{i,i^{\prime}}$ generically
represent the pairwise similarities for treatment subsets $(i,i^{\prime})$
that can be based on DDT, hierarchical clustering or empirical pairwise
Pearson correlation. For three-way comparisons, let
$\widehat{\Sigma}^{(b)}_{i,i^{\prime},i^{\prime\prime}}$ generically represent
the three-way similarities for treatment subset
$(i,i^{\prime},i^{\prime\prime})$ that can be based on DDT, or hierarchical
clustering.

We assess the goodness of recovery by computing
$\sqrt{\sum_{i,i^{\prime}}(\widehat{\Sigma}^{(b)}_{i,i^{\prime}}-\Sigma^{{\mathcal{T}}_{0}}_{i,i^{\prime}})^{2}}$,
the Frobenious norm of the matrix in recovering the entire
$\Sigma^{{\mathcal{T}}_{0}}$. We compute
$\max_{i,i^{\prime},i^{\prime\prime}}|\widehat{\Sigma}^{(b)}_{i,i^{\prime},i^{\prime\prime}}-\Sigma^{{\mathcal{T}}_{0}}_{i,i^{\prime},i^{\prime\prime}}|$,
the max-norm of the matrix in recovering the true three-way similarities. For
a given method and treatment subset $\mathcal{A}$, the above procedure results
in $B$ values, the distribution of which can be compared across methods;
smaller values indicate better recovery of the true similarities.

Alternatively, for each method and each treatment subset, we also compute the
Pearson correlation between the estimated similarities and the true branching
times across replicates for pairwise or three-way treatment subsets:
$\widehat{\text{Cor}}\left((\widehat{\Sigma}^{(b)}_{i,i^{\prime}},\Sigma^{{\mathcal{T}}_{0}}_{i,i^{\prime}}),b=1,\ldots,B=50\right)$,
for treatments $i<i^{\prime}$ and
$\widehat{\text{Cor}}\left((\widehat{\Sigma}^{(b)}_{i,i^{\prime},i^{\prime\prime}},\Sigma^{{\mathcal{T}}_{0}}_{i,i^{\prime},i^{\prime\prime}}),b=1,\ldots,B=50\right)$,
for treatments $i<i^{\prime}<i^{\prime\prime}$. We refer to this metric as
“Correlation of correlations” (the latter uses the fact that the entries in
the true $\Sigma^{{\mathcal{T}}_{0}}$ being correlations; see Equation (5));
higher values indicate better recovery of the true similarities.

##### Simulation results.

We observe that DDT better estimates the treatment similarities even under
misspecified models. In particular, under scenarios where the true data
generating mechanisms deviate from the fitted DDT model assumptions (ii-iv),
the DDT captures the true pairwise and three-way treatment similarities the
best by higher values in correlation of correlations (left panels, Figure 4)
and lower matrix/array distances (right panels, Figure 4). In particular, the
fitted DDT with divergence function $a(t)=c/(1-t)$ under Scenario i, ii and
iii performed similarly well indicating the relative insensitivity to the DDT
modeling assumptions with respect to divergence function and the tree
generative model. Under Scenario iv where the marginal likelihood assumption
deviates from Gaussian with heavier tails, the similarity estimates from all
methods deteriorate relative to Scenarios i-iii. Comparing between methods,
the similarities derived from hierarchical clustering with single linkage is
comparable to DDT model when evaluated by correlation of correlation, but
worse than DDT when evaluated by the matrix norm.

Figure 4: Simulation studies for comparing the quality of estimated treatment
similarities based on DDT, hierarchical clustering, and empirical Pearson
correlation. Two performance metrics are used: (Left) Correlation of
correlation (higher values are better); (Right) Matrix distances with
Frobenius norm for pairwise similarity and max norm for three-way similarity
(lower values are better). DDT captures both true pairwise (upper panels) and
three-way (lower panels) similarity best under four levels of misspecification
scenarios.

###  4.2 Simulation II: Comparison with Single-Stage MCMC Algorithms

We have also conducted extensive simulation studies that focus on the
computational aspect of the proposed algorithms and demonstrate the advantage
of the proposed two-stage algorithm in producing higher quality posterior
samples of the unknown tree than classical single-stage MCMC algorithms. In
particular, we demonstrate that the proposed algorithm produces (i) MAP trees
that are closer to the true tree than alternatives (hierarchical clustering,
single-stage MH with default hierarchical clustering or the true tree at
initialization) and (ii) more accurate estimation of pairwise treatment
similarities compared to single-stage MCMC algorithms. See Supplementary
Material Section S5 for further details.

##### Additional simulations and sensitivity analyses.

Aside from the simulations above focusing on the tree structure and the
divergence time, Supplementary Material S4 offers additional details for
Euclidean parameters including the parameter inference and algorithm
diagnostics. In particular, we empirically show that current ${\bm{S}}^{(c)}$
and $S^{(\sigma^{2})}$ outperform other candidate summary statistics in terms
of bias in Section S4.1. In Section S4.2, we present additional simulation
results that demonstrate that the two-stage algorithm (i) enjoys stable
effective sample size (ESS) for $(c,\sigma^{2})$; (ii) leads to similar or
better inference on $(c,\sigma^{2})$, as ascertained using credible intervals.
Section S4.3 checks the convergence of MH and the goodness of fit for ABC .

##  5 Treatment Trees in Cancer using PDX Data

###  5.1 Dataset Overview and Key Scientific Questions

We leverage a recently collated PDX dataset from the Novartis Institutes for
BioMedical Research - PDX Encyclopedia [NIBR-PDXE, (Gao et al., 2015)] that
interrogated multiple targeted therapies across different cancers and
established that PDX systems provide a more accurate measure of the response
of a population of patients than traditional preclinical models. Briefly, the
NIBR-PDXE consists of $>1,000$ PDX lines across a range of human cancers and
uses a $1\times 1\times 1$ design (one animal per PDX model per treatment);
i.e., each PDX line from a given patient was treated simultaneously with
multiple treatments allowing for direct assessments of treatment hierarchies
and responses. In this paper, we focus on our analyses on a subset of PDX
lines with complete responses across five common human cancers: Breast cancer
(BRCA), Cutaneous Melanoma (CM, skin cancer), Colorectal cancer (CRC), Non-
small Cell Lung Carcinoma (NSCLC), and Pancreatic Ductal Adenocarcinoma
(PDAC). After re-scaling data and missing data imputation, different numbers
of treatments, $I$, and PDX models, $J$, presented in the five cancers were,
$(I,J)$: BRCA, $(20,38)$; CRC, $(20,40)$; CM, $(14,32)$; NSCLC, $(21,25)$; and
PDAC, $(20,36)$. (See Supplementary Material Table S7 for treatment names and
Section S6.1 for details of pre-processing procedures.)

In our analysis, we used the best average response (BAR) as the main response,
by taking the untreated group as the reference group and using the tumor size
difference before and after administration of the treatment(s) following
Rashid et al. (2020). Positive values of BAR indicate the treatment(s) shrunk
the tumor more than the untreated group with higher values indicative of
(higher) treatment efficacy. The treatments included both drugs administered
individually with established mechanisms (referred to as “monotherapy”) and
multiple drugs combined with potentially unknown synergistic effects (referred
to as “combination therapy”). Our key scientific questions were as follows:
(a) identify plausible biological mechanisms that characterize treatment
responses for monotherapies within and between cancers; (b) evaluate the
effectiveness of combination therapies based on biological mechanisms. Due to
a potentially better outcome and lower resistance, combination therapy with
synergistic mechanism is highly desirable (Bayat Mokhtari et al., 2017).

DDT model setup. For all analyses we followed the setup in the Section 4.1 and
obtained $N^{\textrm{syn}}=600,000$ synthetic datasets from the ABC algorithm
(Section 3.1.1) with prior $c\sim{\sf Gamma}(2,2)$ and $1/\sigma^{2}\sim{\sf
Gamma}(1,1)$ and took the first 0.5% ($d=0.5\%$) closest data in terms of
$\bm{S}^{(c)}$ and $S^{(\sigma^{2})}$. We calculated the posterior median of
$(c,\sigma^{2})$ as described in Section 3.2. For the second-stage MH, we ran
five chains of the two-stage algorithm with $(c,\sigma^{2})$ fixed at the
posterior median by 10,000 iterations and discarded the first 9,000 trees,
which resulted in 5,000 posterior trees in total. Finally, we calculated the
R${}_{\textrm{x}}$-tree (MAP) and iPCP based on 5,000 posterior trees for all
subsequent analyses and interpretations. All computations were divided on
multiple different CPUs (see the Supplementary Table S5 for the full list of
CPUs). For the BRCA data with $I=20$ and $J=38$, we divided the ABC stage into
$34$ compute cores with a total of $141$ CPU hours and maximum $4.7$ hours in
real time. For the MH stage and the single-stage MCMC, we split the
computation on $5$ compute cores with a total of $8.6$ and $12$ CPU hours, and
a maximum $1.7$ and $2.5$ hours in real time, respectively.

Our results are organized as follows: we provide a summary of the
R${}_{\textrm{x}}$-tree estimation and treatment clusters in Section 5.2
followed by specific biological and translational interpretations in Sections
5.3 and 5.4 for monotherapy and combination therapy, respectively. Additional
results can be accessed and visualized using our companion R-shiny application
(see Supplementary Material Section S6.4 for details).

###  5.2 R${}_{\textrm{x}}$-tree Estimation and Treatment Clusters

We focus our discussion on three cancers: BRCA, CRC and CM here – see
Supplementary Materials Section S6.3 for NSCLC and PDAC. In Figure 5,
R${}_{\textrm{x}}$-tree, pairwise iPCP and (scaled) Pearson correlation are
shown in the left, middle and right panels, respectively. Focusing on the left
two panels, we observe that the R${}_{\textrm{x}}$-tree and the pairwise iPCP
matrix show the similar clustering patterns. For example, three combination
therapies in CM form a tight subtree and are labeled by a box in the
R${}_{\textrm{x}}$-tree of Figure 5 and a block with higher values of iPCP
among three combination therapies also shows up in the corresponding iPCP
matrix with a box labeled. In our analysis, the treatments predominantly
target six oncogenic pathways that are closely related to the cell
proliferation and cell cycle: (i) phosphoinositide 3-kinases, PI3K; (ii)
mitogen-activated protein kinases, MAPK; (iii) cyclin-dependent kinases, CDK;
(iv) murine double minute 2, MDM2; (v) janus kinase, JAK; (vi)
serine/threonine-protein kinase B-Raf, BRAF. We label targeting pathways above
for monotherapies with solid dots and further group PI3K, MAPK and CDK due to
the common downstream mechanisms (e.g., Repetto et al., 2018; Kurtzeborn et
al., 2019). Roughly, the R${}_{\textrm{x}}$-tree from our model clusters
monotherapies targeting oncogenic processes above and largely agrees with
common and established biology mechanisms. For example, all PI3K-MAPK-CDK
inhibitors (solid square) belong to a tighter subtree in three cancers; two
MDM2 monotherapies (solid triangle) are closest in both BRCA and CRC. While
visual inspection of the MAP R${}_{\textrm{x}}$-tree agrees with known
biology, iPCP further quantifies the similarity by assimilating the
information across multiple trees from our MCMC samples. For the ensuing
interpretations in Sections 5.3 and 5.4, we focus on iPCP and verify our model
through monotherapies with known biology, since our a priori hypothesis is
that monotherapies that share the same downstream pathways should exhibit
higher iPCP values. Furthermore, we extend our work to identify combination
therapies with synergy and discover several combination therapies for each
cancer.

Figure 5: The R${}_{\textrm{x}}$-tree and iPCP for breast cancer (BRCA, top
row), colorectal cancer (CRC, middle row) and melanoma (CM, lower row). Three
panels in each row represent: (left) estimated R${}_{\textrm{x}}$-tree (MAP);
distinct external target pathway information is shown in distinct shapes for
groups of treatments on the leaves; (middle) estimated pairwise iPCP, i.e.,
the posterior mean divergence time for pairs of entities on the leaves (see
the result paragraph for definition for any subset of entities); (right)
scaled Pearson correlation for each pair of treatments. Note that the MAP
visualizes the hierarchy among treatments; the iPCP is not calculated based on
the MAP, but based on posterior tree samples (see definition in Section 3.2)

###  5.3 Biological Mechanisms in Monotherapy

Our estimation procedure exhibits a high level of concordance between known
biological mechanisms and established monotherapies for multiple key
signalling pathways. From the R${}_{\textrm{x}}$-tree in Figure 5, aside from
the oncogenic process (solid dots) introduced above, monotherapies also target
receptors (hollow circles) or other non-kinase targets (e.g. tubulin;
crosses). We summarize our key findings and interpretations for each of
signaling pathways and their regulatory axes, namely, PI3K-MAPK-CDK and MDM2
from cell cycle regulatory pathways, human epidermal growth factor receptor 3
(ERBB3) from receptor pathways, and tubulin from non-kinase pathways along
with their implications in monotherapy across different cancers.

PI3K-MAPK-CDK inhibitors. For treatments targeting PI3K, MAPK and CDK,
treatments have the same target share high iPCP. In the NIBR-PDXE dataset,
three PI3K inhibitors (BKM120, BYL719 and CLR457), two MAPK inhibitors
(binimetinib and CKX620) and one CDK inhibitor (LEE011) were tested, but
different cancers contain different numbers of treatments. Specifically, all
three PI3K inhibitors present in BRCA and CRC, but only BKM120 is tested in
CM; CRC contains two MAPK inhibitors while BRCA and CM only have binimetinib;
LEE011 is tested in all three cancers. In Figure 6, BKM120, BYL719 and CLR457
share high pairwise iPCPs (box (1)) and all target PI3K for BRCA and CRC
(BRCA, (BKM120, CLR457): 0.8986, (BKM120, BYL719): 0.8002, (BYL719, CLR457):
0.8002; CRC, (BKM120, CLR457): 0.7555, (BKM120, BYL719): 0.8041, (BYL719,
CLR457): 0.7597); MAPK (box (2)) inhibitors, binimetinib and CKX620, show a
high pairwise iPCP in CRC (0.7792). Asides from the pairwise iPCPs, our model
also suggests high multi-way iPCPs among PI3K inhibitors in BRCA (0.8002) and
CRC (0.7513). Among these inhibitors, PI3K inhibitor of BYL719 was approved by
FDA for breast cancer; MAPK inhibitor of binimetinib was approved by FDA for
BRAF mutant melanoma in combination with encorafenib; and CDK inhibitor of
LEE011 was approved for breast cancer.

Our model suggests treatments targeting different pathways also share high
iPCP values across different cancers. Monotherapies targeting different cell
cycle regulatory pathways (PI3K, MAPK and CDK) exhibit high iPCPs. CDK
inhibitor, LEE011, and MAPK inhibitors share high pairwise iPCP values in BRCA
((LEE011, binimetinib): 0.7709), CRC ((LEE011, binimetinib): 0.8617, (LEE011,
CKX620): 0.7820) and CM ((LEE011, binimetinib): 0.8210) in the Figure 6 with
box (3). High iPCP among MAPK and CDK inhibitors agree with biology, since it
is known that CDK and MAPK collaboratively regulate downstream pathways such
as Ste5 (Repetto et al., 2018). High pairwise iPCP values between PI3K and
MAPK inhibitors were observed in box (3) in the Figure 6. Specifically, our
model suggests high pairwise iPCPs as follows: (i) BRCA, (binimetinib,
BKM120): 0.7427, (binimetinib, BYL719): 0.7441, (binimetinib, CLR457):
0.7427)); (ii) CRC, (binimetinib, BKM120): 0.7374, (binimetinib, BYL719):
0.7388, (binimetinib, CLR457): 0.7541, (CKX620, BKM120): 0.7366, (CKX620,
BYL719): 0.7357, (CKX620, CLR457): 0.7676)); (iii) CM, (binimetinib, BKM120):
0.8882. Aside from the pairwise iPCPs above, high multi-way iPCPs in BRCA
(0.7422), CRC (0.7300) and CM (0.8882) also show the similar information. From
the existing literature, both PI3K and MAPK can be induced by ERBB3
phosphorylation (Balko et al., 2012) and it is not surprising to see high
iPCPs between PI3K and MAPK inhibitors.

ERBB3 and tubulin inhibitors. Our model also found high iPCP values among
ERBB3, tubulin and PI3K-MAPK-CDK inhibitors in BRCA. ERRB3 inhibitor, LJM716,
exhibits high pairwise iPCP values with PI3K (BKM120: 0.7501, BYL719: 0.7513,
CLR457: 0.7500), MAPK (binimetinib: 0.7811), CDK (LEE011: 0.7847) and tubulin
(paclitaxel: 0.7505) inhibitors in the Figure 6 with box (5). Since PI3K and
MAPK are downstream pathways of ERBB3 (Balko et al., 2012) and CDK works
closely with PI3K and MAPK (Kurtzeborn et al., 2019; Repetto et al., 2018),
high iPCPs between ERBB3 inhibitor and PI3K-MAPK-CDK inhibitors are not
surprising. For ERBB3 and tubulin, ERBB3 is a critical regulator of
microtubule assembly (Wu et al., 2021) and tubulin plays an important role in
building microtubules. Since microtubules form the skeletons of cells and are
essential for cell division (Gunning et al., 2015; Haider et al., 2019),
tubulin inhibitor, paclitaxel, kills cancer cell by interfering cell division
and is an FDA-approved treatment. In congruence with the above results,
tubulin inhibitor paclitaxel also shares high iPCPs with PI3K (BKM120: 0.8076,
BYL719: 0.8063, CLR457: 0.8076), MAPK (binimetinib: 0.7433), CDK (LEE011:
0.7587) and ERBB3 (LJM716: 0.7505) in the Figure 6 with box (5). In addition,
another CDK4 inhibitor BPT also inhibits tubulin (Mahale et al., 2015) and
PI3K inhibitor BKM120 inhibits the formation of microtubule (Bohnacker et al.,
2017). Both offer additional reasons for high iPCP between tubulin and
PI3K-MAPK-CDK inhibitors.

MDM2 inhibitors. We found two drugs: CGM097 and HDM201 share high iPCP values
in BRCA (0.8365) and CRC (0.7860) in the Figure 6 with box (4). Since CGM097
and HDM201 target the same pathway, MDM2, high iPCPs suggest a high similarity
between CGM097 and HDM201 and show consistent results between our model and
underlying biological mechanism. MDM2 negatively regulates the tumor
suppressor, p53 (Zhao et al., 2014) and if MDM2 is suppressed by inhibitors,
p53 is able to prevent tumor formation. Both CGM097 and HDM201 entered phase I
clinical trial (Konopleva et al., 2020) for wild-type p53 solid tumors and
leukemia, respectively.

Figure 6: Bar plot of iPCPs for pairs of combination therapies (red bars) and
pairs of monotherapies (green bars): (A) breast cancer, (B) melanoma and (C)
colorectal cancer. The bar plots are sorted by the iPCP values (high to low);
pairs of treatments are shown only if the estimated iPCP is greater than
$0.7$. Monotherapies have different known targets which are listed in the
bottom-right table (see Section 5.3 for more details and discussion on
monotherapies).

###  5.4 Implications in Combination Therapy

Based on the concordance between the monotherapy and biology mechanism, we
further investigate combination therapies to identify mechanisms with
synergistic effect. In NIBR-PDXE, $21$ combination therapies were tested and
only one of them includes three monotherapies (BYL719 + cetuximab +
encorafenib in CRC) and the rest contain two monotherapies. Out of $21$
combination therapies, only three do not target any cell cycle (PI3K, MAPK,
CDK, MDM2, JAK and BRAF) pathways (see Supplementary Material Table S8 for the
full list of combination therapies). From the R${}_{\textrm{x}}$-tree in
Figure 5, combination therapies tend to form a tighter subtree and are closer
to monotherapies targeting PI3K-MAPK-CDK, which implies that the mechanisms
under combination therapies are similar to each other and are closer to the
PI3K-MAPK-CDK pathways. We identified several combination therapies with known
synergistic effects and provide a brief description for each of the cancers in
the following paragraphs.

Breast cancer. Four combination therapies were tested in BRCA and three
therapies targeting PI3K-MAPK-CDK (BYL719 + LJM716, BYL719 + LEE011 and LEE011
+ everolimus) form a subtree in R${}_{\textrm{x}}$-tree with a high three-way
iPCP (0.8719). Among these combination therapies, PI3K-CDK inhibitor, BYL719 +
LEE011, suggests a possible synergistic regulation (Vora et al., 2014; Bonelli
et al., 2017; Yuan et al., 2019). Based on the high iPCP between BYL719 +
LEE011 and the rest two therapies, we suggest synergistic effect for
combination therapies targeting PI3K-ERBB3 (BYL719 + LJM716), and CDK-MTOR
(LEE011 + everolimus) for future investigation.

Colorectal cancer. Our model suggests a high three-way iPCP (0.7437) among
PI3K-EGFR (BYL719 + cetuximab), PI3K-EGFR-BRAF (BYL719 + cetuximab +
encorafenib) and PI3K-ERBB3 (BYL719 + LJM716) inhibitors. Since the triple
therapy (BYL719 + cetuximab + encorafenib) enters the phase I clinical trial
with synergy (Geel et al., 2014), our model proposes the potential synergistic
effect for PI3K-ERBB3 based on iPCP for future investigation. Of note, we
found a modest iPCP (0.6280) between the FDA-approved combination therapy
EGFR-BRAF (cetuximab + encorafenib) and PI3K-EGFR-BRAF (BYL719 + cetuximab +
encorafenib) and the modest iPCP can be explained by an additional drug-drug
interaction between BYL719 and encorafenib in triple-combined therapy (van
Geel et al., 2017).

Melanoma. In NIBR-PDXE, three combination therapies were tested in CM, and all
of them consist one monotherapy targeting PI3K-MAPK-CDK and the other one
targeting BRAF. A tight subtree is observed in the R${}_{\textrm{x}}$-tree and
our model also suggests a high iPCP (0.9222) among three combination
therapies. Since PI3K, MAPK and CDK work closely and share a high iPCP
(0.8204) among monotherapies in CM, a high iPCP (0.9222) among three
combination therapies is not surprising. Since two combination therapies of
BRAF-MAPK (dabrafenib + trametinib and encorafenib + binimetinib) are approved
by FDA for BRAF-mutant metastatic melanoma (Dummer et al., 2018a, b; Robert et
al., 2019), we recommend the synergy for BRAF-PI3K (encorafenib + BKM120) and
BRAF-CDK (encorafenib + LEE011) inhibitors.

##### Comparison to alternative approaches.

Unlike the probabilistic generative modeling approach proposed in this paper,
standard distance-based agglomerative hierarchical clustering and Pearson
correlation can also be applied to the PDX data to estimate the similarity.
However, simple pairwise similarities can be potentially noisy and the
uncertainty in the estimation is not fully incorporated due to the absence of
a generative model. As we showed in the Section 4.1 (Simulation I) that
agglomerative hierarchical clustering and the Pearson correlation leads to
inferior recovery of the true branching times and the true tree structure
under different data generating mechanisms mimicking the real data. As further
evidence, we compute pairwise similarities based on Pearson correlation (other
distance metrics show similar patterns) in the right panel of Figure 5. By
mapping the original Pearson correlation $\rho\in[-1,1]$ through a linear
function $\frac{\rho+1}{2}$, we make the range of iPCP and Pearson correlation
comparable. We observe that pairwise iPCP estimated through the DDT model is
less noisy than Pearson correlation. For example, both iPCP and Pearson
correlation in CM show higher similarities among combination therapy framed by
a box, but iPCP exhibits a clearer pattern than Pearson correlation.

##  6 Summary and Discussion

In translational oncology research, PDX studies have emerged as a unique study
design that evaluates multiple treatments when applied to samples from the
same human tumor implanted into genetically identical mice. PDX systems are
promising tools for large-scale screening to evaluate a large number of FDA-
approved and novel cancer therapies. However, there remain scientific
questions concerning how distinct treatments may be synergistic in inducing
similar efficacious responses, and how to identify promising subsets of
treatments for further clinical evaluation. To this end, in this paper, we
propose a probabilistic framework to learn treatment trees
(R${}_{\textrm{x}}$-trees) from PDX data to identify promising treatment
combinations and plausible biological mechanisms that confer synergistic
effect(s). In particular, in a Bayesian framework based on the Dirichlet
Diffusion Tree, we estimate a maximum a posteriori rooted binary tree with the
treatments on the leaves and propose a posterior uncertainty-aware similarity
measure (iPCP) for any subset of treatments. The divergence times of the DDT
encode the tree topology and are profitably interpreted within the context of
an underlying plausible biological mechanism of treatment actions.

From the class of probabilistic models with an unknown tree structure
component, we have chosen the DDT mainly owing to the availability of a
closed-form marginal likelihood that directly links the tree topological
structure to the covariance structure of the observed PDX data, which
additionally decouples the Euclidean and tree parameters; to the best of our
knowledge this method has not been proposed or explored hitherto for the DDT.
The decoupling leads to efficient posterior inference via a two-stage
algorithm that confers several advantages. The algorithm generates posterior
samples of Euclidean parameters through approximate Bayesian computation and
passes the posterior medians to a second stage classical Metropolis-Hastings
algorithm for sampling from the conditional posterior distribution of the tree
given all other quantities. Through simulation studies, we show that the
proposed two-stage algorithm generates better posterior tree samples and
captures the true similarity among treatments better than alternatives such as
single-stage MCMC and naive Pearson correlations. The posterior samples of
trees are summarized by iPCP, which we propose to measure the empirical
mechanistic similarity for multiple treatments incorporating uncertainty.

Using the proposed methodology on NIBR-PDXE data, we estimate
R${}_{\textrm{x}}$-trees and iPCPs for five cancers. Among the monotherapies,
iPCP is highly concordant with known biology across different cancers. For
example, CGM097 and HDM201 show a high iPCP value among treatments in breast
cancer and melanoma, which corroborates known mechanisms, since both
monotherapies target the same biological pathway, MDM2, and have been further
validated in recent clinical trials (Konopleva et al., 2020). The proposed
iPCP can also suggest improvements upon an existing combination therapy. We
first identify a combination therapy with known synergy (not based on the our
data) and then determine which additional therapies (monotherapies or
combination therapies) have high iPCPs when considered together with the
existing combination therapy. Based on the NIBR-PDXE data, for each cancer, we
suggest potential synergies between PI3K-ERBB3 and CDK-MTOR for breast cancer,
PI3K-ERBB3 for colorectal cancer, and BRAF-PI3K and BRAF-CDK for melanoma that
could be potentially explored in future translational studies.

The present analysis shows the promise in guiding the potential treatment
strategies by building trees upon NIBR-PDXE dataset using treatment responses,
but assume independent patients without using the underlying patient-specific
genomic information that is available in the NIBR-PDXE. By including patient-
specific genomic information, we may further improve our ability to identify
synergistic treatments that may be specific to a subset of patients. One
approach to utilize genomic information could be to extend the DDT model to
incorporate patient-specific genomic information in the mean structure or the
column covariance of the marginal likelihood of Equation (4). Models with non-
Gaussian marginal likelihood and non-binary treatment tree in principle can be
defined by considering generative tree models based on general diffusion
processes (Heaukulani et al., 2014; Knowles and Ghahramani, 2015). Both
extensions raise significant, non-trivial methodological and computation
issues (e.g., deriving tractable likelihoods; finding low-dimensional summary
statistics for new parameters) and constitute the foundation for future work.

##### Code and data availability

We also provide a general purpose code in R that accompanies this manuscript
along with all the necessary documentation and datasets required to replicate
our results (see https://github.com/bayesrx/RxTree). Furthermore, to aid
access and visualization of the results, we have also developed an R-shiny
application (see Supplementary Material Section S6.4).

## References

  * Ferlay et al. [2020] J Ferlay, M Ervik, F Lam, M Colombet, L Mery, M Piñeros, A Znaor, I Soerjomataram, and F Bray.  _Global Cancer Observatory: Cancer Today_.  Lyon, France: International Agency for Research on Cancer, 2020.  Available from: https://gco.iarc.fr/today, accessed 05.28.2021.
  * Dagogo-Jack and Shaw [2018] I. Dagogo-Jack and A. T. Shaw.  Tumour heterogeneity and resistance to cancer therapies.  _Nat Rev Clin Oncol_ , 15(2):81–94, 02 2018.
  * Groisberg and Subbiah [2021] R. Groisberg and V. Subbiah.  Combination therapies for precision oncology: the ultimate whack-a-mole game.  _Clin Cancer Res_ , 27(10):2672–2674, May 2021\.
  * Sawyers [2013] C. L. Sawyers.  Perspective: combined forces.  _Nature_ , 498(7455):S7, Jun 2013.
  * Bayat Mokhtari et al. [2017] R. Bayat Mokhtari, T. S. Homayouni, N. Baluch, E. Morgatskaya, S. Kumar, B. Das, and H. Yeger.  Combination therapy in combating cancer.  _Oncotarget_ , 8(23):38022–38043, Jun 2017.
  * Sun et al. [2016] W. Sun, P. E. Sanderson, and W. Zheng.  Drug combination therapy increases successful drug repositioning.  _Drug Discov Today_ , 21(7):1189–1195, 07 2016\.
  * Tentler et al. [2012] J. J. Tentler, A. C. Tan, C. D. Weekes, A. Jimeno, S. Leong, T. M. Pitts, J. J. Arcaroli, W. A. Messersmith, and S. G. Eckhardt.  Patient-derived tumour xenografts as models for oncology drug development.  _Nat Rev Clin Oncol_ , 9(6):338–350, Apr 2012\.
  * Bhimani et al. [2020] J. Bhimani, K. Ball, and J. Stebbing.  Patient-derived xenograft models-the future of personalised cancer treatment.  _Br J Cancer_ , 122(5):601–602, 03 2020.
  * Hidalgo et al. [2014] M. Hidalgo, F. Amant, A. V. Biankin, E. Budinská, A. T. Byrne, C. Caldas, R. B. Clarke, S. de Jong, J. Jonkers, G. M. Mælandsmo, S. Roman-Roman, J. Seoane, L. Trusolino, and A. Villanueva.  Patient-derived xenograft models: an emerging platform for translational cancer research.  _Cancer Discov_ , 4(9):998–1013, Sep 2014.
  * Lai et al. [2017] Y. Lai, X. Wei, S. Lin, L. Qin, L. Cheng, and P. Li.  Current status and perspectives of patient-derived xenograft models in cancer research.  _J Hematol Oncol_ , 10(1):106, 05 2017.
  * Yoshida [2020] Go J. Yoshida.  Applications of patient-derived tumor xenograft models and tumor organoids.  _Journal of Hematology & Oncology_, 13(1):4, Jan 2020.  ISSN 1756-8722.  doi:10.1186/s13045-019-0829-z.  URL https://doi.org/10.1186/s13045-019-0829-z.
  * Topp et al. [2014] Monique D. Topp, Lynne Hartley, Michele Cook, Valerie Heong, Emma Boehm, Lauren McShane, Jan Pyman, Orla McNally, Sumitra Ananda, Marisol Harrell, Dariush Etemadmoghadam, Laura Galletta, Kathryn Alsop, Gillian Mitchell, Stephen B. Fox, Jeffrey B. Kerr, Karla J. Hutt, Scott H. Kaufmann, Australian Ovarian Cancer Study, Elizabeth M. Swisher, David D. Bowtell, Matthew J. Wakefield, and Clare L. Scott.  Molecular correlates of platinum response in human high-grade serous ovarian cancer patient-derived xenografts.  _Molecular Oncology_ , 8(3):656–668, 2014.  ISSN 1574-7891.  doi:https://doi.org/10.1016/j.molonc.2014.01.008.  URL https://www.sciencedirect.com/science/article/pii/S1574789114000209.
  * Nunes et al. [2015] Manoel Nunes, Patricia Vrignaud, Sophie Vacher, Sophie Richon, Astrid Lièvre, Wulfran Cacheux, Louis-Bastien Weiswald, Gerald Massonnet, Sophie Chateau-Joubert, André Nicolas, Colette Dib, Weidong Zhang, James Watters, Donald Bergstrom, Sergio Roman-Roman, Ivan Bièche, and Virginie Dangles-Marie.  Evaluating patient-derived colorectal cancer xenografts as preclinical models by comparison with patient clinical data.  _Cancer Research_ , 75(8):1560–1566, 2015.  ISSN 0008-5472.  doi:10.1158/0008-5472.CAN-14-1590.  URL https://cancerres.aacrjournals.org/content/75/8/1560.
  * Clohessy and Pandolfi [2015] J. G. Clohessy and P. P. Pandolfi.  Mouse hospital and co-clinical trial project–from bench to bedside.  _Nat Rev Clin Oncol_ , 12(8):491–498, Aug 2015\.
  * Lunardi and Pandolfi [2015] A. Lunardi and P. P. Pandolfi.  A co-clinical platform to accelerate cancer treatment optimization.  _Trends Mol Med_ , 21(1):1–5, Jan 2015.
  * Grant et al. [2010] R.L. Grant, A.B. Combs, and D. Acosta.  Experimental models for the investigation of toxicological mechanisms.  In Charlene A. McQueen, editor, _Comprehensive Toxicology (Second Edition)_ , pages 203–224. Elsevier, Oxford, second edition edition, 2010\.  ISBN 978-0-08-046884-6.  doi:https://doi.org/10.1016/B978-0-08-046884-6.00110-X.  URL https://www.sciencedirect.com/science/article/pii/B978008046884600110X.
  * Krumbach et al. [2011] R. Krumbach, J. Schüler, M. Hofmann, T. Giesemann, H. H. Fiebig, and T. Beckers.  Primary resistance to cetuximab in a panel of patient-derived tumour xenograft models: activation of MET as one mechanism for drug resistance.  _Eur J Cancer_ , 47(8):1231–1243, May 2011.
  * Sokal and Rohlf [1962] Robert R. Sokal and F. James Rohlf.  The comparison of dendrograms by objective methods.  _Taxon_ , 11(2):33–40, 1962.  ISSN 00400262.  URL http://www.jstor.org/stable/1217208.
  * Narayan et al. [2020] R. S. Narayan, P. Molenaar, J. Teng, F. M. G. Cornelissen, I. Roelofs, R. Menezes, R. Dik, T. Lagerweij, Y. Broersma, N. Petersen, J. A. Marin Soto, E. Brands, P. van Kuiken, M. C. Lecca, K. J. Lenos, S. G. J. G. In ’t Veld, W. van Wieringen, F. F. Lang, E. Sulman, R. Verhaak, B. G. Baumert, L. J. A. Stalpers, L. Vermeulen, C. Watts, D. Bailey, B. J. Slotman, R. Versteeg, D. Noske, P. Sminia, B. A. Tannous, T. Wurdinger, J. Koster, and B. A. Westerman.  A cancer drug atlas enables synergistic targeting of independent drug vulnerabilities.  _Nat Commun_ , 11(1):2935, 06 2020.
  * Neal [2003] Radford Neal.  Density modeling and clustering using Dirichlet diffusion trees.  _Bayesian Statistics_ , 7:619–629, 2003.
  * Gao et al. [2015] H. Gao, J. M. Korn, S. Ferretti, J. E. Monahan, Y. Wang, M. Singh, C. Zhang, C. Schnell, G. Yang, Y. Zhang, O. A. Balbin, S. Barbe, H. Cai, F. Casey, S. Chatterjee, D. Y. Chiang, S. Chuai, S. M. Cogan, S. D. Collins, E. Dammassa, N. Ebel, M. Embry, J. Green, A. Kauffmann, C. Kowal, R. J. Leary, J. Lehar, Y. Liang, A. Loo, E. Lorenzana, E. Robert McDonald, M. E. McLaughlin, J. Merkin, R. Meyer, T. L. Naylor, M. Patawaran, A. Reddy, C. Röelli, D. A. Ruddy, F. Salangsang, F. Santacroce, A. P. Singh, Y. Tang, W. Tinetto, S. Tobler, R. Velazquez, K. Venkatesan, F. Von Arx, H. Q. Wang, Z. Wang, M. Wiesmann, D. Wyss, F. Xu, H. Bitter, P. Atadja, E. Lees, F. Hofmann, E. Li, N. Keen, R. Cozens, M. R. Jensen, N. K. Pryer, J. A. Williams, and W. R. Sellers.  High-throughput screening using patient-derived tumor xenografts to predict clinical trial drug response.  _Nat. Med._ , 21(11):1318–1325, Nov 2015.
  * Konopleva et al. [2020] M. Konopleva, G. Martinelli, N. Daver, C. Papayannidis, A. Wei, B. Higgins, M. Ott, J. Mascarenhas, and M. Andreeff.  MDM2 inhibition: an important step forward in cancer therapy.  _Leukemia_ , 34(11):2858–2874, 11 2020.
  * Lapointe and Legendre [1991] François-Joseph Lapointe and Pierre Legendre.  The generation of random ultrametric matrices representing dendrograms.  _Journal of Classification_ , 8(2):177–200, Dec 1991.  ISSN 1432-1343.  doi:10.1007/BF02616238.  URL https://doi.org/10.1007/BF02616238.
  * McCullagh [2006] P. McCullagh.  Structured covariance matrices in multivariate regression models.  Technical report, Department of Statistics, University of Chicago, 2006\.
  * Bravo et al. [2009] H. C. Bravo, S. Wright, K. H. Eng, S. Keles, and G. Wahba.  Estimating tree-structured covariance matrices via mixed-integer programming.  _J Mach Learn Res_ , 5:41–48, 2009.
  * Knowles and Ghahramani [2015] D. A. Knowles and Z. Ghahramani.  Pitman-Yor diffusion trees for Bayesian hierarchical clustering.  _IEEE Trans Pattern Anal Mach Intell_ , 37(2):271–289, Feb 2015.
  * Turner et al. [2013] B. M. Turner, P. B. Sederberg, S. D. Brown, and M. Steyvers.  A method for efficiently sampling from distributions with correlated dimensions.  _Psychol Methods_ , 18(3):368–384, Sep 2013.
  * Knowles et al. [2011] David A. Knowles, Jurgen Van Gael, and Zoubin Ghahramani.  Message passing algorithms for dirichlet diffusion trees.  _International Conference on Machine Learning (ICML)_ , 2011.
  * Liu [1994] Jun S. Liu.  The collapsed Gibbs sampler in Bayesian computations with applications to a gene regulation problem.  _Journal of the American Statistical Association_ , 89(427):958–966, 1994.  ISSN 01621459.  URL http://www.jstor.org/stable/2290921.
  * Sisson et al. [2018] Scott A Sisson, Yanan Fan, and Mark Beaumont.  _Handbook of approximate Bayesian computation_.  CRC Press, 2018.
  * Casella and Berger [2001] George Casella and Roger Berger.  _Statistical Inference_.  Duxbury Resource Center, June 2001.  ISBN 0534243126.
  * Mathai [1980] A.M. Mathai.  Moments of the trace of a noncentral Wishart matrix.  _Communications in Statistics - Theory and Methods_ , 9(8):795–801, 1980.  doi:10.1080/03610928008827921.  URL https://doi.org/10.1080/03610928008827921.
  * Cardona et al. [2013] G. Cardona, A. Mir, F. Rosselló, L. Rotger, and D. Sánchez.  Cophenetic metrics for phylogenetic trees, after Sokal and Rohlf.  _BMC Bioinformatics_ , 14:3, Jan 2013.
  * Rashid et al. [2020] Naim U. Rashid, Daniel J. Luckett, Jingxiang Chen, Michael T. Lawson, Longshaokan Wang, Yunshu Zhang, Eric B. Laber, Yufeng Liu, Jen Jen Yeh, Donglin Zeng, and Michael R. Kosorok.  High-dimensional precision medicine from patient-derived xenografts.  _Journal of the American Statistical Association_ , 0(0):1–15, 2020.  doi:10.1080/01621459.2020.1828091.  URL https://doi.org/10.1080/01621459.2020.1828091.
  * Repetto et al. [2018] M. V. Repetto, M. J. Winters, A. Bush, W. Reiter, D. M. Hollenstein, G. Ammerer, P. M. Pryciak, and A. Colman-Lerner.  CDK and MAPK synergistically regulate signaling dynamics via a shared multi-site phosphorylation region on the scaffold protein Ste5.  _Mol Cell_ , 69(6):938–952, 03 2018.
  * Kurtzeborn et al. [2019] K. Kurtzeborn, H. N. Kwon, and S. Kuure.  MAPK/ERK Signaling in regulation of renal differentiation.  _Int J Mol Sci_ , 20(7), Apr 2019.
  * Balko et al. [2012] J. M. Balko, T. W. Miller, M. M. Morrison, K. Hutchinson, C. Young, C. Rinehart, V. Sánchez, D. Jee, K. Polyak, A. Prat, C. M. Perou, C. L. Arteaga, and R. S. Cook.  The receptor tyrosine kinase ErbB3 maintains the balance between luminal and basal breast epithelium.  _Proc Natl Acad Sci U S A_ , 109(1):221–226, Jan 2012.
  * Wu et al. [2021] L. Wu, M. R. Islam, J. Lee, H. Takase, S. Guo, A. M. Andrews, T. P. Buzhdygan, J. Mathew, W. Li, K. Arai, E. H. Lo, S. H. Ramirez, and J. Lok.  ErbB3 is a critical regulator of cytoskeletal dynamics in brain microvascular endothelial cells: implications for vascular remodeling and blood-brain-barrier modulation.  _J Cereb Blood Flow Metab_ , page 271678X20984976, Feb 2021.
  * Gunning et al. [2015] P. W. Gunning, U. Ghoshdastider, S. Whitaker, D. Popp, and R. C. Robinson.  The evolution of compositionally and functionally distinct actin filaments.  _J Cell Sci_ , 128(11):2009–2019, Jun 2015.
  * Haider et al. [2019] K. Haider, S. Rahaman, M. S. Yar, and A. Kamal.  Tubulin inhibitors as novel anticancer agents: an overview on patents (2013-2018).  _Expert Opin Ther Pat_ , 29(8):623–641, Aug 2019\.
  * Mahale et al. [2015] S. Mahale, S. B. Bharate, S. Manda, P. Joshi, P. R. Jenkins, R. A. Vishwakarma, and B. Chaudhuri.  Antitumour potential of BPT: a dual inhibitor of CDK4 and tubulin polymerization.  _Cell Death Dis_ , 6:e1743, May 2015.
  * Bohnacker et al. [2017] T. Bohnacker, A. E. Prota, F. Beaufils, J. E. Burke, A. Melone, A. J. Inglis, D. Rageot, A. M. Sele, V. Cmiljanovic, N. Cmiljanovic, K. Bargsten, A. Aher, A. Akhmanova, J. F. Díaz, D. Fabbro, M. Zvelebil, R. L. Williams, M. O. Steinmetz, and M. P. Wymann.  Deconvolution of Buparlisib’s mechanism of action defines specific PI3K and tubulin inhibitors for therapeutic intervention.  _Nat Commun_ , 8:14683, 03 2017.
  * Zhao et al. [2014] Y. Zhao, H. Yu, and W. Hu.  The regulation of MDM2 oncogene and its impact on human cancers.  _Acta Biochim Biophys Sin (Shanghai)_ , 46(3):180–189, Mar 2014.
  * Vora et al. [2014] S. R. Vora, D. Juric, N. Kim, M. Mino-Kenudson, T. Huynh, C. Costa, E. L. Lockerman, S. F. Pollack, M. Liu, X. Li, J. Lehar, M. Wiesmann, M. Wartmann, Y. Chen, Z. A. Cao, M. Pinzon-Ortiz, S. Kim, R. Schlegel, A. Huang, and J. A. Engelman.  CDK 4/6 inhibitors sensitize PIK3CA mutant breast cancer to PI3K inhibitors.  _Cancer Cell_ , 26(1):136–149, Jul 2014.
  * Bonelli et al. [2017] M. A. Bonelli, G. Digiacomo, C. Fumarola, R. Alfieri, F. Quaini, A. Falco, D. Madeddu, S. La Monica, D. Cretella, A. Ravelli, P. Ulivi, M. Tebaldi, D. Calistri, A. Delmonte, L. Ampollini, P. Carbognani, M. Tiseo, A. Cavazzoni, and P. G. Petronini.  Combined inhibition of CDK4/6 and PI3K/AKT/mTOR pathways induces a synergistic anti-tumor effect in malignant pleural mesothelioma cells.  _Neoplasia_ , 19(8):637–648, Aug 2017.
  * Yuan et al. [2019] Y. Yuan, W. Wen, S. E. Yost, Q. Xing, J. Yan, E. S. Han, J. Mortimer, and J. H. Yim.  Combination therapy with BYL719 and LEE011 is synergistic and causes a greater suppression of p-S6 in triple negative breast cancer.  _Sci Rep_ , 9(1):7509, 05 2019.
  * Geel et al. [2014] Robin Van Geel, Elena Elez, Johanna C. Bendell, Jason Edward Faris, Martijn P. J. K. Lolkema, Ferry Eskens, Anna Spreafico, Petr Kavan, Jean-Pierre Delord, Martin H. Schuler, Zev A. Wainberg, Yasuhide Yamada, Takayuki Yoshino, Tim Demuth, Emin Avsar, Arkendu Chatterjee, Peijuan Zhu, Rene Bernards, Josep Tabernero, and Jan HM Schellens.  Phase I study of the selective BRAFV600 inhibitor encorafenib (LGX818) combined with cetuximab and with or without the α-specific PI3K inhibitor BYL719 in patients with advanced BRAF-mutant colorectal cancer.  _Journal of Clinical Oncology_ , 32(15_suppl):3514–3514, 2014.  doi:10.1200/jco.2014.32.15_suppl.3514.  URL https://doi.org/10.1200/jco.2014.32.15_suppl.3514.
  * van Geel et al. [2017] R. M. J. M. van Geel, J. Tabernero, E. Elez, J. C. Bendell, A. Spreafico, M. Schuler, T. Yoshino, J. P. Delord, Y. Yamada, M. P. Lolkema, J. E. Faris, F. A. L. M. Eskens, S. Sharma, R. Yaeger, H. J. Lenz, Z. A. Wainberg, E. Avsar, A. Chatterjee, S. Jaeger, E. Tan, K. Maharry, T. Demuth, and J. H. M. Schellens.  A phase Ib dose-escalation study of encorafenib and cetuximab with or without alpelisib in metastatic BRAF-mutant colorectal cancer.  _Cancer Discov_ , 7(6):610–619, 06 2017.
  * Dummer et al. [2018a] R. Dummer, P. A. Ascierto, H. J. Gogas, A. Arance, M. Mandala, G. Liszkay, C. Garbe, D. Schadendorf, I. Krajsova, R. Gutzmer, V. Chiarion-Sileni, C. Dutriaux, J. W. B. de Groot, N. Yamazaki, C. Loquai, L. A. Moutouh-de Parseval, M. D. Pickard, V. Sandor, C. Robert, and K. T. Flaherty.  Encorafenib plus binimetinib versus vemurafenib or encorafenib in patients with BRAF-mutant melanoma (COLUMBUS): a multicentre, open-label, randomised phase 3 trial.  _Lancet Oncol_ , 19(5):603–615, 05 2018a.
  * Dummer et al. [2018b] R. Dummer, P. A. Ascierto, H. J. Gogas, A. Arance, M. Mandala, G. Liszkay, C. Garbe, D. Schadendorf, I. Krajsova, R. Gutzmer, V. Chiarion Sileni, C. Dutriaux, J. W. B. de Groot, N. Yamazaki, C. Loquai, L. A. Moutouh-de Parseval, M. D. Pickard, V. Sandor, C. Robert, and K. T. Flaherty.  Overall survival in patients with BRAF-mutant melanoma receiving encorafenib plus binimetinib versus vemurafenib or encorafenib (COLUMBUS): a multicentre, open-label, randomised, phase 3 trial.  _Lancet Oncol_ , 19(10):1315–1327, 10 2018b.
  * Robert et al. [2019] C. Robert, J. J. Grob, D. Stroyakovskiy, B. Karaszewska, A. Hauschild, E. Levchenko, V. Chiarion Sileni, J. Schachter, C. Garbe, I. Bondarenko, H. Gogas, M. Mandalá, J. B. A. G. Haanen, C. Lebbé, A. Mackiewicz, P. Rutkowski, P. D. Nathan, A. Ribas, M. A. Davies, K. T. Flaherty, P. Burgess, M. Tan, E. Gasal, M. Voi, D. Schadendorf, and G. V. Long.  Five-year outcomes with dabrafenib plus trametinib in metastatic melanoma.  _N Engl J Med_ , 381(7):626–636, 08 2019.
  * Heaukulani et al. [2014] Creighton Heaukulani, David A. Knowles, and Zoubin Ghahramani.  Beta diffusion trees.  In _Proceedings of the 31st International Conference on International Conference on Machine Learning - Volume 32_ , ICML’14, page II–1809–II–1817. JMLR.org, 2014.