# Signalling entropy: a novel network-theoretical framework for systems
analysis and interpretation of functional omic data

Andrew E. Teschendorff  Peter Sollich  Reimer Kuehn  CAS-MPG Partner Institute
for Computational Biology, Chinese Academy of Sciences, Shanghai Institute for
Biological Sciences, 320 Yue Yang Road, Shanghai 200031, China.  Statistical
Cancer Genomics, Paul O’Gorman Building, UCL Cancer Institute, University
College London, London WC1E 6BT, UK.  Department of Mathematics, King’s
College London, London WC2R 2LS, UK.

###### Abstract

A key challenge in systems biology is the elucidation of the underlying
principles, or fundamental laws, which determine the cellular phenotype.
Understanding how these fundamental principles are altered in diseases like
cancer is important for translating basic scientific knowledge into clinical
advances. While significant progress is being made, with the identification of
novel drug targets and treatments by means of systems biological methods, our
fundamental systems level understanding of why certain treatments succeed and
others fail is still lacking. We here advocate a novel methodological
framework for systems analysis and interpretation of molecular omic data,
which is based on statistical mechanical principles. Specifically, we propose
the notion of cellular signalling entropy (or uncertainty), as a novel means
of analysing and interpreting omic data, and more fundamentally, as a means of
elucidating systems-level principles underlying basic biology and disease. We
describe the power of signalling entropy to discriminate cells according to
differentiation potential and cancer status. We further argue the case for an
empirical cellular entropy-robustness correlation theorem and demonstrate its
existence in cancer cell line drug sensitivity data. Specifically, we find
that high signalling entropy correlates with drug resistance and further
describe how entropy could be used to identify the achilles heels of cancer
cells. In summary, signalling entropy is a deep and powerful concept, based on
rigorous statistical mechanical principles, which, with improved data quality
and coverage, will allow a much deeper understanding of the systems biological
principles underlying normal and disease physiology.

###### keywords:

entropy; network; signalling; genomics; drug resistance; cancer;
differentiation; stem cell

††journal: Methods

##  1 Introduction

Recent advances in biotechnology are allowing us to measure cellular
properties at an unprecedented level of detail [1]. For instance, it is now
possible to routinely measure various molecular entities (e.g. DNA
methylation, mRNA and protein expression, SNPs) genome-wide in hundreds if not
thousands of cellular specimens [2]. In addition, other molecular data
detailing interactions between proteins or between transcription factors and
regulatory DNA elements are growing at a rapid pace [1]. All these types of
data are now widely referred to collectively as “omic” data. The complexity
and high-dimensional nature of this omic data presents a daunting challenge to
those wishing to analyse and interpret the data [1]. The difficulty of
analysing omic data is further compounded by the inherent complexity of
cellular systems. Cells are prime examples of organized complex systems,
capable of highly stable and predictable behaviour, yet an understanding of
how this deterministic behaviour emerges from what is a highly complex and
probabilistic pattern of dynamic interactions between numerous intra and
extracellular components, still eludes us [1]. Thus, elucidating the systems-
biological laws or principles dictating cellular phenotypes is also key for an
improved analysis and interpretation of omic data. Furthermore, important
biological phenomena such as cellular differentiation are fundamentally
altered in diseases like cancer [3]. Hence, an attempt to understand how
cellular properties emerge at a systems level from the properties seen at the
individual gene level is not only a key endeavour for the systems biology
community, but also for those wanting to translate basic insights into
effective medical advances [4, 5].
It is now well accepted that at a fundamental level most biological systems
are best modeled in terms of spatial interactions between specific entities
(e.g. neurons in the case of the brain), which may or may not be dynamically
changing in time [6, 7]. It therefore seems natural to also use the
mathematical and physical framework of networks to help us analyse and
interpret omic data at a systems level [8, 9]. Indeed, the cellular phenotype
is determined to a large extent by the precise pattern of molecular
interactions taking place in the cell, i.e. a molecular interaction network
[10]. Although this network has spatial and dynamic dimensions which at
present remain largely unexplored due to technological or logistical
limitations, there is already a growing number of examples where network-based
analysis strategies have been instrumental [11, 12, 13]. For instance, a
deeper understanding of why sustained treatment with EGFR inhibitors can lead
to dramatic sensitization of cancer cell lines to cytotoxic agents was
possible thanks to a systems approach [13]. Another study used reverse
engineering network approaches to identify and validate drug targets in
glioblastoma multiforme, to be further tested in clinical trials [14]. What is
key to appreciate here is that these successes have been achieved in spite of
noisy and incomplete data, suggesting that there are simple yet deep systems
biological principles underlying cellular biology that we can already probe
and exploit with current technology and data. Thus, with future improvements
in data quality and coverage, network-based analysis frameworks will play an
ever increasing and important role in systems biology, specially at the level
of systems analysis and interpretation [4]. Therefore, it is also imperative
to develop novel, more powerful, network-theoretical methods for the systems
analysis of omic data.
In adopting a network’s perspective for the analysis and interpretation of
omic data, there are in principle two different (but not mutually exclusive)
approaches one can take. One possibility is to infer (i.e. reverse engineer)
the networks from genome-wide data [15]. Most of these applications have done
this in the context of gene expression data, with the earliest approaches
using clustering or co-expression to postulate gene interdependencies [16].
Partial correlations and Graphical Gaussian Models have proved useful as a
means of further refining correlation networks by allowing one to infer the
more likely direct interactions while simultaneously also filtering out those
which are more likely to be indirect [17]. These methods remain popular and
continue to be studied and improved upon [18, 19]. Other methods have drawn on
advanced concepts from information theory, for instance ARACNe (“Algorithm for
the Reconstruction of Accurate Cellular Networks”) has been shown to be
successful in infering regulatory networks in B-cells [15].
In stark contrast to reverse engineering methods, another class of algorithms
have used structural biological networks from the outset, using these as
scaffolds to integrate with omic data. Specifically, by using a structural
network one can sparsify the correlation networks inferred from reverse-
engineering approaches, thus providing another means of filtering out
correlations that are more likely to be indirect [9]. Besides, integration
with a structural network automatically provides an improved framework for
biological interpretation [10, 20, 21, 22]. The structural networks themselves
are typically derived from large databases, which detail literature curated
experimentally verified interactions, including interactions derived from
Yeast 2 Hybrid screens (Y2H) [23]. The main example is that of protein protein
interaction (PPI) maps, which have been generated using a number of different
complementary experimental and in-silico approaches, and merging these maps
from these different sources together has been shown to be an effective
approach in generating more comprehensive high-confidence interaction networks
[24, 25]. PPI networks have been used mainly as a means of integrating and
analysing gene expression data (see e.g. [20, 26, 27]). More recently, this
approach has also been successfully applied in the DNA methylation context,
for instance it has been shown that epigenetic changes associated with age
often target specific gene modules and signalling pathways [28].
Another class of methods that have used structural networks, PPIs in
particular, have integrated them with gene expression data to define an
approximation to the underlying signaling dynamics on the network, thus
allowing more in-depth exploration of the interplay between network topology
and gene expression. Typically, these studies have used the notion of random
walks on weighted graphs where the weights are constructed from differential
expression statistics, and where the aim is to identify nodes (genes) in the
network which may be important in dictating the signaling flows within the
pathological state. For instance, among these random walk methods is NetRank,
a modification of the Google PageRank algorithm, which was able to identify
novel, robust, network based biomarkers for survival time in various cancers
[29, 30]. Other random walk based approaches, aimed at identifying causal
drivers of specific phenotypes (e.g. expression or cancer), have modeled
signal transduction between genes in the causal and phenotypic layers as flows
in an electric circuit diagram, an elegant formulation capable of not only
identifying the likely causal genes but also of tracing the key pathways of
information flow or dysregulation [31, 32]. Random walk theory has also been
employed in the development of differential network methodologies. An example
is NetWalk [33], which is similar to NetRank but allows differential signaling
fluxes to be inferred. This approach was successful in identifying and
validating the glucose metabolic pathway as a key determinant of lapatinib
resistance in ERBB2 positive breast cancer patients [34].
Another important concept to have emerged recently is that of network rewiring
[35, 36, 37]. This refers to the changes in interaction patterns that
accompany changes in the cellular phenotype. Network rewiring embodies the
concept that it is the changes in the interaction patterns, and not just the
changes in absolute gene expression or protein activity, that are the main
determinants of the cellular phenotype. That network rewiring may be key to
understanding cellular phenotypes was most convincingly demonstrated in a
differential epistasis mapping study conducted in yeast cells exposed to a DNA
damaging agent [35]. Specifically, what this study demonstrated is that
responses to perturbations or cellular stresses are best understood in terms
of the specific rewiring of protein complexes and functional modules. Thus,
this conceptual framework of network rewiring may apply equally well to the
genetic perturbations and cellular stresses underlying disease pathologies
like cancer.
In this article we advocate a network-theoretical framework based on
statistical mechanical principles and more specifically on the notion of
signalling entropy [8, 9]. This theoretical framework integrates gene
expression (but in principle also other functional data) with a PPI network,
merging existing concepts such as signaling dynamics (i.e. random walks) and
network rewiring with that of signalling entropy. In previous work we have
shown how signalling entropy (i) provides a proxy to the elevation in
Waddington’s epigenetic landscape, correlating with a cell’s differentiation
potential [38], (ii) how it can be used to identify signaling pathways and
nodes important in differentiation and cancer [38, 8, 9], and (iii) how it
predicts two cancer system-omic hallmarks: (a) cancer is characterised by an
increase in signalling entropy and (b) local signaling entropy changes anti-
correlate with differential gene expression [8, 9]. Here, we present and unify
the different signaling entropy measures used previously and further explore a
novel application of signalling entropy to understanding drug sensitivity
profiles in cancer cell lines. Specifically, we first use simulated data to
justify the existence of an entropy-robustness theorem, and subsequently
provide empirical evidence for this theorem by demonstrating that increases in
local signalling entropy correlate with drug resistance (robustness). We
further show the importance of network topology in dictating the signalling
entropy changes underlying drug response. In addition, we provide R-functions
implementing the entropy rate calculation and ranking of genes according to
differential entropy, all freely available from
sourceforge.net/projects/signalentropy/files/ .

##  2 Materials and Methods

###  2.1 Basic rationale and motivation for Signalling Entropy: understanding
systems biology through uncertainty

Loosely defined, entropy of a system, refers to a measure of the disorder,
randomness or uncertainty of processes underlying the system’s state. In the
context of a single cell, signalling entropy will reflect the amount of
overall disorder, randomness or uncertainty in how information, i.e.
signaling, is passed on in the molecular interaction network. At a fundamental
level, all signaling in a cell is probabilistic, determined in part by the
relative cellular concentrations of the interacting molecules. Hence, in
discussing cellular signalling entropy, it is useful to picture a molecular
interaction network in which edges represent possible interactions and with
the edge weights reflecting the relative probabilities of interaction
(Fig.1A). Thus, an interaction network in a high-entropy state is
characterised by signaling interaction probabilities that are all fairly
similar in value, whereas a low-entropy state will be characterised by
specific signalling interactions possessing much higher weights (Fig.1A).

Figure 1: Signalling Entropy: understanding systems biology through
uncertainty. A) A caricature model of a cellular interaction network with edge
widths/color indicating the relative probabilities of interaction. On the left
and right, we depict states of high and low signalling entropy, respectively.
At the cellular population level, this translates into samples of high and low
intra-sample heterogeneity, respectively. B) Signalling entropy correlates
with pluripotency as demonstrated in our previous work (Banerji et al 2013).
The pluripotent state is a highly promiscuous signaling state, generating high
intra-sample heterogeneity, and allowing the stem cell population to
differentiate into any cell type. In contrast, in a terminally differentiated
state, signaling is almost deterministic, reflecting activation of very
specific pathways in the majority of cells, leading to a highly homogeneous
and differentiated cell population. Thus, signalling entropy defines the
height in Waddington’s differentiation landscape.

Why would this type of signalling entropy, loosely defined as the amount of
uncertainty in the signaling interaction patterns, be useful to systems
biology? One way to think of signalling entropy is as representing signaling
promiscuity, which has been proposed as a key systems feature underlying the
pluripotent or multipotent capacity of cells (Fig.1B) [39, 40, 41, 42].
Indeed, it has been suggested that pluripotency is an intrinsic statistical
mechanical property, best defined at the cellular population level [40].
Specifically, it has been demonstrated that pluripotent stem cells exhibit
remarkable heterogeneity in gene expression levels, including well known stem
cell markers such as NANOG [40]. It is also well known that a large number of
genes, many encoding transcription factors, exhibit low-levels of expression
in stem cells, yet simultaneously are being kept in a poised chromatin state,
allowing immediate activation if this were required [43]. Thus, in a
pluripotent stem cell like state, signal transduction is in a highly
egalitarian and, thus, promiscuous state, i.e. a state of high signalling
entropy. Conversely, differentiation leads, by necessity, to activation of
specific transcription factors and pathways and thus to a lowering in the
uncertainty of signaling patterns, and thus to a lowering of entropy. We
recently demonstrated, using gene expression profiles of over 800 samples,
comprising cells of all major stages of differentiation, including human
embryonic stem cells (hESCs), induced pluripotent stem cells (iPSCs),
multipotent cell types (e.g. hematopoietic stem cells (HSCs)), and terminally
differentiated cells within these respective lineages, that signalling entropy
not only correlates with differentiation potential but that it provides a
highly quantitative measure of potency [38]. Indeed, we showed that signalling
entropy provides a reasonably good approximation to the energy potential in
Waddington’s epigenetic landscape [38].
Here we decided to explore the concept of signalling entropy in relation to
cellular robustness and specifically to drug resistance in cancer. That
signalling entropy may be informative of a cell’s robustness is a proposal
that stems from a general (but unproven) theorem, first proposed by Manke and
Demetrius [44, 45, 46]: namely, that a system’s entropy and robustness are
correlated. Mathematically, this can be expressed as $\Delta S\Delta R>0$,
which states that a positive change in a system’s entropy (i.e. $\Delta S>0$)
must lead to an increase in robustness ($\Delta R>0$). Now, cells are fairly
robust entities, having evolved the capacity to buffer the intra-and-
extracellular stresses and noise which they are constantly exposed to [47,
39]. Much of this overall stability and robustness likely stems from the
topological features of the underlying signaling and regulatory networks, for
instance features such as scale-freeness and hierarchical modularity, which
are thought to have emerged through natural evolutionary processes such as
gene duplication [10, 48, 49]. However, another key feature which contributes
to cellular robustness is cross-talk and signalling pathway redundancy [47].
Pathway redundancy refers to a situation where a cell has the choice of
transmitting signals via two or more possible routes. In the language of
statistical mechanics, this corresponds to a state of high uncertainty (or
entropy) in signaling. High signalling entropy could thus underpin a cell’s
robustness to perturbations, suggesting that a cell’s entropy and robustness
may indeed be correlated (Fig.2A). Consistent with this, pathway redundancy is
also well recognized to be a key feature underlying drug resistance of cancer
cells [50].

Figure 2: Signalling entropy and cellular robustness: A) Signalling entropy
ought to correlate with cellular robustness. The inequality encapsulates this
information by stating that a decrease in signalling entropy (i.e. if $\Delta
S<0$), then the system’s robustness $R$ must also decrease, i.e. $\Delta R<0$,
so that the product $\Delta S\Delta R>0$. Observe how in the low entropy
state, random removal of edges through e.g. inactivating mutations, can lead
to deactivation of a key signaling pathway connecting a given start and end
nodes (shown in orange). In the high entropy state, the same perturbations do
not prevent signal transduction between the two orange nodes. B) Depicted are
the effects of two major forms of cancer perturbation. In the upper panel,
inactivation (typically of tumour suppressors), leads to underexpression and a
corresponding loss of correlations/interactions with neighbors in the PPI
network. This is tantamount to a state of increased entropy and drug
intervention is unlikely to be effective. In the lower panel, we depict the
case of an oncogene, which is overexpressed in cancer. This overexpression
leads to activation of a specific oncogenic pathway which results in oncogene
addiction and increased sensitivity to targeted drug intervention. Thus, local
signalling entropy and robustness (as determined by response to a drug), may
also correlate locally.

Further supporting the notion that entropy and robustness may be correlated,
we previously showed that (i) cancer is characterised by a global increase in
signalling entropy compared to its respective normal tissue [9], in line with
the observation that cancer cells are specially robust to generic
perturbations, and (ii) that a gene’s differential entropy between normal and
cancer tissue is generally speaking anticorrelated with its differential
expression [9], consistent with the view that cancer cells are specially
sensitive to drug interventions that target overexpressed oncogenes, a
phenomenon known as oncogene addiction (Fig.2B) [3]. Interestingly, the
observations that differential entropy and differential expression are anti-
correlated and that cancer is characterised globally by an increase in entropy
[9], are also consistent with the prevalent view that most driver mutations
are indeed inactivating, targeting tumor suppressor genes [51, 52, 53]. Hence,
based on all of these observations and insights, we posited that signalling
entropy could also prove useful as a means of predicting drug resistance in
cancer cells.

###  2.2 The Signalling Entropy Method: definitions and construction

Briefly, we review the definitions and construction of signalling entropy as
used in our previous studies [8, 9, 38]. The construction relies on a
comprehensive and high-confidence PPI network which is integrated with gene
expression data [24, 9] (see Appendix A). Briefly, the PPI is used as a
scaffold, and edge weights are constructed from the gene expression data to
approximate the interaction or signaling probabilities between the
corresponding proteins in the PPI. Thus, the gene expression data provides the
biological context in which to modulate the PPI interaction probabilities. To
compute signalling entropy requires the estimation of a stochastic matrix,
reflecting these interaction probablities over the network.
The construction of the stochastic matrix can proceed in two different ways.
In the earliest studies we used a construction which was done at the level of
phenotypes [8, 9]. Under this model, edge weights $w_{ij}$ between proteins
$i$ and $j$ were constructed from the correlation between the expression
levels of the corresponding genes $i$ and $j$, as assessed over independent
samples all representing the same phenotype. Estimating the correlations over
independent samples, all within the same phenotype, can be viewed as
representing a multifactorial perturbation experiment, with e.g. genetic
differences between individuals mimicking specific perturbations, and thus
allowing influences to be inferred. Thus, this approach hinges on the
assumption that highly correlated genes, whose coding proteins interact, are
more likely to be interacting in the given phenotype than two genes which do
not correlate. The use of a PPI as a scaffold is important to filter out
significant correlations which only result from indirect influences. The
correlations themselves can be defined in many different ways, for instance,
using simple Pearson correlations or non-linear measures such as Mutual
Information. For example, one way to define the weights is as
$w_{ij}=\frac{1}{2}(1+c_{ij})$ with $c_{ij}$ describing the Pearson
correlation coefficient between genes $i$ and $j$. This definition guarantees
positivity, but also treats positive and negative correlations differently,
which makes biological sense because activating and inhibitory interactions
normally have completely distinct consequences on downstream signalling. Thus,
the above weight scheme treats zero or insignificant correlations as
intermediate, consistent with the view that an absent interaction is neither
inhibitory nor activating. However, other choices of weights are possible:
e.g. $w_{ij}=|c_{ij}|$, which treats negative and positive correlations on an
equal footing. Once edge weights are defined as above, these are then
normalised to define the stochastic matrix $p_{ij}$ over the network,

$p_{ij}=\frac{w_{ij}}{\sum_{k\in\mathcal{N}_{i}}w_{ik}},$

with $\mathcal{N}_{i}$ denoting the PPI neighbors of gene $i$. Thus, $p_{ij}$
is the probability of interaction between genes $i$ and $j$, and as required,
$\sum_{j}{p_{ij}}=1$. However, there is no requirement for $p_{ij}$ to be
doubly stochastic, i.e. $P$ is in general not a symmetric matrix. Hence, edges
are bi-directional with different weights labeling the probability of signal
transduction from $i$ to $j$, and that from $j$ to $i$ ($p_{ij}\neq p_{ji}$).
An alternative to the above construction of the stochastic matrix is to invoke
the mass action principle, i.e. one now assumes that the probability of
interaction in a given sample is proportional to the product of expression
values of the corresponding genes in that sample [38]. Thus, the PPI is again
used as a scaffold to only allow interactions supported by the PPI network,
but the weights are defined using the mass action principle, as

$w_{ij}\propto E_{i}E_{j}$

where $E_{i}$ denotes the normalised expression intensity value (or normalised
RNA-Seq read count) of gene $i$. An important advantage of this construction
is that the stochastic matrix is now sample specific, as the expression values
are unique to each sample.
Given a stochastic matrix, $p_{ij}$, constructed using one of the two methods
above, one can now define a local Shannon entropy for each gene $i$ as

$\tilde{S}_{i}=-\frac{1}{\log{k_{i}}}\sum_{k\in\mathcal{N}_{i}}p_{ik}\log{p_{ik}},$

where $k_{i}$ denotes the degree of gene $i$ in the PPI network. The
normalisation is optional but ensures that this local Shannon entropy is
normalised between 0 and 1. Clearly, if only one weight is non-zero, then the
entropy attains its minimal value (0), representing a state of determinism or
lowest uncertainty. Conversely, if all edges emanating from $i$ carry the same
weight, the entropy is maximal (1), representing the case of highly
promiscuous signaling. In principle, local Shannon entropies can thus be
compared between phenotypes to identify genes where there are changes in the
uncertainty of signaling. In the case where entropies are estimated at the
phenotype level, jackknife approaches can be used to derive corresponding
differential entropy statistics [9]. Deriving statistics is important because
node degree has a dramatic influence on the entropy variance, with high degree
nodes exhibiting significantly lower variability in absolute terms, which
nevertheless could be highly significant [9]. In the case where entropies are
estimated at the level of individual samples, ordinary statistical tests (e.g.
rank sum tests) can be used to derive sensible P-values, assuming of course
that enough samples exist within the phenotypes being compared.
In addition to the local entropy, it is also of interest to consider
statistical properties of the distribution of local entropies, for instance
their average. Comparing the average of local entropies between two phenotypes
would correspond to a comparison of non-equilibrium entropy rates. To see
this, consider the formal definition of the entropy rate $SR$ [54, 55], i.e.

$SR=\sum_{i=1}^{n}\pi_{i}S_{i},$

where $\pi_{i}$ is the stationary distribution (equivalently the left
eigenvector with unit eigenvalue) of $P$ (i.e. $\pi P=\pi$), and where now

$S_{i}=-\sum_{k\in\mathcal{N}_{i}}p_{ik}\log{p_{ik}}.$

Note that the entropy rate $SR$ is an equilibrium entropy since it involves
the stationary distribution of the random walker. As such, the entropy rate
also depends on the global topology of the network. Thus, the entropy rate is
a weighted average of the local unnormalized entropies, with the weights
specified by the stationary distribution. It follows that comparing the
unweighted averages of local entropies reflects a comparison of a non-
equilibrium entropy rate since the stationary distribution is never used.

Figure 3: Entropy rate in normal and cancer tissue: A) Boxplots of sample
specific entropy rates comparing normal liver and liver cancer samples.
Expression data set is the one used in [38]. B) As A), but also shown are the
sample specific entropy rates obtained by randomly permuting the gene
expression values over the network. Note how the entropy rates for the normal
and cancer states are significantly reduced upon permutation and are no longer
highly discriminative between normal and cancer. C) ROC curves and associated
AUC normal-cancer discriminatory values for the unpermuted and permutated
cases depicted in A) and B).

###  2.3 The importance of the integrated weighted network in estimating
signalling entropy

The entropy rate constructed using the mass action principle is sample
specific. We previously demonstrated that this entropy rate was highly
discriminative of the differentiation potential of samples within a
developmental lineage, as well as being highly discriminative of normal and
cancer tissue [38] (Fig.3A). Since the entropy rate takes as input a PPI
network and a sample’s genome-wide expression profile, the significance of the
resulting entropy values, as well as that of the difference between entropy
rates, also needs to be assessed relative to a null distribution in which the
putative information content between network and gene expression values is
non-existent. In fact, since the weights in the network determine the entropy
rate and these weights are dependent on both the specific network nodes and
their respective gene expression profiles, it is natural to assess the
significance of the entropy rate by “destroying” the mutual information
between the network nodes and their gene expression profiles, for instance by
randomising (i.e. permuting) the gene expression profiles over the network.
Thus, under this randomisation, the topological properties of the network
remain fixed, but the weights are redefined. Application of this randomisation
procedure to the normal/cancer expression set considered previously [38]
(Fig.3A) shows that the discriminatory potential is significantly reduced upon
permuting the gene expression values over the network (Fig.3B-C). Importantly,
we observe that the entropy rate is much higher in the normal and cancer
states compared to the rates obtained upon randomisation of the gene
expression profiles (Fig.3B), indicating that both normal and cancer states
are characterised by a higher level of signaling promiscuity compared to a
network with random weights. That the discrimination between normal and cancer
is significantly reduced in the randomly weighted network further demonstrates
that there is substantial mutual information between the PPI network and the
gene expression profiles, thus justifying the signalling entropy approach.

###  2.4 Signalling entropy R-package: input, output and code availability

A vignette/manual and user-friendly R-scripts that allow computation of the
entropy rate is available at the following url:
sourceforge.net/projects/signalentropy/files/. Here we briefly describe the
salient aspects of this package:
The input: The main R-script provided (CompSR) takes as input a user-specified
PPI network, and a genome-wide expression vector representing the gene
expression profile of a sample. It is assumed that this has been generated
using either Affy, Illumina or RNA-Sequencing. In principle one ought to use
as gene expression value the normalised unlogged intensity (Affy/Illu) or RNA-
seq count, since this is what should be proportional to the number of RNA
transcripts in the sample. However, in practice we advise taking the log-
transformed normalised value since the log-transformation provides a better
compromise between proportionality and regularisation, i.e. some
regularisation is advisable since the underlying kinetic reactions are also
regular.
The output: The R-functions provided in the package then allow the user to
estimate the global entropy rate for a given sample, as well as the local
normalised entropies for each gene in the integrated network. If a phenotype
is specified then genes can be ranked according to the correlation strength of
their local entropy to the phenotype. Thus, the signalling entropy method
allows us to assess (i) if the overall levels of signalling promiscuity is
different between phenotypes, which could be important, for instance, to
compare the pluripotent capacity of iPSCs generated via different protocols or
to assess the stem cell nature of putative cancer stem cells [38], and (ii) to
rank genes according to differential entropy between phenotypes, allowing key
signalling genes associated with differentiation, metastasis or cancer to be
identified [8, 9, 38].

##  3 Results

###  3.1 Signalling entropy and cellular robustness

Our previous observation that signalling entropy is increased in cancer [9],
and that one of the key characteristics of cancer cells is their robustness to
intervention and environmental stresses, suggested to us that high cellular
signalling entropy may be a defining feature of a cell that is robust to
general perturbations. In addition, cancer cells often exhibit the phenomenon
of oncogene addiction, whereby they become overly reliant on the activation of
a specific oncogenic pathway, rendering them less robust to targeted
intervention [3]. Since oncogenes are normally overexpressed in cancer, it
follows by the second cancer system-omic hallmark [9], that their lower
signalling entropy may underpin their increased sensitivity to targeted drugs
(Fig.2). Based on these insights, we posited that cellular signalling entropy
may be correlated to the overall cellular system’s robustness to
perturbations.
In order to explore the relation between signalling entropy and robustness in
a general context, it is useful to consider another network property of the
stochastic matrix, namely the global mixing rate. This mixing rate is defined
formally as the inverse of the number of timesteps for which one has to evolve
a random walk on the graph so that its probability distribution is close to
the stationary distribution, independently of the starting node (Appendix B).
This is a notion that is closer to that of robustness or resilience as defined
by Demetrius and Manke [44, 45, 46], allowing the mixing rate to be viewed as
a crude proxy of a system’s overall robustness to generic perturbations.
Furthermore, the global mixing rate can be easily estimated as

$\mu_{R}=-\log{\mathrm{SLEM}}$ (1)

where $SLEM$ is the second largest (right) eigenvalue modulus of the
stochastic matrix (Appendix B). Thus, we first asked how this global mixing
rate is related to the signalling entropy rate.

Figure 4: Correlation between global entropy and mixing rates. A) Plotted is
the entropy rate against the mixing rate for Erdos-Renyi graphs of 100 nodes
and average degree 10. The light blue lines show the results over 100
different network realisations, with the dark blue line representing the
ensemble average. B) As A) but for connected subnetworks of average size 500
nodes, generated by random subsampling of 1000 nodes from the full PPI network
of 8038 nodes. As in A), a range of edge weight distributions were considered
reflecting variable entropy rates. The light blue lines show the results over
100 different realisations, with the dark blue line representing the ensemble
average. C) Scatterplot of the global entropy versus mixing rates for the 488
Cancer Cell-Line Encyclopedia (CCLE) samples. Spearman rank correlation
coefficient (SCC) and associated P-value are given (left panel). Wilcoxon rank
sum test P-value between high and low mixing rate groups (as defined by
tertiles) (right panel).

For a regular network of degree $d$ it can be easily shown that the entropy
rate $SR=\log{d}$, whilst results on graph theory also show that for
sufficiently large regular graphs, $\mu_{R}\propto\log{d}$ [56]. Hence, at
least for regular networks a direct correlation exists. It follows that for
non-regular graphs with tight degree distributions, e.g. Erdös-Renyi (ER)
graphs, the entropy and mixing rates should also be approximately correlated.
Indeed, using an empirical computational approach to evaluate the entropy and
mixing rates for ER graphs with variable entropy rates, we were able to
confirm this correlation (Appendix C, Fig.4A). Next, we wanted to investigate
if this relationship also holds for networks with more realistic topologies
than ER graphs. Hence, we generated connected networks on the order of 500
nodes by random subsampling 1000 nodes from our large protein interaction
network ($\sim 8000$ nodes) followed by extraction of the maximally connected
component (Appendix D). We verified that these networks possessed approximate
scale-free topologies with clustering coefficients which were also
significantly higher than for ER graphs. As before, for each generated
network, stochastic matrices of variable entropy rates were defined.
Signalling entropy and mixing rates were then estimated for each of these
networks, and subsequently averaged over an ensemble of such graphs. As with
the random Poisson (ER) graphs, average mixing and entropy rates were highly
correlated (Fig.4B).
Having demonstrated a direct correlation between these two measures on
simulated data, we next asked if such a correlation could also be present in
the full protein interaction networks and with the stochastic matrices derived
from real expression data. Thus, we computed the global entropy and mixing
rates for 488 cancer cell lines from the Cancer Cell Line Encylopedia (CCLE)
(Appendix E) [57]. Remarkably, despite the dynamic ranges of both entropy and
mixing rates being much smaller (Fig.4C) compared to those of the artificially
induced stochastic matrices (c.f Fig.4A-B), we were still able to detect a
statistically significant correlation between entropy and mixing rates, as
estimated across the 488 cell lines (Fig.4C). Thus, all these results support
the view that global entropy and mixing rates are correlated, albeit only in
an average/ensemble sense.

###  3.2 Local signalling entropy predicts drug sensitivity

In the case of realistic expression and PPI data, the observed correlation
between entropy and mixing rates was statistically significant but weak
(Fig.4C). This could be easily attributed to the fact that in real biological
networks, the global mixing rate is a very poor measure of cellular
robustness. In fact, it is well known that cellular robustness is highly
sensitive to which genes undergo the perturbation. For instance, in mice some
genes are embryonically lethal, whereas others are not [10]. Robustness itself
also admits many different definitions. Because of this, we decided to
investigate signalling entropy in relation to other more objective measures of
cellular robustness. One such measure is drug sensitivity, for instance, IC50
values, which measure the amount of drug dose required to inhibit cell
proliferation by 50$\%$. According to this measure, a cell that is insensitive
to drug treatment is highly robust to that particular treatment. Since most
drugs target specific genes, we decided to explore the relation, if any,
between the local signalling entropy of drug targets and their associated drug
sensitivity measures. Specifically, we hypothesized that since local entropy
provides a proxy for local pathway redundancy, that it would correlate with
drug resistance. To test this, we first computed for each of the 8038 genes in
the PPI network its local signalling entropy in each of the 488 CCLE cancer
cell-lines. To identify associations between the 24 drug sensitivity profiles
and the 8038 local network entropies, we computed non-linear rank (Spearman)
correlation coefficients across the 488 cancer cell-lines, resulting in
$24\times 8038$ correlations and associated P-values. We observed that there
were many more significant associations than could be accounted for by random
chance (Fig.5A), with the overall strength of association very similar to that
seen between gene expression and drug sensitivity (Fig.5B).

Figure 5: Anti-correlation between local signalling entropy and drug
sensitivity. A) Histogram of Spearman rank correlation P-values between drug
sensitivities (n=24) and local signalling entropies (n=8038 genes), as
computed over the 488 CCLE cell-lines. B) As A) but for gene expression
instead of signalling entropy. C) Scatterplot of Spearman rank Correlation
Coefficient (SCC) between local signalling entropy (sigS) and drug sensitivity
(DS) against $-log_{10}$P-value for each of 134 drug gene target pairs. D) As
C) but for gene expression instead of local entropy. In C) & D), we give the
distribution of significant positive and negative correlations and the
associated Binomial test P-value. E) Drug target gene pairs ranked according
to negative SCC (cyan color) between signalling entropy and drug sensitivity.
Only pairs where at least one of entropy or gene expression were significantly
associated are shown. Upper panels show the SCC (cyan=strong negative SCC,
white=zero or non-significant SCC, magenta=strong positive SCC), while lower
panels show the corresponding P-values with the last row showing the P-value
from the Likelihood Ratio Test (LRT) assessing the added predictive value of
signalling entropy over gene expression. The darker the tones of red the more
significant the P-values are, whilst white indicates non-significance. F) A
subset of E), with pairs now ranked according to the LRT P-value.

One would expect the targets of specific drugs to be highly informative of the
sensitivity of the corresponding drugs. We used the CancerResource [58] to
identify targets of the 24 drugs considered here, and found a total of 154
drug-target pairs. For 134 of these pairs we could compute a P-value of
association between local entropy and drug sensitivity with 76 pairs (i.e.
57$\%$) exhibiting a significant association (Fig.5C). This was similar to the
proportion (54$\%$) of observed significant associations between gene
expression and drug sensitivity (Fig.5D). However, interestingly, only 42 of
these significant drug-target pairs were in common between the 76 obtained
using signalling entropy and the 72 obtained using gene expression.
Importantly, while the significant associations between gene expression and
drug sensitivity involved preferentially positive correlations, in the case of
signalling entropy most of the significant correlations were negative
(Fig.5C-D), exactly in line with our hypothesis that high entropy implies drug
resistance. Thus, as expected, cell-lines with highly expressed drug targets
were more sensitive to treatment by the corresponding drug, but equally, drug
targets exhibiting a relatively low signalling entropy were also predictive of
drug sensitivity.
To formally demonstrate that local signalling entropy adds predictive power
over gene expression, we considered bi-variate regression models including
both the target’s gene expression as well as the target’s signalling entropy.
Using such bivariate models and likelihood ratio tests we found that in the
majority of cases where signalling entropy was significantly associated with
drug sensitivity that it did so independently of gene expression, adding
predictive value (Fig.5E). Top ranked drug-target pairs where signalling
entropy added most predictive value over gene expression included
Topotecan/TP53 and Paclitaxel/MYC (Fig.5F).
To further demonstrate that the observed associations between local signalling
entropy and drug sensitivity are statistically meaningful, we conducted a
control test, in which we replaced in the multivariate model the signalling
entropy of the target with a non-local entropy computed by randomly replacing
the PPI neighbours of the target with other “far-away” genes in the network.
For a considerable fraction (41$\%$) of drug-target pairs, the original
multivariate models including the local entropy constituted better predictive
models than those with the non-local entropy (false discovery rate $<0.25$),
indicating that the observed associations are driven, at least partly, by the
network structure.

###  3.3 High signalling entropy of intra-cellular signaling hubs is a
hallmark of drug resistance

Figure 6: High signalling entropy of intra-cellular hubs confers drug
resistance: Upper panel plots the topological degree of drug targets (x-axis)
against the Spearman rank Correlation Coefficient (SCC) between its local
signalling entropy and drug sensitivity, as assessed over the CCLE samples.
EC=target annotated as extra-cellular, MR=target annotated as membrane
receptor, IC=target annotated as intra-cellular. Left lower panel shows the
difference in SCC values between the IC and EC+MR targets. Wilcoxon rank sum
test P-value given. Right lower panel shows the difference in SCC values
between the IC and EC+MR targets, where now the SCC value were computed
between gene expression and drug sensitivity. Wilcoxon rank sum test P-value
given.

Among the drug-target pairs for which signalling entropy was a significant
predictor of drug sensitivity, we observed a striking non-linear association
with the topological degree of the targets in the network (Fig.6A). In fact,
for hubs in the network, most of which encode nodes located in the
intracellular signaling hierarchy, high signalling entropy was exclusively
associated with drug resistance (negative SCC). Examples included well-known
intracellular signalling hubs like HDAC1, HDAC2, AKT1, TP53, STAT3, MYC). Some
intracellular non-hubs (e.g. CASP9, BCL2L1, BIRC3) also exhibited negative
correlations between signalling entropy and drug sensitivity. Among targets
for which high signalling entropy predicted drug sensitivity, we observed
several membrane receptors (e.g ERBB2, ERBB3, EGFR, MET) and growth factors
(e.g HBEGF, EGF, TGFA). Given that the correlation coefficients were estimated
across samples (cell-lines) and that the underlying network topology is
unchanged between samples, the observed non-linear relation between the
directionality of the correlation and node degree is a highly non-trivial
finding. We also observed a clear dependence on the main signaling domain of
the target, with intracellular hubs showing preferential anti-correlations, in
contrast to growth factors and membrane receptors which exhibited positive and
negative correlations in equal proportion (Fig.6B). Thus, we can conclude from
these analyses that cancer associated changes to the interaction patterns of
intra-cellular hubs are key determinants of drug resistance. In particular, in
cancers where the local entropy at these hubs is increased, as a result of
increased promiscuous signaling, drugs targeting these hubs are less likely to
be effective.

##  4 Discussion and Conclusions

Here we have advocated a fairly novel methodological framework, based on the
notion of signalling entropy, to help analyze and interpret functional omic
data sets. The method uses a structural network, i.e. a PPI network, from the
outset, and integrates this with gene expression data, using local and global
signalling entropy measures to estimate the amount of uncertainty in the
network signaling patterns. We made the case as to why uncertainty or entropy
might be a key concept to consider when analysing biological data. In previous
work [9, 38], we showed how signalling entropy can be used to estimate the
differentiation potential of a cellular sample in the context of normal
differentiation processes, as well as demonstrating that signalling entropy
also provides a highly accurate discriminator of cancer phenotypes.
In this study we focused on a novel application of signalling entropy to
understanding cellular robustness in the context of cancer drug sensitivity
screens. Our main findings are that (i) local signalling entropy measures add
predictive value over models that only use gene expression, (ii) that the
local entropy of drug targets generally correlates positively with drug
resistance, and (iii) that increased local entropy of intra-cellular hubs in
cancer cells is a key hallmark of drug resistance.
These results are consistent and suggestive of an underlying entropy-
robustness correlation theorem, as envisaged by previous authors [44]. Here,
we provided additional empirical justification for such a theorem, using both
simulated as well as real data, and using drug sensitivity measures as proxies
for local robustness measures. A more detailed theoretical analysis of local
mixing and entropy rates and incorporation of additional information (e.g.
phosphorylation states of kinases, protein expression,..etc) when estimating
entropies on real data, will undoubtedly lead to further improvements in our
systems-level understanding of how entropy/uncertainty dictates cellular
phenotypes. From a practical perspective, we have already shown in previous
work [38] how local network entropies could be used to identify key signaling
pathways in differentiation. It will therefore be interesting in future to
apply the signalling entropy framework in more detail to identify specific
signaling nodes/pathways underlying drug sensitivity/resistance.
Our results have also confirmed the importance of network topology (e.g.
hubness and therefore scale-freeness) in dictating drug resistance patterns.
Thus, it will be interesting to continue to explore the deep relation between
topological features such as scale-freeness and hierarchical modularity in
relation to the gene expression patterns seen in normal and disease
physiology. It is entirely plausible that, although our current data and
network models are only mere caricatures of the real biological networks, that
underlying fundamental systems biology principles characterising cellular
phenotypes can still be gleaned from this data. Indeed, our discovery that
increased signalling entropy correlates with drug resistance demonstrates that
such fundamental principles can already be inferred from existing data
resources.
It will also be of interest to consider other potential applications of the
signaling entropy method. For instance, one application could be to the
identification of functional driver aberrations in cancer. This would first
use epigenomic (e.g. DNA methylation and histone modification profiles) and
genomic information (SNPs, CNVs) together with matched gene or protein
expression data to identify functional epigenetic/genetic aberrations in
cancer. Signalling entropy would subsequently be used as a means of
identifying those aberrations which also cause a fundamental rewiring of the
network. With multi-dimensional matched omic data readily available from the
TCGA/ICGC, this represents another potentially important application of the
signalling entropy method. Another important future application of the
signaling entropy method would be to single-cell data, which is poised to
become ever more important [59]. So far, all signaling entropy analyses have
been performed on cell populations, yet single-cell analysis will be necessary
to disentangle the entropies at the single-cell and population-cell levels.
In summary, we envisage that signalling entropy will become a key concept in
future analyses and interpretation of biological data.

## Acknowledgements

AET is supported by the Chinese Academy of Sciences, Shanghai Institute for
Biological Sciences and the Max-Planck Gesellshaft. RK and PS acknowledge
funding under FP7/2007-2013/grant agreement nr. 290038.

##  Appendix A The protein protein interaction (PPI) network

We downloaded the complete human protein interaction network from Pathway
Commons (www.pathwaycommons.org) (Jun.2012) [24], which brings together
protein interactions from several distinct sources. We built a protein protein
interaction (PPI) network from integrating the following sources: the Human
Protein Reference Database (HPRD) [23], the National Cancer Institute Nature
Pathway Interaction Database (NCI-PID) (pid.nci.nih.gov), the Interactome
(Intact) http://www.ebi.ac.uk/intact/ and the Molecular Interaction Database
(MINT) http://mint.bio.uniroma2.it/mint/. Protein interactions in this network
include physical stable interactions such as those defining protein complexes,
as well as transient interactions such as post-translational modifications and
enzymatic reactions found in signal transduction pathways, including 20 highly
curated immune and cancer signaling pathways from NetPath (www.netpath.org)
[60]. We focused on non-redundant interactions, only included nodes with an
Entrez gene ID annotation and focused on the maximally conntected component,
resulting in a connected network of 10,720 nodes (unique Entrez IDs) and
152,889 documented interactions.

##  Appendix B Mixing rate results for a general random walk on a connected
graph

Suppose that we have a connected graph with an ergodic Markov chain defined on
it, given by a stochastic matrix $P$ with stationary distribution $\pi$ (i.e.
$\pi P=\pi$). We further assume that the detailed balance equation,
$\pi_{i}p_{ij}=\pi_{j}p_{ji}$ holds. Defining the diagonal matrix
$\Pi\equiv\mathrm{diag}(\pi_{1},...,\pi_{N})$, the detailed balance equation
can be rewritten as $\Pi P=P^{T}\Pi$. That the matrix $P$ is stochastic means
that each row of $P$ sums to 1. Equivalently, the vector with all unit
entries, 1, is a right eigenvector of $P$, i.e. $P1=1$. Note also that the
stationary distribution $\pi$ is a left eigenvector of $P$ with eigenvalue 1.
The Perron-Frobenius theorem further implies that all other eigenvalues are
less than 1 in magnitude. If detailed balance holds, all eigenvalues are also
real.
The global mixing rate can be defined by considering the rate at which the
node visitation probabilities of a random walker approaches that of the
stationary distribution, independently of the starting position. Formally, if
we let $Q_{i}(t)$ denote the probability that at time $t$ we find the walker
at node $i$, then the mixing rate, $\mu_{R}$, is defined by [61]

$\mu_{R}=\lim_{t\rightarrow\infty}\sup\max_{i}{|Q_{i}(t)-\pi_{i}|^{1/t}}.$

Denoting by $Q(t)$ the column vector with elements $Q_{i}(t)$, one can write

$Q(t)=(P^{t})^{T}Q(0).$

To determine the elements, $Q_{i}(t)$, of this vector, it is convenient to
first introduce the matrix $M=\Pi^{\frac{1}{2}}P\Pi^{-\frac{1}{2}}$. This is
because $M^{t}=\Pi^{\frac{1}{2}}P^{t}\Pi^{-\frac{1}{2}}$, and so $P^{t}$ can
be rewritten in terms of $M$, but also because $M$ satisfies the following
lemma:

######  Lemma 1.

$M$ has the same eigenvalues as $P$ and if $u_{a}$ is an eigenvector of $M$,
then $r_{a}=\Pi^{-\frac{1}{2}}u_{a}$ and $l_{a}=\Pi^{\frac{1}{2}}u_{a}$ are
right and left eigenvectors of $P$.

###### Proof of Lemma 1.

Suppose that $Mu_{a}=\lambda_{a}u_{a}$. It then follows that

$P(\Pi^{-\frac{1}{2}}u_{a})=\lambda_{a}(\Pi^{-\frac{1}{2}}u_{a}).$

In the case of the left-eigenvector, multiply $Mu_{a}=\lambda_{a}u_{a}$ from
the left with $\Pi^{\frac{1}{2}}$. Then,

$\Pi P\Pi^{-1}(\Pi^{\frac{1}{2}}u_{a})=\lambda_{a}(\Pi^{\frac{1}{2}}u_{a}).$

Detailed balance implies that $\Pi P=P^{T}\Pi$, so
$P^{T}l_{a}=\lambda_{a}l_{a}$. Taking the transpose of this implies that
$l_{a}$ is indeed a left-eigenvector of $P$. ∎

The significance of the above lemma becomes clear in light of the detailed
balance equation, which implies that $M=M^{T}$, and so $M$ and $M^{t}$ can be
orthogonally diagonalized. In particular, we can express $Q_{i}(t)$ in terms
of the eigenvalue decomposition of $M^{t}$, as

$Q_{i}(t)=\sum_{a}{q_{a}|\lambda_{a}|^{t}\tilde{u}_{ai}\pi_{i}^{1/2}}$

, where $q_{a}=\sum_{j}{\tilde{u}_{aj}\pi_{j}^{-1/2}Q_{j}(0)}$ and where
$\tilde{u}_{a}$ is the $a$’th unit norm eigenvector of $M$. Since
$\Pi^{\frac{1}{2}}1$ is the top unit norm eigenvector of $M$ (using previous
lemma), which has an eigenvalue of 1, it follows that $q_{1}=1$ and hence that

$Q_{i}(t)=\pi_{i}+\sum_{a\geq
2}{q_{a}|\lambda_{a}|^{t}\tilde{u}_{ai}\pi_{i}^{1/2}}.$

It follows that

$|Q_{i}(t)-\pi_{i}|=\sum_{a\geq
2}{q_{a}|\lambda_{a}|^{t}\tilde{u}_{ai}\pi_{i}^{1/2}}.$

Since $1=|\lambda_{1}|\geq|\lambda_{2}|\geq|\lambda_{3}|\dots$, we can
conclude that as $t\rightarrow\infty$, the rate at which $Q_{i}$ approaches
the stationary value $\pi_{i}$ is determined by the modulus of the second
largest eigenvalue (the Second Largest Eigenvalue Modulus-SLEM). The global
mixing rate $\mu_{R}$ can thus be estimated as

$\mu_{R}\approx-\log{|\lambda_{2}|}=-\log{\mathrm{SLEM}}$

.

##  Appendix C Entropy and mixing rates in simulated weighted Erdös-Renyi
graphs

For large regular graphs of degree $d$, the mixing rate is proportional to
$\log{d}$ [56] and thus directly proportional to the entropy rate
($SR=\log{d}$ for a regular graph of degree $d$). By extrapolation, we can
thus reasonably assume that for any sufficiently large graph with a tight
degree distribution, such as random Erdos-Renyi (ER) graphs, that the entropy
and mixing rates will also be correlated, albeit perhaps only in an average
ensemble sense. The analytical demonstration of this is beyond the scope of
this work. Hence, we adopt a computational empirical approach to see if the
entropy and mixing rates may indeed be correlated. In detail (and without loss
of generality concerning the end result) we considered ER graphs of size 100
nodes and average degree 10 (other values for these parameters gave similar
answers). We built an ensemble of 100 distinct ER graphs, and for each of
these we constructed a family of weighted networks, parameterised by a
parameter $\epsilon$ which controls the level of signalling entropy.
Specifically, we gradually shifted the weight distribution around each node
$i$ (with degree $d_{i}$), in such a way that $p_{ij}=\epsilon/d_{i}$ for
$j\neq k$ and $p_{ij}=1-\frac{\epsilon}{d_{i}}(d_{i}-1)$ for $j=k$, with
$0\leq\epsilon\leq 1$ and with $k$ labeling a randomly chosen neighbor of node
$i$. Thus, $\epsilon$ is a parameter that directly controls the uncertainty in
the information flow, i.e. the entropy rate, since for $\epsilon=1$ we have
that $p_{ij}=A_{ij}/d_{i}$ ($A_{ij}$ is the symmetric adjacency matrix of the
graph), whilst for $\epsilon=0$, $p_{ij}=\delta_{ik}$, i.e. the information
flow from node $i$ can only proceed along one node (node $k$). Thus, one would
expect the entropy rate to decrease as $\epsilon$ is decreased to zero. For
each value of $\epsilon$ we thus have an ensemble of 100 ER-graphs, for each
of which the entropy and mixing rates can be computed. Finally, at each value
of $\epsilon$ the entropy and mixing rates are averaged over the ensemble.

##  Appendix D Entropy and mixing rates in simulated weighted subgraphs of a
PPI network

The analysis described above was performed also for maximally connected
subnetworks generated from the underlying realistic PPI network described
earlier. Specifically, we randomly subsampled 1000 nodes from the 8038 node
PPI network, and extracted the maximally connected subnetwork, which resulted
(on average) in a subnetwork of approximately 500 nodes. A family of
stochastic matrices of variable entropy rates were constructed as explained
above and for each resulting weighted network we estimated the entropy and
mixing rates. Finally, ensemble averages over 100 different realisations were
computed.

##  Appendix E The Cancer Cell Line Encyclopedia (CCLE) data

We used the gene expression data and drug sensitivity profiles as provided in
the previous publication [57]. Briefly, integration of the gene expression
data with our PPI network resulted in a maximally connected component
consisting of 8038 genes/proteins. There were 488 cell-lines with gene
expression and drug sensitivity profiles for 24 drugs. As a measure of drug
response we used the Activity Area [57] since this measure gave optimal
results when correlating drug response to gene expression levels of well
established drug targets.

## References

  * [1] L. Hood, L. Rowen, The human genome project: big science transforms biology and medicine, Genome Med 5 (9) (2013) 79.
  * [2] N. TCGA, Integrated genomic analyses of ovarian carcinoma, Nature 474 (7353) (2011) 609–615.
  * [3] D. Hanahan, R. A. Weinberg, Hallmarks of cancer: the next generation, Cell 144 (5) (2011) 646–674.
  * [4] A. L. Barabasi, N. Gulbahce, J. Loscalzo, Network medicine: a network-based approach to human disease, Nat Rev Genet 12 (1) (2011) 56–68.
  * [5] P. Creixell, E. M. Schoof, J. T. Erler, R. Linding, Navigating cancer network attractors for tumor-specific therapy, Nat. Biotechnol. 30 (9) (2012) 842–848.
  * [6] G. Weng, U. S. Bhalla, R. Iyengar, Complexity in biological signaling systems, Science 284 (5411) (1999) 92–96.
  * [7] A. L. Barabasi, Network science, Philos Trans A Math Phys Eng Sci 371 (1987) (2013) 20120375.
  * [8] A. E. Teschendorff, S. Severini, Increased entropy of signal transduction in the cancer metastasis phenotype, BMC Syst Biol 4 (2010) 104.
  * [9] J. West, G. Bianconi, S. Severini, A. E. Teschendorff, Differential network entropy reveals cancer system hallmarks, Sci Rep 2 (2012) 802.
  * [10] A. L. Barabasi, Z. N. Oltvai, Network biology: understanding the cell’s functional organization, Nat Rev Genet 5 (2) (2004) 101–113.
  * [11] X. Robin, P. Creixell, O. Radetskaya, C. Costa Santini, J. Longden, R. Linding, Personalized network-based treatments in oncology, Clin. Pharmacol. Ther.
  * [12] J. T. Erler, R. Linding, Network medicine strikes a blow against breast cancer, Cell 149 (4) (2012) 731–733.
  * [13] M. J. Lee, A. S. Ye, A. K. Gardino, A. M. Heijink, P. K. Sorger, G. MacBeath, M. B. Yaffe, Sequential application of anticancer drugs enhances cell death by rewiring apoptotic signaling networks, Cell 149 (4) (2012) 780–794.
  * [14] M. S. Carro, W. K. Lim, M. J. Alvarez, R. J. Bollo, X. Zhao, E. Y. Snyder, E. P. Sulman, S. L. Anne, F. Doetsch, H. Colman, A. Lasorella, K. Aldape, A. Califano, A. Iavarone, The transcriptional network for mesenchymal transformation of brain tumours, Nature 463 (7279) (2010) 318–325.
  * [15] C. Lefebvre, G. Rieckhof, A. Califano, Reverse engineering human regulatory networks, Wiley Interdiscip Rev Syst Biol Med 4 (4) (2012) 311–25.
  * [16] E. Segal, N. Friedman, D. Koller, A. Regev, A module map showing conditional activity of expression modules in cancer, Nat Genet 36 (10) (2004) 1090–1098.
  * [17] R. Opgen-Rhein, K. Strimmer, From correlation to causation networks: a simple approximate learning algorithm, BMC Syst Biol 1 (2007) 37.
  * [18] B. Barzel, A. L. Barabasi, Network link prediction by global silencing of indirect correlations, Nat Biotechnol 31 (8) (2013) 720–5.
  * [19] S. Feizi, D. Marbach, M. Medard, M. Kellis, Network deconvolution as a general method to distinguish direct dependencies in, Nat Biotechnol 31 (8) (2013) 726–33.
  * [20] H. Y. Chuang, E. Lee, Y. T. Liu, D. Lee, T. Ideker, Network-based classification of breast cancer metastasis, Mol Syst Biol 3 (2007) 140.
  * [21] J. Dutkowski, T. Ideker, Protein networks as logic functions in development and cancer, PLoS Comput Biol 7 (9) (2011) e1002180.
  * [22] K. Mitra, A. R. Carvunis, S. K. Ramesh, T. Ideker, Integrative approaches for finding modular structure in biological networks, Nat Rev Genet 14 (10) (2013) 719–32.
  * [23] T. S. Prasad, K. Kandasamy, A. Pandey, Human protein reference database and human proteinpedia as discovery tools for systems biology, Methods Mol Biol 577 (2009) 67–79.
  * [24] E. G. Cerami, B. E. Gross, E. Demir, I. Rodchenkov, O. Babur, N. Anwar, N. Schultz, G. D. Bader, C. Sander, Pathway commons, a web resource for biological pathway data, Nucleic Acids Res 39 (Database) (2011) D685–D690.
  * [25] M. Vidal, M. E. Cusick, A. L. Barabasi, Interactome networks and human disease, Cell 144 (6) (2011) 986–998.
  * [26] I. Ulitsky, R. Shamir, Identification of functional modules using network topology and high-throughput data, BMC Syst Biol 1 (2007) 8.
  * [27] I. W. Taylor, R. Linding, D. Warde-Farley, Y. Liu, C. Pesquita, D. Faria, S. Bull, T. Pawson, Q. Morris, J. L. Wrana, Dynamic modularity in protein interaction networks predicts breast cancer outcome, Nat Biotechnol 27 (2) (2009) 199–204.
  * [28] J. West, S. Beck, X. Wang, A. E. Teschendorff, An integrative network algorithm identifies age-associated differential methylation interactome hotspots targeting stem-cell differentiation pathways, Sci Rep 3 (2013) 1630.
  * [29] J. Roy, C. Winter, Z. Isik, M. Schroeder, Network information improves cancer outcome prediction, Brief Bioinform.
  * [30] C. Winter, G. Kristiansen, S. Kersting, J. Roy, D. Aust, T. Knosel, P. Rummele, B. Jahnke, V. Hentrich, F. Ruckert, M. Niedergethmann, W. Weichert, M. Bahra, H. J. Schlitt, U. Settmacher, H. Friess, M. Buchler, H. D. Saeger, M. Schroeder, C. Pilarsky, R. Grutzmann, Google goes cancer: improving outcome prediction for cancer patients by, PLoS Comput Biol 8 (5) (2012) e1002511.
  * [31] S. Suthram, A. Beyer, R. M. Karp, Y. Eldar, T. Ideker, eqed: an efficient method for interpreting eqtl associations using protein, Mol Syst Biol 4 (2008) 162.
  * [32] Y. A. Kim, S. Wuchty, T. M. Przytycka, Identifying causal genes and dysregulated pathways in complex diseases, PLoS Comput Biol 7 (3) (2011) e1001095.
  * [33] K. Komurov, M. A. White, P. T. Ram, Use of data-biased random walks on graphs for the retrieval of context-specific networks from genomic data, PLoS Comput Biol 6 (8).
  * [34] K. Komurov, J. T. Tseng, M. Muller, E. G. Seviour, T. J. Moss, L. Yang, D. Nagrath, P. T. Ram, The glucose-deprivation network counteracts lapatinib-induced toxicity in, Mol Syst Biol 8 (2012) 596.
  * [35] S. Bandyopadhyay, M. Mehta, D. Kuo, M. K. Sung, R. Chuang, E. J. Jaehnig, B. Bodenmiller, K. Licon, W. Copeland, M. Shales, D. Fiedler, J. Dutkowski, A. Guénolé, H. van Attikum, K. M. Shokat, R. D. Kolodner, W. K. Huh, R. Aebersold, M. C. Keogh, N. J. Krogan, T. Ideker, Rewiring of genetic networks in response to dna damage, Science 330 (6009) (2010) 1385–1389.
  * [36] A. Califano, Rewiring makes the difference, Mol Syst Biol 7 (2011) 463.
  * [37] T. Ideker, N. J. Krogan, Differential network biology, Mol Syst Biol 8 (2012) 565\.
  * [38] C. R. Banerji, D. Miranda-Saavedra, S. Severini, M. Widschwendter, T. Enver, J. X. Zhou, A. E. Teschendorff, Cellular network entropy as the energy potential in waddington’s differentiation landscape, Sci Rep 3 (2013) 3039.
  * [39] H. H. Chang, M. Hemberg, M. Barahona, D. E. Ingber, S. Huang, Transcriptome-wide noise controls lineage choice in mammalian progenitor cells, Nature 453 (7194) (2008) 544–547.
  * [40] B. D. Macarthur, I. R. Lemischka, Statistical mechanics of pluripotency, Cell 154 (3) (2013) 484–489.
  * [41] C. Furusawa, K. Kaneko, Chaotic expression dynamics implies pluripotency: when theory and experiment meet, Biol Direct. 15 (4) (2009) 17.
  * [42] C. Furusawa, K. Kaneko, A dynamical systems view of stem cell biology, Science 338 (6104) (2012) 215–7.
  * [43] T. I. Lee, R. G. Jenner, L. A. Boyer, M. G. Guenther, S. S. Levine, R. M. Kumar, B. Chevalier, S. E. Johnstone, M. F. Cole, K. Isono, H. Koseki, T. Fuchikami, K. Abe, H. L. Murray, J. P. Zucker, B. Yuan, G. W. Bell, E. Herbolsheimer, N. M. Hannett, K. Sun, D. T. Odom, A. P. Otte, T. L. Volkert, D. P. Bartel, D. A. Melton, D. K. Gifford, R. Jaenisch, R. A. Young, Control of developmental regulators by polycomb in human embryonic stem cells, Cell 125 (2) (2006) 301–313.
  * [44] L. Demetrius, V. M. Grundlach, G. Ochs, Complexity and demographic stability in population models, Theo Pop Biol 65 (2004) 211–225.
  * [45] L. Demetrius, T. Manke, Robustness and network evolution-an entropic principle, Physica A 346 (2005) 682–696.
  * [46] T. Manke, L. Demetrius, M. Vingron, An entropic characterization of protein interaction networks and cellular robustness, J R Soc Interface 3 (11) (2006) 843–850.
  * [47] J. Stelling, U. Sauer, Z. Szallasi, F. J. Doyle, J. Doyle, Robustness of cellular functions, Cell 118 (6) (2004) 675–685.
  * [48] J. D. Han, N. Bertin, T. Hao, D. S. Goldberg, G. F. Berriz, L. V. Zhang, D. Dupuy, A. J. Walhout, M. E. Cusick, F. P. Roth, M. Vidal, Evidence for dynamically organized modularity in the yeast protein-protein interaction network, Nature 430 (6995) (2004) 88–93.
  * [49] J. Li, Z. Yuan, Z. Zhang, The cellular robustness by gene redundancy in budding yeast, PLoS Genet 6 (11) (2010) e1001187.
  * [50] J. A. Engelman, L. Chen, X. Tan, K. Crosby, A. R. Guimaraes, R. Upadhyay, M. Maira, K. McNamara, S. A. Perera, Y. Song, L. R. Chirieac, R. Kaur, A. Lightbown, J. Simendinger, T. Li, R. F. Padera, C. Garcia-Echeverria, R. Weissleder, U. Mahmood, L. C. Cantley, K. K. Wong, Effective use of PI3K and MEK inhibitors to treat mutant Kras G12D and PIK3CA H1047R murine lung cancers, Nat. Med. 14 (12) (2008) 1351–1356.
  * [51] T. Sjoblom, S. Jones, L. D. Wood, D. W. Parsons, J. Lin, T. D. Barber, D. Mandelker, R. J. Leary, J. Ptak, N. Silliman, S. Szabo, P. Buckhaults, C. Farrell, P. Meeh, S. D. Markowitz, J. Willis, D. Dawson, J. K. Willson, A. F. Gazdar, J. Hartigan, L. Wu, C. Liu, G. Parmigiani, B. H. Park, K. E. Bachman, N. Papadopoulos, B. Vogelstein, K. W. Kinzler, V. E. Velculescu, The consensus coding sequences of human breast and colorectal cancers, Science 314 (5797) (2006) 268–274.
  * [52] L. D. Wood, D. W. Parsons, S. Jones, J. Lin, T. Sjöblom, R. J. Leary, D. Shen, S. M. Boca, T. Barber, J. Ptak, N. Silliman, S. Szabo, Z. Dezso, V. Ustyanksky, T. Nikolskaya, Y. Nikolsky, R. Karchin, P. A. Wilson, J. S. Kaminker, Z. Zhang, R. Croshaw, J. Willis, D. Dawson, M. Shipitsin, J. K. Willson, S. Sukumar, K. Polyak, B. H. Park, C. L. Pethiyagoda, P. V. Pant, D. G. Ballinger, A. B. Sparks, J. Hartigan, D. R. Smith, E. Suh, N. Papadopoulos, P. Buckhaults, S. D. Markowitz, G. Parmigiani, K. W. Kinzler, V. E. Velculescu, B. Vogelstein, The genomic landscapes of human breast and colorectal cancers, Science 318 (5853) (2007) 1108–1113.
  * [53] B. Vogelstein, N. Papadopoulos, V. E. Velculescu, S. Zhou, L. A. Diaz, K. W. Kinzler, Cancer genome landscapes, Science 339 (6127) (2013) 1546–58.
  * [54] V. Latora, M. Baranger, Kolmogorov-sinai entropy rate versus physical entropy, Phys Rev Lett 82 (3) (1999) 520–524.
  * [55] J. Gomez-Gardenes, V. Latora, Entropy rate of diffusion processes on complex networks, Phys Rev E Stat Nonlin Soft Matter Phys 78 (6) (2008) 065102.
  * [56] J. Friedman, On the second eigenvalue and random walks in random d-regular graphs, Combinatorica 11 (1991) 331–262.
  * [57] J. Barretina, G. Caponigro, N. Stransky, K. Venkatesan, A. A. Margolin, S. Kim, C. J. Wilson, J. Lehár, G. V. Kryukov, D. Sonkin, A. Reddy, M. Liu, L. Murray, M. F. Berger, J. E. Monahan, P. Morais, J. Meltzer, A. Korejwa, J. Jané-Valbuena, F. A. Mapa, J. Thibault, E. Bric-Furlong, P. Raman, A. Shipway, I. H. Engels, J. Cheng, G. K. Yu, J. Yu, P. Aspesi Jr, M. de Silva, K. Jagtap, M. D. Jones, L. Wang, C. Hatton, E. Palescandolo, S. Gupta, S. Mahan, C. Sougnez, R. C. Onofrio, T. Liefeld, L. MacConaill, W. Winckler, M. Reich, N. Li, J. P. Mesirov, S. B. Gabriel, G. Getz, K. Ardlie, V. Chan, V. E. Myer, B. L. Weber, J. Porter, M. Warmuth, P. Finan, J. L. Harris, M. Meyerson, T. R. Golub, M. P. Morrissey, W. R. Sellers, R. Schlegel, L. A. Garraway, The cancer cell line encyclopedia enables predictive modelling of anticancer drug sensitivity, Nature 483 (7391) (2012) 603–607.
  * [58] J. Ahmed, T. Meinel, M. Dunkel, M. S. Murgueitio, R. Adams, C. Blasse, A. Eckert, S. Preissner, R. Preissner, Cancerresource: a comprehensive database of cancer-relevant proteins and compound interactions supported by experimental knowledge, Nucleic Acids Res 39 (Database) (2011) D960–D967.
  * [59] M. R. Speicher, Single-cell analysis: toward the clinic, Genome Med 5 (8) (2013) 74.
  * [60] K. Kandasamy, S. S. Mohan, R. Raju, S. Keerthikumar, G. S. Kumar, A. K. Venugopal, D. Telikicherla, J. D. Navarro, S. Mathivanan, C. Pecquet, S. K. Gollapudi, S. G. Tattikota, S. Mohan, H. Padhukasahasram, Y. Subbannayya, R. Goel, H. K. Jacob, J. Zhong, R. Sekhar, V. Nanjappa, L. Balakrishnan, R. Subbaiah, Y. L. Ramachandra, B. A. Rahiman, T. S. Prasad, J. X. Lin, J. C. Houtman, S. Desiderio, J. C. Renauld, S. N. Constantinescu, O. Ohara, T. Hirano, M. Kubo, S. Singh, P. Khatri, S. Draghici, G. D. Bader, C. Sander, W. J. Leonard, A. Pandey, Netpath: a public resource of curated signal transduction pathways, Genome Biol 11 (1) (2010) R3.
  * [61] L. Lovasz, Random walks on graphs: A survey, in: Combinatorics, Paul Erdos is Eighty, Vol. 2 of Mathematical Studies, Bolyai Society, Keszthely, Hungary, 1993, pp. 1–46.

## Figures


# Optimizing Warfarin Dosing Using Contextual Bandit: An Offline Policy
Learning and Evaluation Method

Yong Huang1,∗, Charles A. Downs3, and Amir M. Rahmani1,2,4 1Dept. of Computer
Science, University of California, Irvine. 2School of Nursing, University of
California, Irvine. 3School of Nursing and Health Sciences, University of
Miami. 4Institute for Future Health (IFH), University of California, Irvine.
(*correspondence email<EMAIL_ADDRESS>
###### Abstract

Warfarin, an anticoagulant medication, is formulated to prevent and address
conditions associated with abnormal blood clotting, making it one of the most
prescribed drugs globally. However, determining the suitable dosage remains
challenging due to individual response variations, and prescribing an
incorrect dosage may lead to severe consequences. Contextual bandit and
reinforcement learning have shown promise in addressing this issue. Given the
wide availability of observational data and safety concerns of decision-making
in healthcare, we focused on using exclusively observational data from
historical policies as demonstrations to derive new policies; we utilized
offline policy learning and evaluation in a contextual bandit setting to
establish the optimal personalized dosage strategy. Our learned policies
surpassed these baseline approaches without genotype inputs, even when given a
suboptimal demonstration, showcasing promising application potential.

##  I INTRODUCTION

Warfarin, introduced in the 1950s, is designed to prevent and manage
conditions related to abnormal blood clotting, such as deep vein thrombosis,
pulmonary embolism, and stroke [1]. It is estimated to be prescribed over 11
million times in the United States in 2021 [2]. Despite its extensive use,
determining the appropriate dosage poses a challenge due to variations in
individual responses. Factors like age, weight, genetics, and liver and kidney
function, among others, play pivotal roles in the drug’s metabolism, making
identification of the optimal dosage challenging [3].

Choosing an incorrect warfarin dosage can have dangerous consequences, leading
to life-threatening excessive bleeding. Consequently, there is a growing
interest in developing policies for determining the optimal dose. In practice,
clinicians typically initiate treatment with an initial dosage, then closely
monitor the patient’s responses and make dosing adjustments. This iterative
process continues until the optimal therapeutic dosage is identified,
typically taking weeks [4].

Ideally, identifying the correct dose without requiring an extensive
exploration procedure would save critical time. Numerous efforts have been
proposed to predict the optimal initial dosage. For instance, the Warfarin
Clinical Dosing Algorithm (WCDA) [5], a linear model, considers variables like
weight, age, and medications. Similarly, the Warfarin Pharmacogenetic Dosing
Algorithm (WPDA) is another linear model that introduces additional genotype
inputs [5].

These simplistic linear models, however, often fall short of capturing the
intricate relationships between variables, resulting in only marginally
satisfactory success rates for predicting the correct therapeutic dose.
Furthermore, the WPDA algorithms necessitate genotype inputs for more accurate
predictions, posing challenges in their application to a broader population
where genotype information may be unavailable. This underscores the need to
develop more advanced approaches to assist clinicians in decision-making when
determining the appropriate dosage.

Recent advancements in contextual bandit and reinforcement learning have
demonstrated promising results in enhancing decision-making across various
applications [6, 7], and several pilot studies have attempted to use online
approaches for optimizing warfarin dosing, achieving promising results [8, 9].
However, a critical limiting factor of online approaches is their dependence
on generating new data through active exploration. In other words, online
learning approaches typically require generating data from the policy being
improved in the learning process, rendering them less applicable in domains
like healthcare, where safety and ethics are critical, and building
simulations can be challenging.

Alternatively, there is a growing interest in healthcare for offline
reinforcement learning and bandits [10]. Unlike online approaches, offline
methods leverage historical observational data only, typically collected by
another policy, and attempt to derive a new policy from it. Offline approaches
lend themselves to healthcare applications where historical observational data
is widely available. The offline learning nature frees us from building a
simulation environment for exploration, which is challenging and complicated
to extend to more complex decision-making problems, such as sepsis management.
Furthermore, recent theoretical work has demonstrated that learning a policy
merely from observational data that outperforms the demonstration is possible
[11].

In this study, we frame the Warfarin dosing problem within offline contextual
bandit, aiming to maximize expected rewards by identifying a high-quality
dosage policy. We utilize historical observational data from baseline
policies, employing offline learning algorithms to derive improved policies.
Given the limited research on Off-policy evaluation (OPE) methods, we assess
results using three representative OPE estimators. Our learned policy
outperforms demonstration policies, even when the baseline policy is
suboptimal (e.g., random dosage policy). This highlights the effectiveness of
our approach to Warfarin dosing.

Our contributions are three folds:

  * •

To the best of our knowledge, this work is the first to propose a machine-
learning dosing model in an offline fashion.

  * •

Our novel policy is capable of either matching or surpassing the effectiveness
of existing clinical baselines, and it achieves this without requiring
genotype data, thus enhancing its scalability for practical implementation.

  * •

We have explored the validity and reliability of the OPE methods within
clinical settings, offering practical empirical evidence to guide the
application and selection of these estimators.

##  II Methods

###  II-A Dataset

We utilize a publicly accessible patient dataset curated by the
Pharmacogenomics Knowledge Base (PharmGKB) [12]. This dataset comprises
information from 5700 patients who underwent warfarin treatment across the
globe. Each patient’s profile in this dataset includes various features such
as demographics (e.g., gender, race), background details (e.g., height,
weight, medical history), and phenotypic and genotypic information.
Additionally, the dataset includes the actual patient-specific optimal
warfarin doses determined through the physician-guided dose adjustment process
over several weeks. We categorized the therapeutic doses into low (less than
21 mg/week), medium (21-49 mg/week), and high (more than 49 mg/week) following
conventions [5].

###  II-B Problem Formulation

We formulate the problem of optimal warfarin dosing under the contextual
bandit framework. The Dataset $\mathcal{D}$ consists of $\mathcal{N}$
patients, and for each patient, we observe its feature vector
$\mathcal{X}\in\mathcal{R}^{d}$. This represents the available knowledge about
the patient, also known as the context in bandit literature, and it is used to
help determine actions $\mathcal{A}$ chosen by the policy $f$, which has
access to $\mathcal{K}$ actions/arms where the action represents the warfarin
dosage to provide to the patient. As defined in the previous section, we have
$\mathcal{K}=3$ arms corresponding to low, medium, and high warfarin doses. If
the algorithm identifies the correct dosage for the patient, the reward
$\mathcal{R}$ is 1. Otherwise, a reward of 0 is received. The objective is to
optimize the policy such that the action selected by the policy has a maximum
expected reward $E[\mathcal{R}]$. To demonstrate better how offline contextual
bandit works, we present the high-level idea of the workflow in Figure 1.

Figure 1: Workflow of offline learning and evaluation, an essential distinct
between contextual bandit and supervised learning is that the ground truth
optimal action is not revealed to learning and evaluation algorithms, making
it close to real-world decision-making problems where the outcome associated
with the optimal action may be counterfactual and unavailable in observational
data.

###  II-C Offline Policy Learning

We implemented two offline policy learning (OPL) algorithms, the Offset Tree
and the doubly robust estimator [13, 14], to learn the optimal policy. The
offset tree algorithm simplifies the original problem, a censored label
classification problem, to binary importance-weighted classification, which
can be implemented using any binary classifier, such as logistic regression.
Conversely, the doubly robust (DR) estimators require two estimates: one for
the reward, providing an estimate of the reward each arm will give, and
another for the probability, offering an estimate of the probability that the
policy collecting the data assigned to each arm chosen for each observation,
which is used to weight the observed rewards. As long as one of the estimates
is well-specified, the DR estimator is unbiased and, hence, doubly robust.

###  II-D Offline Policy Evaluation

Although the true counterfactual outcome of every possible action is known
under our specific reward formulation, making it possible to evaluate the
performance of the proposed policy directly, it is worthwhile to investigate
it through the lens of OPE. In a more general scenario, the counterfactual
outcome is typically unobservable, and OPE becomes the only tool to quantify
policy quality. Given that limited work has been done to assess these OPE
methods empirically [15], the Warfarin problem provides an excellent
opportunity to help us understand how well OPE works by comparing the
estimates with the oracle ground truth. Therefore, we perform OPE to estimate
the expected reward of the proposed policy and compare these estimates with
the ground truth expected reward.

At a high level, Off-policy evaluation (OPE) offers a statistical framework
for estimating the value function of a specific target policy, typically
utilizing information from another policy that generated the data. We
implemented three policy evaluation methods:

  * •

Rejection Sampling [16]: Reject sampling is an unbiased estimator that
estimates the performance of a policy by collecting some data on which actions
are chosen at random, a reject sampling scheme is performed later on the
policy to be evaluated, namely, if the action agrees with the observation from
random action, it is kept and it is rejected otherwise.

  * •

Doubly Robust [14]: Doubly Robust in evaluation works the same way as
described in the learning section.

  * •

Normalized Capped Importance Sampling [17] Normalized Capped Importance
Sampling estimates rewards of arm choices of a policy from data collected from
the demonstration policy, making corrections by weighting rewards according to
the difference between the estimations of the evaluation policy and
demonstration policy over the actions that were chosen.

##  III EXPERIMENTS AND Discussion

We pre-processed the dataset, which involved normalizing continuous variables
and imputing missing values as necessary. For instance, missing weights were
imputed by using the average value across the dataset. The dataset was then
divided into an 80 percent training set and a 20 percent test set.

We implemented three baseline policies. The first policy is a random dosage
policy where every patient is assigned a random dosage of warfarin. The other
two policies are WCDA and WPDA algorithms. The rewards associated with the
actions selected by these policies were computed. For each baseline policy,
the data tuples consisting of $\langle X,A,R\rangle$ were used in the training
of offline learning algorithms to derive new policies. Subsequently, we
calculated the expected rewards of these learned new policies on the test set
and compared them with the expected rewards of baseline policies that were
used as demonstration, the OPL algorithm and OPE methods are implemented using
contextual bandit library in python [18].

To account for the randomness introduced by the train-test split, we performed
the splitting process 30 times with different random seeds and reported the
results of these 30 experiments. The offline learning results are presented in
Figure 2; we compared the oracle reward on the left-out test set of two OPL
methods against the baseline demonstration which they are learned from. Our
newly learned policies outperformed their respective demonstrations by a large
margin, notably when presented with the random policy data, a suboptimal
policy, OPL algorithms still managed to learned a policy that matches other
baselines. Furthermore, both offset tree and doubly robust performance are
similar regardless what observational data they learned from, we speculate
that this may be the most ideal results we could achieve given the available
features collected. Importantly, the OPL algorithms showcased superior
performance without requiring genotype inputs, making them more practical for
deployment.

Figure 2: The expected reward of thirty experiments on test sets is presented in a boxplot. In each subfigure, offset tree and doubly robust estimator learn from a corresponding old policy. Table I: OPE estimation of rewards and oracle rewards of new policies learned from random policy demonstration OPE | OT learned policy |  | DR learned policy |
---|---|---|---|---
Reject Sampling | 0.622 (0.545, 0.660) |  | 0.625 (0.578, 0.660) |
DR | 0.619 (0.552, 0.661) |  | 0.625 (0.577, 0.660) |
NCIS | 0.331 (0.308, 0.376) |  | 0.331 (0.308, 0.376) |
oracle | 0.624 (0.598, 0.650) |  | 0.631 (0.612, 0.662) |
Table II: OPE estimation of rewards and oracle rewards of new policy learned from WCDA policy demonstration OPE | OT learned policy |  | DR learned policy |
---|---|---|---|---
Reject Sampling | 0.650 (0.623, 0.676) |  | 0.639 (0.618, 0.656) |
DR | 0.651 (0.600, 0.890) |  | 0.633 (0.610, 0.653) |
NCIS | 0.601 (0.581, 0.625) |  | 0.601 (0.581, 0.625) |
oracle | 0.627 (0.601, 0.656) |  | 0.623 (0.606 ,0.639) |
Table III: OPE estimation of rewards and oracle rewards of new policy learned from WPDA policy demonstration OPE | OT learned policy |  | DR learned policy |
---|---|---|---|---
Reject Sampling | 0.652 (0.633, 0.677) |  | 0.639 (0.620, 0.663) |
DR | 0.737 (0.444, 1.000) |  | 0.631 (0.611, 0.656) |
NCIS | 0.562 (0.545, 0.581) |  | 0.562 (0.545, 0.581) |
oracle | 0.626 (0.602, 0.659) |  | 0.629 (0.609, 0.648) |

Another task of interest in this work is to examine the empirical performance
of OPE methods on newly proposed policies. We implemented a logistic
regression to acquire the distribution over actions generated by the new
policy, which is used to compute reward estimates in doubly robust OPE and
NCIS OPE, and Table I, II,, III detail our experimental results on OPEs. We
presented mean rewards estimation from the 30 experiments with lower and upper
bounds and the actual oracle rewards. It is worth noting that the NCIS appears
to deviate the most from the ground truth, while reject sampling is
empirically the best choice. A potential reason for more significant errors in
NCIS estimation may stem from the inaccurate importance weights we calculated
using logistic regression; reject sampling typically works well when the
sample size is large enough. Additionally, doubly robust methods work best
with continuous actions, and the same dependency on importance weight
computation makes it less practical for real-world applications.

In conclusion, we demonstrated the feasibility of deriving an improved policy
that outperforms its demonstrations using historical observational data. We
empirically evaluated three representative OPE methods in a bandit setting. A
primary limitation of this study is that the feature space of collected
observational data may not be sufficiently representative, potentially
hindering the OPL from achieving optimal performance due to unobserved
confounding factors. Additionally, the empirical OPE evaluation results may
not generalize to other use cases, especially considering the variability in
the properties of these estimators in reinforcement learning settings. We
intend to address these issues in future research.

## References

  * [1] A. E. Rettie and G. Tai, “The pharmocogenomics of warfarin,” _Molecular interventions_ , vol. 6, no. 4, p. 223, 2006.
  * [2] P. Sean P. Kane. [Online]. Available: http://clincalc.com/DrugStats
  * [3] M. Kuruvilla and C. Gurk-Turner, “A review of warfarin dosing and monitoring,” in _Baylor University Medical Center Proceedings_ , vol. 14, no. 3. Taylor & Francis, 2001, pp. 305–306.
  * [4] D. A. Garcia and M. J. Schwartz, “Warfarin therapy: tips and tools for better control,” _Journal of Family Practice_ , vol. 60, no. 2, p. 70, 2011.
  * [5] I. W. P. Consortium, “Estimation of the warfarin dose with clinical and pharmacogenetic data,” _New England Journal of Medicine_ , vol. 360, no. 8, pp. 753–764, 2009.
  * [6] Y. Huang, R. Cao, and A. Rahmani, “Reinforcement learning for sepsis treatment: A continuous action space solution,” in _Machine Learning for Healthcare Conference_. PMLR, 2022, pp. 631–647.
  * [7] A. Raghu, M. Komorowski _et al._ , “Deep reinforcement learning for sepsis treatment,” _arXiv preprint arXiv:1711.09602_ , 2017.
  * [8] A. Vats, “Estimation of warfarin dosage with reinforcement learning,” _arXiv preprint arXiv:2109.07564_ , 2021.
  * [9] H. Xiao, “Online learning to estimate warfarin dose with contextual linear bandits,” _arXiv preprint arXiv:1907.05496_ , 2019.
  * [10] O. Gottesman, F. Johansson _et al._ , “Evaluating reinforcement learning algorithms in observational health settings,” _arXiv preprint arXiv:1805.12298_ , 2018.
  * [11] P. Rashidinejad, B. Zhu _et al._ , “Bridging offline reinforcement learning and imitation learning: A tale of pessimism,” _Advances in Neural Information Processing Systems_ , vol. 34, pp. 11 702–11 716, 2021.
  * [12] J. M. Barbarino, M. Whirl-Carrillo _et al._ , “Pharmgkb: a worldwide resource for pharmacogenomic information,” _Wiley Interdisciplinary Reviews: Systems Biology and Medicine_ , vol. 10, no. 4, p. e1417, 2018.
  * [13] A. Beygelzimer and J. Langford, “The offset tree for learning with partial labels,” in _Proceedings of the 15th ACM SIGKDD international conference on Knowledge discovery and data mining_ , 2009, pp. 129–138.
  * [14] M. Dudík, J. Langford, and L. Li, “Doubly robust policy evaluation and learning,” _arXiv preprint arXiv:1103.4601_ , 2011.
  * [15] C. Voloshin, H. M. Le _et al._ , “Empirical study of off-policy policy evaluation for reinforcement learning,” _arXiv preprint arXiv:1911.06854_ , 2019.
  * [16] L. Li, W. Chu _et al._ , “A contextual-bandit approach to personalized news article recommendation,” in _Proceedings of the 19th international conference on World wide web_ , 2010, pp. 661–670.
  * [17] A. Gilotte, C. Calauzènes _et al._ , “Offline a/b testing for recommender systems,” in _Proceedings of the Eleventh ACM International Conference on Web Search and Data Mining_ , 2018, pp. 198–206.
  * [18] D. Cortes, “Adapting multi-armed bandits policies to contextual bandits scenarios,” _arXiv preprint arXiv:1811.04383_ , 2018.