11institutetext: Department of Computer Science and Information Systems
Teesside University
Middlesbrough, UK

# Integrating splice-isoform expression into genome-scale models characterizes
breast cancer metabolism

Claudio Angione

###### Abstract

Motivation: Despite being often perceived as the main contributors to cell
fate and physiology, genes alone cannot predict cellular phenotype. During the
process of gene expression, $95\%$ of human genes can code for multiple
proteins due to alternative splicing. While most splice variants of a gene
carry the same function, variants within some key genes can have remarkably
different roles. To bridge the gap between genotype and phenotype, condition-
and tissue-specific models of metabolism have been constructed. However,
current metabolic models only include information at the gene level.
Consequently, as recently acknowledged by the scientific community, common
situations where changes in splice-isoform expression levels alter the
metabolic outcome cannot be modeled.
Results: We here propose GEMsplice, the first method for the incorporation of
splice-isoform expression data into genome-scale metabolic models. Using
GEMsplice, we make full use of RNA-Seq quantitative expression profiles to
predict, for the first time, the effects of splice isoform-level changes in
the metabolism of $1455$ patients with $31$ different breast cancer types. We
validate GEMsplice by generating cancer-versus-normal predictions on metabolic
pathways, and by comparing with gene-level approaches and available literature
on pathways affected by breast cancer. GEMsplice is freely available for
academic use at https://github.com/GEMsplice/GEMsplice_code. Compared to
state-of-the-art methods, we anticipate that GEMsplice will enable for the
first time computational analyses at transcript level with splice-isoform
resolution.

This is a pre-copyedited, author-produced version of an article accepted for
publication in Bioinformatics following peer review. The version of record “C.
Angione, ”Integrating splice-isoform expression into genome-scale models
characterizes breast cancer metabolism”, Bioinformatics, btx562,
34(3):494-501, 2018” is available online at:
https://academic.oup.com/bioinformatics/article-abstract/34/3/494/4107935.

##  1 Introduction

The increasing availability of multi-omic datasets has rapidly moved the focus
of many research efforts from data collection to finding effective methods for
data interpretation and analysis. In terms of biomedical results, this step is
therefore currently considered more limiting than data collection itself [1].

Genes and their expression have been the subject of the vast majority of
research studies in computational biology and bioinformatics. However, it is
the interaction of genes, proteins, reactions, and metabolites (different
omics) that shapes the behavior of a cell. When analyzing biological models,
the classical pathway-based perspective has been replaced, in the last 20
years, by a network-based approach, therefore bypassing parametrization and
need for kinetic data, and leading to the generation of genome-scale metabolic
models. These models include thousands of biochemical reactions, often the
full set of reactions known for a given organism, allowing for prediction of
cell phenotype. In this regard, metabolism is the only biological system that
can be fully modeled at genome-scale, and also the closest to the phenotype.

Metabolism is now considered as a driver, rather than a marker, of cancer
onset and proliferation. In breast cancer, metabolic heterogeneity is one of
the causes of poor clinical outcome: although being classified as a single
disease, more than $20$ types of breast cancer exist. A genome-scale analysis
of cancer metabolism captures many effects that could not be identified using
standard data and gene expression analysis. Furthermore, a wide set of tools
for the incorporation of cancer omic data into these models have been
developed, making them suitable for integrating and interpreting the great
amount of data that is being gathered on cancer metabolic alterations [2].

As mentioned above, reduced cost for data gathering and analysis has recently
yielded a rapid increase in the amount of available omics data. RNA-Seq is now
widely used to produce high-throughput data in more detail compared to
microarrays, with a more accurate estimation of transcript levels. As a
result, splice isoform expression levels are now becoming available in cancer
studies, where the same gene can code for multiple proteins due to alternative
RNA splicing before translation. Alterations of specific splice-isoform
expression in some genes can constitute a biomarker for cancer metabolism. In
humans, alternative splicing affects $95\%$ of multiexon genes [3].

However, to date, functional information included in most metabolic models
refers to genes or proteins only. The idea exploited here is that splice
isoform data can be readily integrated and used in conjunction with annotated
genome-scale models in order to constrain and refine existing models. Although
Recon1 [4] was the only human metabolic model that introduced isoform-level
annotations for some genes through bibliographic research, it adopted a custom
annotation for isoforms, and therefore did not allow any mapping to known
identifiers from public databases. For these reasons, such annotations have
been subsequently lost and then simply ignored, mapping transcriptomic data
with gene-level resolution only [5]. Consequently, expression data at the
splice-isoform level has been neglected or simply averaged within the same
gene to approximate the expression at the gene level.

Nevertheless, especially in human metabolic models, the incorporation of
splice isoforms is key to understanding complex diseases like cancer. In this
regard, pyruvate kinase (PK) is a striking example. In fact, the switch to the
second isoform of pyruvate kinase (PKM2) is considered essential for cancer
growth [6]. Conversely, the switch from PKM2 to the main pyruvate kinase
isoform (PKM1) is able to reverse the Warburg effect and could therefore
constitute a therapeutic target. In general, these crucial isoform-level
events, e.g. the switch to minor isoforms, which get overexpressed compared to
the major isoform of a gene, cannot be captured by current metabolic models,
as highlighted in a number of recent reviews [7, 8, 5, 9].

Here we propose GEMsplice (genome-scale metabolic modeling with splice-isoform
integration), the first method for incorporating RNA-Seq data at the splice-
isoform level into a metabolic model. The GEMsplice pipeline exploits, for the
first time, the full potential of the next-generation sequencing technology in
the context of genome-scale metabolic reconstructions. As a result, it enables
more accurate predictions of human metabolic behavior and cancer metabolism.
We show that GEMsplice compares favorably to existing tools for integration of
transcriptomics data at the gene level only. We also validate GEMsplice by
building breast cancer-versus-normal genome-scale models and by comparing
predictions with available results on metabolic pathways affected by breast
cancer. The full GEMsplice pipeline is illustrated in Figure 1. GEMsplice is
made freely available in Matlab/Octave at
https://github.com/GEMsplice/GEMsplice_code, and is fully compatible with the
COBRA 2.0 Toolbox [10].

Figure 1: GEMsplice incorporates RNA-Seq data into genome-scale metabolic
models at the splice-isoform level. Starting from a metabolic model of breast
cancer metabolism [11] (A), gene expression (B) and transcript level
information (C) are incorporated into the model. As a result, for the first
time, we can exploit the full potential of next-generation sequencing in the
context of genome-scale metabolic reconstructions. A set of phenotype-specific
RNA-Seq transcript expression levels in a variety of breast cancer types and
stages from the Cancer RNA-Seq Nexus dataset [12], including data from TCGA,
GEO and SRA (D), are then mapped onto the model using constraint-based
modeling (E). Cancer-specific metabolic models are finally generated and
investigated using multilevel linear programming (F), leading to phenotype
prediction for different types of breast cancer (G).

### Beyond the Warburg effect: metabolism is a key player in cancer formation
and progression

Known molecular markers of cancer can be associated with two main classes: (i)
cell proliferation and (ii) tissue remodeling [13]. Cancer cells show higher
rate of proliferation than normal cells in the tissue from which they
originated. The main goal of metabolism in cancer cells is to keep cell
viability and ensure new biomass production by acquiring nutrients from an
environment where nutrients are often scarce. To this end, cancer cells
exhibit a modified metabolism and an increased need for proteins, energy,
nucleotides and lipids. The two main nutrients used by cancer cells are
glucose and glutamine; both support cell growth and are used to build carbon
intermediates, which are employed in a number of processes that build
macromolecules.

The first difference between normal and cancer metabolism was observed in the
1920s by Otto Warburg [14]. Normal cells take up glucose and perform
glycolysis to obtain pyruvate, which is then transferred to the mitochondrion
and used by oxidative phosphorylation (OXPHOS) in the TCA cycle along with
oxygen to efficiently produce ATP. Conversely, proliferating cancer cells
require higher amount of glucose, but after obtaining pyruvate they
preferentially use it to secrete lactate in the cytoplasm. Strikingly, this
rather inefficient energy production pipeline is used also when the cell is
exposed to oxygen. Preferential use of glycolysis allows faster proliferation
but maximizes the secretion of lactic acid, which damages surrounding cells.
This behavior, called Warburg effect or aerobic glycolysis, facilitates
proliferation and migration of cancer cells.

However, a misinterpretation of this phenomenon by Warburg himself led to the
enduring misconception that cancer cells do not use the TCA cycle to produce
ATP. Indeed, it must be noted that the vast majority of cancer cells still use
the mitochondrion and its TCA cycle to produce a fraction of the required ATP.
In fact, in contrast to quiescent cells, glycolysis does not fuel directly the
TCA cycle, but is essentially decoupled from it [15]. The advantage of
converting excess pyruvate into lactate instead of transferring it into the
mitochondrion is that high glycolytic activity can continue without leading to
excessive flux through the electron transport system and consequent
overproduction of reactive oxygen species (ROS), or excess ATP and NADH
generation, which would in turn repress glycolysis. As a result of this
decoupling, useful glycolytic intermediates (e.g., precursors of serine
biosynthesis) can be generated at higher rates.

Pyruvate kinase (PK), the last step in the glycolytic pathway, is responsible
for keeping the balance between the production of pyruvate for the
mitochondrion (high PK activity) and the production of glycolytic
intermediates for biosynthesis (low PK activity). This is achieved through the
preferential expression of PKM1 or PKM2, two splice isoforms of PK whose
activity respectively supports pyruvate or glycolytic intermediates for
biosynthetic processes. In cancer phenotypes, being mostly controlled by the
isoform PKM2, PK and its ATP product are independent of oxygen, therefore
enabling ATP generation and growth during hypoxia (unlike mitochondrial
respiration, which needs oxygen).

Glucose is not the only key nutrient for cancer cells. Proliferating cells can
consume up to ten times more glutamine than any other amino acid. Glutamine is
used as a carbon source for fatty acid synthesis, in biosynthetic pathways
(especially the biosynthesis of nucleotides), and for TCA cycle intermediates
(although as a secondary source, less than glucose). While in normal cells
glucose/glutamine intake depends mainly on extracellular stimuli, mutations
undergone by cancer cells confer the ability to proliferate with a high degree
of independence, by constantly importing amino acids and glucose from the
extracellular environment [16]. Taken together, these studies show that
metabolism, once believed to be mainly a passive indicator of the state of a
cell, is now widely recognized as a key player in cancer formation and
progression.

##  2 Methods

### Incorporating RNA-Seq data into a breast cancer model at the splice-
isoform level

To model splice-isoform data in the breast cancer metabolic model, $31$ RNA-
Seq expression profiles of breast cancer were considered from Cancer RNA-Seq
Nexus [12], which provides phenotype-specific transcript expression levels in
a variety of cancer types and stages. In Cancer RNA-Seq Nexus, data were
processed as follows. For GEO/SRA samples, Bowtie [17] was used to align RNA-
Seq reads to the reference human transcriptome from GENCODE [18]; eXpress [19]
was then used to obtain FPKM isoform abundances. For TCGA samples, starting
from Level 3 RNA-Seq v2 “tau values” quantified through RSEM [20], TPM values
were obtained through multiplication by $10^{6}$.

The expression values, at this point, are measured in FPKM for GEO/SRA, and in
TPM for TCGA samples. FPKMs were then converted to TPMs using the following
formula for all the expression values $j$ in each sample $i$:
$\mbox{TPM}_{ij}=10^{6}\frac{\mbox{FPKM}_{ij}}{\lambda_{i}\sum_{j}{\mbox{FPKM}_{ij}}}$,
where $\sum_{j}{\mbox{FPKM}_{ij}}$ is the sum of the FPKM values of all the
transcripts in the $i$th sample, and $\lambda_{i}$ a sample-specific scaling
factor accounting for the fact that the number of transcripts varies across
the different samples, $\lambda_{i}=$ (number of transcripts in the largest
sample)$/$(number of transcripts in sample $i$). We remark that TPM (and not
FPKM) expression values were used as a starting point to generate the
metabolic outcome across samples because they are proportional to the
abundance and independent of the mean expressed transcript length, therefore
making them more suitable to comparisons across samples [20].

Transcript annotations in Cancer RNA-Seq Nexus consist of Ensembl and UCSC
IDs. These were converted to RefSeq annotations using BioMart [21]. However,
the breast cancer model is annotated only with gene-level Entrez IDs. To
expand these and generate transcript-level IDs with splice-isoform resolution,
RefSeq IDs for $141$ transcripts were retrieved from the SBML source of Recon1
using a custom script, as they were not included in the final Matlab version
of the model (see File S1). These newly generated transcript-level RefSeq IDs
and the gene-level default Entrez IDs of the remaining genes in the model were
used to associate each transcript in the model with the corresponding
expression value in the $31$ RNA-Seq cancer profiles. The profiles were
finally mapped onto the breast cancer metabolic model (see Figure 1, and the
following Methods subsection for details on how the mapping was achieved).

### Phenotype predictions from RNA-Seq data

Breast cancer is a multifactorial disease and, as a result, breast cancer
cells are intrinsically multi-target. To model this behavior and generate
cancer-specific models from Cancer RNA-Seq Nexus, an extended flux balance
analysis (FBA) framework was used, coupled with a multi-omic integration
method and multi-level linear programming. Several approaches have been
proposed and reviewed for the integration of gene expression data into FBA
models [22, 23]. Each reaction in a FBA model is controlled by an associated
gene set, defined through AND/OR Boolean operators between genes. In this
work, METRADE [24] was used to constrain the upper- and lower- limits of each
reaction as a function of the expression level of the gene set controlling the
reaction.

Multiple cellular targets were modeled as a trilevel linear program:

max $\displaystyle h^{\intercal}v$ (1) such that $\displaystyle\text{max
}g^{\intercal}v$ such that $\displaystyle\text{max }f^{\intercal}v,\quad
Sv=0,$ $\displaystyle v^{\text{min}}\varphi(\Theta)\leq v\leq
v^{\text{max}}\varphi(\Theta),$

where $S$ is the stoichiometry matrix of the metabolic reactions in the cell,
$v$ is the vector of reaction flux rates, while $f$, $g$, $h$ are Boolean
vectors of weights selecting the three reactions in $v$ whose flux rate will
be considered as objective. Lower- and upper-limits for the flux rates in $v$
in the unconstrained model are given by the vectors $v^{\text{min}}$ and
$v^{\text{max}}$. The vector $\Theta$ represents the gene set expression of
the reactions associated with the fluxes in $v$. The expression level $\Theta$
of a gene set is defined from the expression levels $\theta(g)$ of its genes.
According to the type of gene set, we define $\Theta(g)=\theta(g)$ for single
genes, $\Theta(g_{1}\wedge g_{2})=\min\\{\theta(g_{1}),\theta(g_{2})\\}$ for
enzymatic complexes, and $\Theta(g_{1}\vee
g_{2})=\max\\{\theta(g_{1}),\theta(g_{2})\\}$ for isozymes. These rules were
applied recursively in case of nested gene sets. To enable transcript-
specificity, METRADE was also applied at the splice-isoform level, with
annotations retrieved as above. The function $\varphi$ maps the expression
level of each gene set to a coefficient for the lower- and upper-limits of the
corresponding reaction, and is defined as

$\varphi(\Theta)=\left[\boldsymbol{1}+\gamma\left|log(\Theta)\right|\right]^{\mbox{sgn}(\Theta-\boldsymbol{1})}.$
(2)

Note that the vector notation was adopted, with the convention $0^{0}=1$ when
some element of $\Theta$ is $1$. The sign operator returns a vector of $\pm 1$
(signs of $\Theta-\boldsymbol{1}$). $\gamma$ models the reliability of the
gene set expression level as an indicator of the rate of production of the
associated enzyme.

One may argue that gene expression data is not a good proxy for protein
abundance and ultimately for flux rates. However, gene expression data is
certainly the omic data with better quality and coverage (almost always
genome-scale). Furthermore, in mammals, mRNA levels can be considered the main
contributors to the overall protein expression level [25, 26]. Positive
correlation between mRNA and protein levels was also recently found in most
normal and cancer cell lines [27]. We remark that, although GEMsplice natively
handles absolute expression values, it can also handle fold-change values
directly, by replacing the logarithmic map (2) with
$\varphi(\Theta)=\Theta^{\gamma}$ (see README file in the source code).
Furthermore, our method is general and not dependent on a particular proxy.
For instance, one may use protein abundances within the same approach proposed
here. All simulations were carried out in Matlab R2016b and Octave 4.0.3, with
the GLPK solver.

#### Modeling cancer biomarkers in multi-level linear programming.

The three objectives of the linear program in Equation (1) are chosen to model
known mechanisms and key players in cancer metabolism. Here we consider: (i)
the biomass (growth rate) reaction, (ii) pyruvate kinase (PK) and (iii)
lactate dehydrogenase (LDH). (Note however that this is fully customizable in
GEMsplice.) In our pipeline, we minimize the maximum allowed flux rate of PK
as the second level of the trilevel linear program, selected through the
vector $g$. As a third (outer) level of our linear program, selected by the
vector $h$, we model the Warburg effect by minimizing the flux through LDH. In
fact, a negative flux rate of this reversible reaction models the production
of lactate from pyruvate. The first level (inner level), governed by the
vector $f$, is assigned to represent the maximization of the flux through the
biomass reaction, modeling the observation that cancer cells grow and divide
faster than normal cells in order to achieve their main goal, proliferation.

##  3 Results

### GEMsplice maps RNA-Seq expression levels and splice isoforms onto a
genome-scale model of breast cancer

GEMsplice maps RNA-Seq expression profiles to the metabolic model of breast
cancer, using for the first time splice isoform annotations. A cohort of
$1455$ patients from Cancer RNA-Seq Nexus was taken into account with $31$
different breast cancer subsets, $13$ of which are invasive carcinomas at
different stages. The subset classification is built from observation of
genotype and phenotype, and can for instance refer to a specific type of stage
of cancer, a disease state, or a particular cell line. To assess differences
in the metabolic response among the different types and stages of breast
cancer, we consider the average of RNA-Seq expression levels within the same
cancer subset, and we finally map these to the phenotypic space of flux rates
in the metabolic model using trilevel optimization (Figure 2, Figure
LABEL:sup_fig:biomass_lactate_pyruvate and File S2; see Methods for details on
how phenotype predictions are achieved from RNA-Seq data).

Figure 2: The $31$ RNA-Seq profiles representing different breast cancer
subsets are mapped to the tridimensional space of biomass, pyruvate kinase and
lactate dehydrogenase. The trilevel linear program (1) is solved after
constraining the breast cancer metabolic model using each RNA-Seq expression
profile from Cancer RNA-Seq Nexus, including splice-isoform expression levels.
The bottom left and both right panels show the three projections onto each
pair of axes. In each of these, the optimal breast cancer profiles are
highlighted and connected for both cancer and normal profiles to identify the
trade-off Pareto optimality.

GEMsplice correctly predicts breast cancer cells to grow faster than normal
healthy cells, while keeping comparable levels of PK activity. Some breast
cancer cells are also predicted to use less PK compared to the healthy
counterpart, leading to accumulation of intermediates in the glycolytic
pathway, therefore fueling the serine biosynthesis pathway [28]. The trade-off
between high and low PK flux might also explain the balance between PKM2
(limiting step, controlling PK flux) and PKM1 (promoting OXPHOS in the
mitochondrion instead). The trade-off between LDH and PK in cancer cells is
also lower than in normal cells in terms of negative flux rate values,
consistent with the widely accepted evidence that cancer cells use the reverse
direction of LDH to produce lactate from pyruvate (Warburg effect). When all
the $31$ breast cancer subsets are considered, the bottom right panel shows a
weak negative correlation between biomass and negative flux of lactate
dehydrogenase (Pearson’s $r=-0.34$, $p\mbox{-value}=7.60\cdot 10^{-2}$,
Spearman’s $\rho=-0.37$, $p\mbox{-value}=5.58\cdot 10^{-2}$). These
predictions are in keeping with the widely accepted fact that the amount of
lactate production is positively correlated with tumor growth. More
specifically, in breast cancer, a positive correlation has been highlighted
between degree of malignancy, degree of mitochondrial structural abnormality
[29], and the most common biomarkers of malignant tumor, e.g. intensive use of
glycolysis and increased production of lactate [30].

### Splice isoform expression-based flux control analysis

To further analyze the cancer-specific models, and to assess the role of the
expression of splice isoforms in the model, we propose a steady-state control
analysis based on transcript specificity. This method builds on the technique
named metabolic control analysis [31], which considers the relationship
between enzyme activity and flux rates. We adapt this analysis to assess the
contribution of each splice isoform to the cellular goals.

As described above, flux rates in our model also depend on the expression
level of splice isoforms. In each sample, a control analysis on flux rates
$v_{i}$ with respect to these expression levels $x_{j}$ involves the
estimation of the relation between fractional changes in the flux rates and
fractional changes in the splice isoform expression. This can be written as a
scaled partial derivative of the form $C=\frac{x}{v}\frac{\partial v}{\partial
x}.$

We here adopt an adjusted control coefficient by shrinking the sample-specific
denominator $v$ towards its average $\bar{v}$ across all the samples. As in
the standard deviation correction for differential expression tests [32], our
correction prevents gross overestimation of control coefficients. Without this
correction, “false positives” would arise when $v_{i}(x_{j})$ is a very small
value, irrespective of the numerator. More formally, to approximate this
partial derivative, we separately take into account each splice isoform
$j=1,...,N$ included in the breast cancer metabolic model, and we evaluate
positive and a negative intermediate control coefficients $C_{ij}^{+}$ and
$C_{ij}^{-}$, defined as

$C_{ij}^{\pm}=\frac{v_{i}(x_{j}\pm\delta)-v_{i}(x_{j})}{(v_{i}(x_{j})+\bar{v_{i}})/2}\bigg{/}\frac{\delta}{x_{j}},\
i=1,...,M,\ j=1,...,N,\\\ $ (3)

where $\delta$ is small enough so that the ratio approximates the derivative;
$v_{i}$ is the flux of interest (in our analysis, we focus on the three
objectives of the linear program, namely biomass, PK and LDH). Then, for each
flux and splice isoform of interest, we finally calculate the overall flux
control coefficient $C_{ij}$ as the maximum variation caused by a positive or
negative perturbation of the isoform expression level:

$C_{ij}=\text{max}(\left|C_{ij}^{+}\right|,\left|C_{ij}^{-}\right|).$ (4)

These transcript- and cancer-specific coefficients evaluate the relative
steady-state change in three pivotal flux rates in breast cancer cells, with
respect to the relative change in the expression level of the transcript (see
File S3). For our simulations, we set the denominator $x_{j}=1$ to avoid
proportionality between the final control coefficient and the expression value
itself, and $\delta=10^{-3}$. In general, the smaller the value chosen for
$\delta$, the more effective $C$ in approximating the scaled derivative of
$v_{i}$. Note that calculating a flux control coefficient for each splice
isoform does not explicitly identify controlling transcripts, but rather
provides an effective way to quantify their influence on the key flux rates in
breast cancer cells.

The $10$ most influential transcripts were selected independently for biomass,
PK and LDH. The union of these three sets of ten transcripts, composed of $14$
transcripts, was chosen to perform further enrichment analysis through PANTHER
[33]. A functional classification of each transcript was obtained using the
“protein class” ontology, which is adapted from the PANTHER/X molecular
function ontology and also includes Gene Ontology (GO) annotations.
Interestingly, the selected transcripts are highly enriched for transmembrane
transport and respiratory electron transport chain.

As shown in Figure 3a, the most effective transcript at controlling the value
of biomass and LDH with a less stringent control on PK is ENST00000330775, a
transcript of glucose-6-phosphate translocase (G6PT). In brain cancer, G6PT is
a key player in transducing intracellular signaling events; modulating its
expression has been proposed as an anti-cancer strategy [34]. Furthermore, our
results suggest that PK and LDH are maximally and simultaneously controlled by
ENST00000591899 and ENST00000378667, transcripts of ubiquinol-cytochrome c
reductase, whose amplification was suggested to correlate with more aggressive
breast cancer [35]. ENST00000507754 and ENST00000327772, transcripts of
Complex I (whose activity is known to regulate breast cancer progression
[36]), can control both PK and LDH without disrupting the growth rate.
Although further experimental investigation is needed on these key
transcripts, our genome-scale method may suggest for the first time
transcript-specific targets for anti-cancer strategies.

### Pathway-based flux analysis

A pathway-based perspective has been often taken in genome-scale models with
the aim of investigating sensitivity analysis [37, 38, 39], and coupled with
Bayesian techniques to detect pathway crosstalks and temporal activation
profiles [40]. To assess the variation in the average flux of each pathway
with respect to the unconstrained breast cancer model, we here compute a
normalized average pathway flux

$d_{i}=\left({\bar{w}^{(i)}-w_{U}^{(i)}}\right)/{w_{U}^{(i)}},\quad
i=1,...,P,$ (5)

where $\bar{w}^{(i)}$ is the average flux in the $i$th pathway across the
different cancer subsets, while $w_{U}^{(i)}$ indicates the flux of the $i$th
pathway in the unconstrained breast cancer model. To account for the tolerance
of the linear solver, average pathway fluxes of less than $10^{-10}$ were
assumed to be zero. Pathways, defined using the KEGG LIGAND database, were
inherited from Recon1.

Figure 3b shows the indicator $d$ plotted for each pathway when computed
across breast cancer cells and normal cells (results for invasive and
unlabeled breast cells are shown in Figure LABEL:sup_fig:bubble_pathways). The
least-square linear regression reveals a positive pathway flux correlation
(Pearson’s $r=0.43$, $p\mbox{-value}=1.82\cdot 10^{-3}$, Spearman’s
$\rho=0.80$, $p\mbox{-value}=2.99\cdot 10^{-12}$). The eleven pathways
indicated in the figure were detected as outliers, with significantly
different behavior between cancer and normal breast cells.

To support our predictions, literature-based evidence of breast cancer
alterations in such pathways was aggregated. Alanine and aspartate are
byproducts of amino acid and glutamine fermentation, an important source for
protein synthesis, especially when the cell lacks oxygen. Aspartate is also a
main contributor to nucleotide and protein synthesis; the malate/aspartate
shuttle pathway translocates electrons for the mitochondrial electron
transport chain to produce ATP [41]. Perturbations in the citric acid cycle
(TCA) were largely expected because ot the known “Warburg” reduction of
glucose metabolism through the mitochondrion. Likewise, the prediction of
altered fatty acid elongation metabolism is in keeping with recent studies
reporting elongation of fatty acids as a marker in breast cancer [42].

GEMsplice also correctly predicts that pyrimidine biosynthesis pathways are
downregulated in normal cells when compared to breast cancer cells [43].
Vitamin A and carbohydrate altered metabolism are also associated with breast
cancer [44, 45]. In breast cancer cells, Inosine monophosphate (IMP)
dehydrogenase inhibitors cause growth reduction and phenotypic alterations
[46]. The prediction of starch and sucrose metabolism confirms previous
genome-wide association studies, which reported significant relation between
altered starch and sucrose pathway and the risk of developing ER-negative
breast cancer [47]. Sphingolipids are responsible for oestrogen signaling and
can control differentiation and proliferation of breast cancer cells [48]. The
glycerophospholipid pathway, also flagged as biomarker by GEMsplice, was
previously found altered in the metabolic phenotype of breast [49]. Finally,
triacylglycerol was found to be a biomarker and a prognostic predictor of poor
clinical outcome in triple negative breast cancer [50].

|  |
---|---|---
(a) | (b) | (c)
Figure 3: Transcript-based metabolic flux and pathway analysis. (a) The
control analysis is applied to the set of $14$ transcripts obtained as a union
of the three sets of ten most influential transcripts with respect to biomass,
pyruvate kinase and lactate dehydrogenase flux rate. Error bars represent the
standard error of the mean $\mbox{\it SE}=\sigma/\sqrt{n}$, where $\sigma$ is
the standard deviation of the measured effect, computed across all breast
cancer subsets, and $n$ represents the number of breast cancer subsets. (b)
The indicator $d$ plotted for each pathway in breast cancer versus normal
breast cells (adjacent to tumor cells). The size of each point represents the
size of the corresponding pathway, quantified as the number of reactions
belonging to it. The $y=x$ line is shown on the plots to highlight outliers.
By definition of $d$, if a pathway lies outside the line, its perturbation in
cancer samples is different from its perturbation in healthy samples. (Both
perturbations are computed with respect to the default breast cell model run
without omic-derived constraints.) Colors are given according to a normalized
error of the computed mean flux in each pathway, namely $\mbox{\it
NE}=\sigma/\sqrt{p}$, where $\sigma$ is the standard deviation of the flux
rate in the pathway, computed across all breast cancer subsets, and $p$
represents the size of the pathway. Note that the average flux in a pathway is
positive or negative depending on the direction in which its reactions take
place. (c) GEMsplice correctly predicts 13 out of 15 pathway
increased/decreased activity (difference “Diff” between upregulated and
downregulated reactions), while GIMME only predicts 7 out of 15. More accurate
predictions are achieved for key cancer biomarkers: glycolysis, glycine,
serine, and threonine metabolism, and valine, leucine and isoleucine
degradation.

### Comparison with gene-level approaches

We here compare the results obtained by GEMsplice with GIMME [51], where
expression levels are mapped to the model with gene-level resolution only.
Present/absent calls for GIMME were obtained by applying zFPKM normalization
[52]. To compare the pathway activity predicted by GEMsplice and GIMME in
cancer conditions, we computed the average absolute value of flux rates in all
nonzero biomass samples. As a benchmark for comparison and validation, we used
the set of metabolic pathways found to be significantly upregulated and
downregulated in unfavorable breast cancer conditions [53]. Pathways with
equal number of upregulated and downregulated reactions were excluded from the
analysis.

GEMsplice correctly predicts $13$ out of $15$ pathway
overexpression/underexpression patterns, while GIMME only predicts $7$ out of
$15$ (Figure 3c, and File S5). Three key pathways widely accepted to be
dysregulated in breast cancer are correctly identified by GEMsplice, but
incorrectly by GIMME: (i) glycolysis, highly upregulated in cancer, shows 3.44
cancer/normal fold change in GEMsplice, while 0.55 in GIMME; (ii) glycine,
serine, and threonine metabolism, upregulated in cancer: 2.54 fold change in
GEMsplice, 0.41 in GIMME; (iii) valine, leucine and isoleucine degradation,
downregulated in cancer, 0.30 fold change in GEMsplice, 2.85 in GIMME.

The exact flux rates predicted for lactate dehydrogenase activity were also
evaluated against the widely-accepted cancer “Warburg” metabotype. Ten cancer
samples are incorrectly predicted by GIMME with zero or negligible
($<10^{-12}$) “Warburg” LDH flux (therefore not switching to aerobic
glycolysis), while only three by GEMsplice. Furthermore, all non-malignant and
normal cells with nonzero growth are correctly predicted by GEMsplice with
lowest activity of LDH, while two out of four nonzero-growth normal or non-
malignant tissues are incorrectly predicted by GIMME with a typical cancer
metabotype, consisting of very high LDH activity. These examples suggest that
integrating splice isoform information with metabolic models improves the
characterization of breast cancer metabolism. With the current fast-paced
improvements in speed and costs of omic profiling, and therefore with
increasing availability of isoform-level model annotations and RNA-Seq data,
we expect such specific differences to increase even further.

##  4 Conclusion

While the rate of acquisition of omics data is rapidly increasing, analyzing
and extracting information through computational tools arguably remains the
main bottleneck in biology. Most studies focus on statistical analysis on gene
expression values to evaluate how they vary across different samples. However,
modeling how the gene expression alterations change metabolic processes at
genome scale provides greater understanding of the phenotypic outcome with
respect to studies involving only transcriptomic data [54]. Given that the
last decade in cancer research has repeatedly shown that cancer is a complex
disease that cannot be studied by narrowing it down to a single gene or
enzyme, a genome-wide metabolic approach seems therefore the right direction
to assess and predict the phenotype of a cancer cell.

We here proposed GEMsplice, a method for linking gene expression and splice
isoform data to genome-scale metabolic models. GEMsplice is the first attempt
at solving one of the main issues of metabolic modeling: as outlined in recent
reviews [9, 7, 8, 5], current methods only allow integration of omics data up
to the gene level, but not with splice-isoform resolution.

The idea is that every single profile can be used to create a profile-specific
model of metabolism that includes splice-isoform annotations. This integrated
model can be readily used to predict the flux rate of any biochemical reaction
included in the metabolic model, in a range of cancer subtypes. Individual
profiles, patients or cells related to the same cancer can be mapped onto this
model, for instance to cross-compare their metabolic behavior within the same
cancer type, instead of cross-comparing cells using transcriptomics data only.

Although it allows mapping cancer-specific transcript expression levels onto
any metabolic network, GEMsplice comes with limitations and room for
improvement. For instance, if reaction-specific information is available, the
map from genes to reaction flux bounds could be chosen in a reaction-specific
manner and tailored to specific environmental or growth conditions. Likewise,
post-transcriptional regulation of expression levels could be readily included
in the model if available. On the other hand, inconsistencies can be used to
improve the model and the associated biological knowledge on metabolic
enzymes. For instance, incorrect predictions of specific flux rates may shed
light on where post-transcriptional regulation take place [13]. A further
development of cancer metabolic studies may involve a model that takes into
account interactions between cancer cells and their environment, e.g.
interactions with supporting cells. This approach will likely capture features
that cannot be observed with models of single cancer cells.

Breast cancer is manifested through a variety of effects that cannot be
reduced to a single feature. The difference in behavior of different cancer
cells shows that reconstructing a generic model of a cancer cell is not a
viable approach [55]. Models should therefore always be created in a tissue-
and stage-specific fashion. We anticipate that GEMsplice will allow for the
first time the generation of such models harnessing the full potential of RNA-
Seq, and will facilitate in silico combinatorial experiments with RNA-Seq
data. In fact, such omic-informed models can now effectively integrate
information on the expression level of splice isoforms, to date largely
ignored.

## Acknowledgements

The author would like to thank Dr Jim Liu for discussions on Cancer RNA-Seq
Nexus, and Dr Syed Haider for discussions on splice isoforms.

## References

  * [1] Zachary D Stephens, Skylar Y Lee, Faraz Faghri, Roy H Campbell, Chengxiang Zhai, Miles J Efron, Ravishankar Iyer, Michael C Schatz, Saurabh Sinha, and Gene E Robinson.  Big data: astronomical or genomical?  PLoS Biol, 13(7):e1002195, 2015.
  * [2] Zhen Qi and Eberhard O Voit.  Inference of cancer mechanisms through computational systems analysis.  Molecular BioSystems, 13(3):489–497, 2017.
  * [3] Qun Pan, Ofer Shai, Leo J Lee, Brendan J Frey, and Benjamin J Blencowe.  Deep surveying of alternative splicing complexity in the human transcriptome by high-throughput sequencing.  Nature genetics, 40(12):1413–1415, 2008.
  * [4] Natalie C Duarte, Scott A Becker, Neema Jamshidi, Ines Thiele, Monica L Mo, Thuy D Vo, Rohith Srivas, and Bernhard Ø Palsson.  Global reconstruction of the human metabolic network based on genomic and bibliomic data.  Proceedings of the National Academy of Sciences, 104(6):1777–1782, 2007.
  * [5] Jae Yong Ryu, Hyun Uk Kim, and Sang Yup Lee.  Reconstruction of genome-scale human metabolic models using omics data.  Integrative Biology, 7(8):859–868, 2015.
  * [6] Nicholas Wong, Diane Ojo, Judy Yan, and Damu Tang.  Pkm2 contributes to cancer metabolism.  Cancer letters, 356(2):184–191, 2015.
  * [7] Thomas Pfau, Maria Pires Pacheco, and Thomas Sauter.  Towards improved genome-scale metabolic network reconstructions: unification, transcript specificity and beyond.  Briefings in bioinformatics, 17(6):1060–1069, 2016.
  * [8] Keren Yizhak, Barbara Chaneton, Eyal Gottlieb, and Eytan Ruppin.  Modeling cancer metabolism on a genome scale.  Molecular systems biology, 11(6):817, 2015.
  * [9] Jun Geng and Jens Nielsen.  In silico analysis of human metabolism: reconstruction, contextualization and application of genome-scale models.  Current Opinion in Systems Biology, 2:28–37, 2017.
  * [10] Jan Schellenberger, Richard Que, Ronan MT Fleming, Ines Thiele, Jeffrey D Orth, Adam M Feist, Daniel C Zielinski, Aarash Bordbar, Nathan E Lewis, Sorena Rahmanian, et al.  Quantitative prediction of cellular metabolism with constraint-based models: the COBRA Toolbox v2.0.  Nature protocols, 6(9):1290–1307, 2011.
  * [11] Livnat Jerby, Lior Wolf, Carsten Denkert, Gideon Y Stein, Mika Hilvo, Matej Oresic, Tamar Geiger, and Eytan Ruppin.  Metabolic associations of reduced proliferation and oxidative stress in advanced breast cancer.  Cancer research, 72(22):5712–5720, 2012.
  * [12] JR Li, CH Sun, W Li, RF Chao, CC Huang, XJ Zhou, and CC Liu.  Cancer RNA-Seq Nexus: a database of phenotype-specific transcriptome profiling in cancer cells.  Nucleic acids research, 44(D1):D944–51, 2016.
  * [13] Elke Katrin Markert and Alexei Vazquez.  Mathematical models of cancer metabolism.  Cancer & metabolism, 3(1):1–13, 2015.
  * [14] Otto Warburg, Franz Wind, and Erwin Negelein.  The metabolism of tumors in the body.  The Journal of general physiology, 8(6):519–530, 1927.
  * [15] Adam D Richardson, Chen Yang, Andrei Osterman, and Jeffrey W Smith.  Central carbon metabolism in the progression of mammary carcinoma.  Breast cancer research and treatment, 110(2):297–307, 2008.
  * [16] Lifang Yang, Sriram Venneti, and Deepak Nagrath.  Glutaminolysis: A hallmark of cancer metabolism.  Annual Review of Biomedical Engineering, 19(1), 2017.
  * [17] Ben Langmead, Cole Trapnell, Mihai Pop, and Steven L Salzberg.  Ultrafast and memory-efficient alignment of short dna sequences to the human genome.  Genome biology, 10(3):R25, 2009.
  * [18] Jennifer Harrow, Adam Frankish, Jose M Gonzalez, Electra Tapanari, Mark Diekhans, Felix Kokocinski, Bronwen L Aken, Daniel Barrell, Amonida Zadissa, Stephen Searle, et al.  Gencode: the reference human genome annotation for the encode project.  Genome research, 22(9):1760–1774, 2012.
  * [19] Adam Roberts and Lior Pachter.  Streaming fragment assignment for real-time analysis of sequencing experiments.  Nature methods, 10(1):71–73, 2013.
  * [20] Bo Li and Colin N Dewey.  RSEM: accurate transcript quantification from rna-seq data with or without a reference genome.  BMC bioinformatics, 12(1):323, 2011.
  * [21] Damian Smedley, Syed Haider, Steffen Durinck, Luca Pandini, Paolo Provero, James Allen, Olivier Arnaiz, Mohammad Hamza Awedh, Richard Baldock, Giulia Barbiera, et al.  The biomart community portal: an innovative alternative to large, centralized data repositories.  Nucleic Acids Research, 43(W1):W589–W598, 2015.
  * [22] Daniel Machado and Markus Herrgård.  Systematic evaluation of methods for integration of transcriptomic data into constraint-based models of metabolism.  PLoS Comput. Biol, 10:e1003580, 2014.
  * [23] Supreeta Vijayakumar, Max Conway, Pietro Lió, and Claudio Angione.  Seeing the wood for the trees: a forest of methods for optimization and omic-network integration in metabolic modelling.  Briefings in Bioinformatics, bbx053, 2017.
  * [24] C Angione and P Lió.  Predictive analytics of environmental adaptability in multi-omic network models.  Scientific Reports, 5:15147, 2015.
  * [25] Jingyi Jessica Li, Peter J Bickel, and Mark D Biggin.  System wide analyses have underestimated protein abundances and the importance of transcription in mammals.  PeerJ, 2:e270, 2014.
  * [26] Marko Jovanovic, Michael S Rooney, Philipp Mertins, Dariusz Przybylski, Nicolas Chevrier, Rahul Satija, Edwin H Rodriguez, Alexander P Fields, Schraga Schwartz, Raktima Raychowdhury, et al.  Dynamic profiling of the protein life cycle in response to pathogens.  Science, 347(6226):1259038, 2015.
  * [27] Idit Kosti, Nishant Jain, Dvir Aran, Atul J Butte, and Marina Sirota.  Cross-tissue analysis of gene and protein expression in normal and cancer tissues.  Scientific reports, 6, 2016.
  * [28] Jason W Locasale, Alexandra R Grassian, Tamar Melman, Costas A Lyssiotis, Katherine R Mattaini, Adam J Bass, Gregory Heffron, Christian M Metallo, Taru Muranen, Hadar Sharfi, et al.  Phosphoglycerate dehydrogenase diverts glycolytic flux and contributes to oncogenesis.  Nature genetics, 43(9):869–874, 2011.
  * [29] RL Elliott, XP Jiang, and JF Head.  Mitochondria organelle transplantation: introduction of normal epithelial mitochondria into human cancer cells inhibits proliferation and increases drug sensitivity.  Breast cancer research and treatment, 136(2):347–354, 2012.
  * [30] Michael J Gonzalez, Jorge R Miranda Massari, Jorge Duconge, Neil H Riordan, Thomas Ichim, Ana I Quintero-Del-Rio, and Norma Ortiz.  The bio-energetic theory of carcinogenesis.  Medical hypotheses, 79(4):433–439, 2012.
  * [31] H Kacser, JA Burns, and DA Fell.  The control of flux.  Biochemical Society Transactions, 23(2):341–366, 1995.
  * [32] Virginia Goss Tusher, Robert Tibshirani, and Gilbert Chu.  Significance analysis of microarrays applied to the ionizing radiation response.  Proceedings of the National Academy of Sciences, 98(9):5116–5121, 2001.
  * [33] Huaiyu Mi, Sagar Poudel, Anushya Muruganujan, John T Casagrande, and Paul D Thomas.  Panther version 10: expanded protein families and functions, and analysis tools.  Nucleic acids research, 44(D1):D336–D342, 2016.
  * [34] Anissa Belkaid, Ian B Copland, Duna Massillon, and Borhane Annabi.  Silencing of the human microsomal glucose-6-phosphate translocase induces glioma cell death: potential new anticancer target for curcumin.  FEBS letters, 580(15):3746–3752, 2006.
  * [35] Yoko Ohashi, Saori J Kaneko, Tommy E Cupples, and S Robert Young.  Ubiquinol cytochrome c reductase (uqcrfs1) gene amplification in primary breast cancer core biopsy samples.  Gynecologic oncology, 93(1):54–58, 2004.
  * [36] Antonio F Santidrian, Akemi Matsuno-Yagi, Melissa Ritland, Byoung B Seo, Sarah E LeBoeuf, Laurie J Gay, Takao Yagi, and Brunhilde Felding-Habermann.  Mitochondrial complex i activity and nad+/nadh balance regulate breast cancer progression.  The Journal of clinical investigation, 123(3):1068, 2013.
  * [37] Edward Kent, Stefan Neumann, Ursula Kummer, and Pedro Mendes.  What can we learn from global sensitivity analysis of biochemical systems?  PloS one, 8(11):e79244, 2013.
  * [38] J. Costanza, G. Carapezza, C. Angione, P. Lió, and G. Nicosia.  Robust design of microbial strains.  Bioinformatics, 28(23):3097–3104, 2012.
  * [39] Max Conway, Claudio Angione, and Pietro Lió.  Iterative multi level calibration of metabolic networks.  Current Bioinformatics, 11(1):93–105, 2016.
  * [40] Claudio Angione, Naruemon Pratanwanich, and Pietro Lió.  A hybrid of metabolic flux analysis and Bayesian factor modeling for multi-omics temporal pathway activation.  ACS Synthetic Biology, 4(8):880–889, 2015.
  * [41] Claire Pecqueur, Lisa Oliver, Kristell Oizel, Lisenn Lalier, and François M Vallette.  Targeting metabolism to induce cell death in cancer cells and cancer stem cells.  International journal of cell biology, 2013, 2013.
  * [42] Yin-Hsun Feng, Wei-Yu Chen, Yu-Hsuan Kuo, Chao-Ling Tung, Chao-Jung Tsao, Ai-Li Shiau, and Chao-Liang Wu.  Elovl6 is a poor prognostic predictor in breast cancer.  Oncology Letters, 12(1):207–212, 2016.
  * [43] Frederic D Sigoillot, Severine M Sigoillot, and Hedeel I Guy.  Breakdown of the regulatory control of pyrimidine biosynthesis in human breast cancer cells.  International journal of cancer, 109(4):491–498, 2004.
  * [44] Elena Doldo, Gaetana Costanza, Sara Agostinelli, Chiara Tarquini, Amedeo Ferlosio, Gaetano Arcuri, Daniela Passeri, Maria Giovanna Scioli, and Augusto Orlandi.  Vitamin a, cancer treatment and prevention: the new role of cellular retinol binding proteins.  BioMed research international, 2015, 2015.
  * [45] Arvindhan Nagarajan, Parmanand Malvi, and Narendra Wajapeyee.  Oncogene-directed alterations in cancer cell metabolism.  Trends in Cancer, 2(7):365–377, 2016.
  * [46] Y Sidi, C Panet, L Wasserman, A Cyjon, A Novogrodsky, and J Nordenberg.  Growth inhibition and induction of phenotypic alterations in mcf-7 breast cancer cells by an imp dehydrogenase inhibitor.  British journal of cancer, 58(1):61, 1988.
  * [47] Jingmei Li, Keith Humphreys, Hatef Darabi, Gustaf Rosin, Ulf Hannelius, Tuomas Heikkinen, Kristiina Aittomäki, Carl Blomqvist, Paul DP Pharoah, Alison M Dunning, et al.  A genome-wide association scan on estrogen receptor-negative breast cancer.  Breast Cancer Research, 12(6):R93, 2010.
  * [48] O Sukocheva and Carol Wadham.  Role of sphingolipids in oestrogen signalling in breast cancer cells: an update.  Journal of Endocrinology, 220(3):R25–R35, 2014.
  * [49] Cristina Cadenas, Sonja Vosbeck, Eva-Maria Hein, Birte Hellwig, Alice Langer, Heiko Hayen, Dennis Franckenstein, Bettina Büttner, Seddik Hammad, Rosemarie Marchan, et al.  Glycerophospholipid profile in oncogene-induced senescence.  Biochimica et Biophysica Acta (BBA)-Molecular and Cell Biology of Lipids, 1821(9):1256–1268, 2012.
  * [50] Danian Dai, Bo Chen, Bin Wang, Hailin Tang, Xing Li, Zhiping Zhao, Xuan Li, Xiaoming Xie, and Weidong Wei.  Pretreatment tg/hdl-c ratio is superior to triacylglycerol level as an independent prognostic factor for the survival of triple negative breast cancer patients.  Journal of Cancer, 7(12):1747, 2016.
  * [51] Scott A Becker and Bernhard Ø Palsson.  Context-specific metabolic networks are consistent with experiments.  PLoS computational biology, 4(5):e1000082, 2008.
  * [52] Traver Hart, H Kiyomi Komori, Sarah LaMere, Katie Podshivalova, and Daniel R Salomon.  Finding the active genes in deep rna-seq gene expression studies.  BMC genomics, 14(1):778, 2013.
  * [53] Gunnar Schramm, Eva-Maria Surmann, Stefan Wiesberg, Marcus Oswald, Gerhard Reinelt, Roland Eils, and Rainer König.  Analyzing the regulation of metabolic pathways in human breast cancer.  BMC medical genomics, 3(1):39, 2010.
  * [54] Claudio Angione, Max Conway, and Pietro Lió.  Multiplex methods provide effective integration of multi-omic data in genome-scale models.  BMC Bioinformatics, 17(4):257, 2016.
  * [55] Pouyan Ghaffari, Adil Mardinoglu, Anna Asplund, Saeed Shoaie, Caroline Kampf, Mathias Uhlen, and Jens Nielsen.  Identifying anti-growth factors for human cancer cell lines through genome-scale metabolic modeling.  Scientific reports, 5, 2015.


\begin@twocolumnfalse

This paper embarks upon exploration of quantum vulnerability analysis. By introducing vulnerability graphs, related to attack graphs, this paper provides background theory and a subsequent method for solving significant cybersecurity problems with quantum computing. The example given is to prioritize patches by expressing the connectivity of various vulnerabilities on a network with a QUBO and then solving this with quantum annealing. Such a solution is then proved to remove all kill-chains (paths to security compromise) on a network. The results demonstrate that the quantum computer's solve time is almost constant compared to the exponential increase in classical solve time for vulnerability graphs of expected real world density. As such, this paper presents a novel example of advantageous quantum vulnerability analysis.

\end@twocolumnfalse]

§ INTRODUCTION

Much of the defensive work carried out by organisations in the cybersecurity space has become a data problem [17]. Whilst this presents several challenges, it also presents numerous opportunities for moving the operating mode of any Security Operations Centre (SOC) from mostly reactive to including predictive activities [3].

Quantum computing has heralded several major claims for the future of computing, particularly surrounding the significant optimisations (theoretically) to be found in quantum computing on specific problem classes, including number factoring, parameter discovery, travelling salesman problems, and many more [13].

This paper presents a method for utilising such quantum advantage in a new approach to the theoretical underpinnings of network-based cybersecurity. By finding a new way of solving some of the data problems in cyber-defence, and improving the response times thereby, the aim is to demonstrate a more efficient move towards prioritising and applying vulnerability patches.

Patch management is a common pain point for any large scaled enterprises [19] or widely distributed systems such as smartphones or IoT devices [5]. Indeed, lack of appropriate patching was indicated as being a central cause for some high profile cybersecurity breaches, such as the infamous Equifax hack [12]. A variety of approaches have been proposed to improve the categorization and management of patches, including deep learning technologies [21].

This paper demonstrates this by analysing vulnerability data on hosts as a bipartite graph. With this, we reason that attacks are made up of `kill-chains', which are themselves comprised of sequences of exploits leveraging vulnerabilities that are coincident - in our model, by being shared on a host. By then creating a process through which we can totally disconnect vulnerabilities from one another, we effectively remove every possible kill-chain. This problem, however, involves a known NP-hard problem, but one which is tractable on a DWave quantum computer (see section <ref>). The full results of various tests of this method can be found in section <ref>, and the code may be found in [2].

By leveraging quantum computation and optimisation methods for vulnerability analysis of this kind, new avenues of optimisation of cybersecurity and related data can be considered. The example presented here is fully worked through to demonstrate where a possible quantum advantage may lie for this kind of analysis on this type of data.

The leveraging of short-time solutions to NP-hard problems that are present in cybersecurity data is a potentially rich vein of exciting possibilities. The fast and efficient resolution of cybersecurity data problems also helps reduce the analysis and reaction times of security teams, thereby tightening their OODA loop and increasing their security level proportionately [16].

§ PREREQUISITES

§.§ Graph Theoretic Pre-requisites

[4]
A graph $G = (V,E)$; $V(G)$ a set of vertices and a set of edges $E(G) \subseteq V \times V$, $E$ is composed of pairs of elements from $V$.

All of the graphs we will consider here are finite simple graphs; the number of edges and vertices is finite, and the graph has only one unique edge between any two vertices, with no edges starting and terminating on the same vertex (simple).

[4]
A path $P$ is a sequence of edges $e_0, e_1, \ldots, e_n \in E$ such that for every $e_{i-1}, e_i, e_{i+1} \in P$, $\exists v_i \in (e_{i-1} \cap e_i) \text{ and } \exists v_{i+1} \in (e_{i} \cap e_{i+1})$. There are two distinguished vertices, one at the start and one at the end, and a path $P$ is a loop iff the starting and terminating edges in the path overlap.

[4]
A weighted graph $G$ is a graph with edges that are triples $(u,v,w)$ such that $u,v \in V$ and weight $w \in \mathbb{N}$. If $w$ is not specified, it is assumed to be 1.

The use of the weight parameter $w$ is such that we can talk about the maximum/minimum `flow' across different paths along a graph.

We shall be using a bipartite graph structure, defined as follows:

[4]
A bipartite graph $G$ is a graph with a partition of $V(G)$ into two sets $A, B$ such that $\forall (a,b) \in E(G)$, $a \in A$ and $b \in B$.

§.§ Graph Algorithms

Various algorithms have been found that provide graph structured information with significant insight into the underlying relationships between different elements in a dataset. The graph algorithm utilised in this paper concerns finding the minimum vertex cover.

A vertex cover for a graph $G$ is a subset $C \subset V(G)$ such that for every edge $(u,v) \in E(G)$ at least one of $u$ or $v$ is in $C$. $C$ is minimal if $V(G) \setminus C$ is as large as possible.

The idea behind a vertex cover is that it is a set of vertices that `touches' every edge in the graph. That is for any vertex cover $C$ for a graph $G$, $V(G) \setminus C \implies E(G) = \emptyset $. An important feature of vertex covers is that finding one is a known NP-hard problem.[8]

Let an independent set $\mathcal{I}$ be a set of vertices on a graph such that no two members of $\mathcal{I}$ are connected. An alternative definition that we use in the code in [2] is that a minimum vector cover of some graph $G$ is the compliment of a maximal independent set $\mathcal{I} \subseteq V(G)$.

§.§ Quantum Annealing and QUBOs for Graph Algorithms

There are many problems that have been found to be easily coded into Quadratic Unconstrained Binary Optimisation (QUBO) problems[Also called `Unconstrained Binary Quadratic Programs', or UBQPs. See [9].][10]. Indeed, DeSimone et al. [18] showed how Ising Hamiltonians can be equivalent to graph problems, an idea that we make use of here.

[18]
A Quadratic Unconstrained Binary Optimisation (QUBO) is the problem of finding, for $x \in \mathbb{B}^n$ the following
$$ \min x^\top Q x$$
for a given upper triangular matrix $Q \in \mathbb{R}^{n \times n}$.

DeSimone et al. [18] showed that QUBO problems can be represented in the form
$$ f_Q(x) = \sum_{i=1}^n \sum_{j=1}^i q_{ij} x_i x_j $$
for coefficients $q_{ij} \in \mathbb{R}$, $x_i, x_j \in \mathbb{B}^n$. This is equivalent to extracting the coefficients from $Q$ in definition <ref>. This is equivalent to Ising Hamiltonians of the form
$$ H(\sigma) = -\sum_{\langle i,j \rangle} J_{ij} \sigma_i \sigma_j - \mu \sum_{j} \sigma_0 \sigma_j $$
for real parameters $J_{ij}$, $\sigma_i$, and $\mu$.

With this equivalence, there is a way to use the adiabatic quantum computing methods of DWave Systems to acquire the vector values for $x$ above such that it can solve the minimum vector cover problem we described in section <ref> [15].

§ NEW APPROACHES TO VULNERABILITY PATCH PRIORITISATION

This section describes a graph theoretic way of analysing vulnerabilities on a computer network called `vulnerability graphs', derived from the notion of `attack graphs' [11]. The aim of this section is to present a new way to prioritise the patching of vulnerabilities by considering their connectivity and solving for this using quantum computation.

Our aim here is to provide a way of identifying the `most well connected' issues by means of describing the dual on our vulnerability graphs.

§.§ Attack Graphs and Kill-chains

`Attack graphs' have featured in some interesting approaches to managing and mitigating security threats ([6], [1], [11]). They provide various ways regarding how to analyse network-oriented vulnerability data that many cybersecurity information sources generate. As pointed out in François-Xavier et al. [6], representations of a Bayesian Attack Model can take the form of network graphs, from which Matthews et al. [11] show that various network characteristics can be extracted. Callahan eg al. go so far as to construct and optimisation model for virtualised networks that incorporate attack graphs to give some guarantee of security.

Whilst the formulation in [11] utilises directed graphs, this paper uses un-directed simple graphs to represent the same data. This is due to the fact that it is not that important to consider directionality for the purposes of this analysis.

A `kill chain' is a multi-stage sequence of events that leads to the compromise of a network [23]. What is apparent is that many of the examples of kill-chains involve sequences of vulnerabilities, with the sequence dependant on the assets that are intersected between these vulnerabilities.

§.§ Vulnerability Graphs

Define a vulnerability graph as follows:

A vulnerability graph $G$ is a bipartite graph where one partition of vertices represents network hosts, and the other represents vulnerabilities. Let the edges of $G$ represent that a given host is affected by some detected vulnerability.

This data structure is a very direct way of representing the results from industry standard vulnerability scanning and penetration test reports [17]. This particular data structure lends itself to analysing the data in vulnerability reports much more efficiently. This is due to the fact that vulnerability reports are asset-first, i.e. they are a list of hosts with sub-lists of what vulnerabilities affect them. However, many of the questions we wish to ask of this dataset are the other way round - we wish for a list of vulnerabilities with sub-lists of affected hosts. To search host-first data to extract all vulnerabilities, it is clear that all the host records must be read. This is computationally expensive, compared with a graph structure.

For this paper, we shall consider kill-chains as sub-sequences of paths through the vulnerability graph. For our purposes, we shall use the following definition:

A kill chain is a sequence of vertices $K = \{ v_1, v_2, \ldots v_n \}$ from the vulnerability partition of a vulnerability graph $\mathcal{V}$ such that for each $v_i, v_j \in K$, there exists at least one host $h \in \mathcal{V}$ with $(v_i,h), (v_j,h) \in E(\mathcal{V})$.

Each part of a kill chain that comprises of movement from one vulnerability to the next will start on some vulnerability vertex, go to some host vertex, and then on to another vulnerability vertex connected to that same host.

This paper shall only consider kill-chains that require vulnerabilities to be on the same host, not adjacent hosts on a given network subnet. However, it is very easy to consider adding nodes to the host partition of the vulnerability graph to encode these visibility relationships for vulnerabilities that are network facing. But to save space, the approach here will not consider this for now.

It is worth pointing out that the lack of any information coded about severity ratings for vulnerabilities, e.g. CVSS scores. This information is not considered here as critical vulnerabilities should always be patched as soon as possible. The aim of these definitions is to consider non-critical paths to compromise as generally a subset of all paths through a vulnerability graph, and solve for these issues by analysing their connectivity.

In short, this process aims to find the lower criticality issues that are widespread and well connected enough to cause more harm later.

It is important not to only consider only the vulnerabilities that we know to be `very bad'. This approach, in line with the quantum speedup we describe later, allows cybersecurity to consider an entire attack surface, not just potentially isolated pain points. This is the logical extension of the common maxim “defence in depth" [7].

§.§ Connectivity Dual Graphs and Vector Covers

We define the following dual graph construction for a vulnerability graph $\mathcal{V}$ with a partition of vulnerability vertices and host vertices:

The dual $\mathcal{D}_{\mathcal{V}}$ is constructed as follows. For each vulnerability vertex $v_i \in V(\mathcal{V})$ for $1 \leq i \leq |V(\mathcal{V})|$:

* Add $v_i$ to $V(\mathcal{D}_{\mathcal{V}})$ if $v_i \notin V(\mathcal{D}_{\mathcal{V}})$.

* Enumerate a list of host vertices $\{ h_0 , h_1 , \ldots \}$ connected to $v_i$.

* Iterating over this list of hosts, for each $v'_i$ connected to each host $h_j$:

* Add $v'_i$ to $V(\mathcal{D}_{\mathcal{V}})$.

* Add $(v_i, v'_i)$ to $E(\mathcal{D}_{\mathcal{V}})$

* If $(v_i, v'_i)$ already exists, add 1 to the weight of that edge.

* Remove $v_i$ from $V(\mathcal{V})$ and continue with $v_{i+1}$.

The dual $\mathcal{D}_{\mathcal{V}}$ represents all of the connections between vulnerabilities on a network. It is weighted, so that priority can be given to vulnerabilities that are more connected than others, by virtue of the weight coding the number of hosts a vulnerability was found on.

§.§ Removing Kill-chains with Vertex Covers

Previous approaches, notably [11], have relied on locating cycles of probabilities in attack graphs. The idea in this paper is to utilise vulnerability graphs to consolidate all possible avenues for compromise, the kill-chains, and then remove all possible attack routes through a network.

Removing the vertices in a vertex cover on $\mathcal{D}_{\mathcal{V}}$ from $V(\mathcal{V})$ will leave $\mathcal{V}$ totally disconnected on the vulnerability partition to itself via the host partition.

Edges on the dual graph $\mathcal{D}_{\mathcal{V}}$ represent an `edge-host-edge' sub-path between one vulnerability and another on $\mathcal{V}$. The vertex cover on $\mathcal{D}_{\mathcal{V}}$ therefore intersects every `vulnerability-host-vulnerability' path on $\mathcal{V}$.

Let $v_i, v_j$ be vulnerabilities and $h_k$ a host in $V(\mathcal{V})$. Removing every host in the vertex cover of $\mathcal{D}_{\mathcal{V}}$ will result in a path of the form $\{(v_i,h_k),(h_k,v_j)\}$ being removed, and the sequence $\{v_i, h_k, v_j\}$ reduced to just one of $\{v_i, h_k\}$ or $\{h_k,v_j\}$. Therefore, when all the vertices are removed, $\mathcal{V}$ will be disconnected from the vulnerability side, as there is no way to get from any vulnerability to any other vulnerability.

The point of this proof is to show that in order to remove all possible kill-chains, we do not need to resolve all vulnerabilities.

A disconnected vulnerability graph can contain no kill-chains.

Consider our definition <ref>. By removing every `vulnerability-host-vulnerability' sub-path in a vulnerability graph $\mathcal{V}$ by means of a minimum vertex cover on $\mathcal{D}_{\mathcal{V}}$, we have removed every possible kill chain $K$ found in the paths of $\mathcal{V}$.

Based on our assumption of a kill chain's structure, this approach will remove every kill chain. Furthermore, it becomes clear that we do not need to patch everything in order to have a significant impact and improvement. We can patch smarter, not harder.

§.§ Considerations

It has been a mainstay for some time that many notions of security rely, in some significant sense, on the reaction times of those tasked with defending assets [16]. Simply put; the shorter our detection and reaction times and the faster we iterate over security data with feedback from new information, the better our security posture will be.

The problem with considering attack graphs at enterprise scale is the sheer size of the datasets that could be involved. Thus, if we can minimise time taken to analyse large vulnerability datasets then we can improve enterprise reaction times to these issues being identified.

By utilising an adiabatic quantum computing setup, this time should be shortened as we have solved an NP-hard problem efficiently. By adding in other work by Pelofske et al [15] we can reduce the time required to find the most at-risk vulnerabilities and patch them with more priority.

However, as pointed out in Di Tizio et al [20], much of the activity in patch management may be relatively superfluous for any attack that is below the complexity of a nation state APT level threat actor.

§.§ Weighted MVC

We recall that in the construction process of the dual in section <ref>. As pointed out in [8], the weighted minimum vector cover (wMVC) problem is reducible to the non-weighted MVC. As shown in Pelofske et al [15], these can also be performed on adiabatic quantum computers by means of coding the problem into a QUBO.

Therefore, by applying a weighted MVC to our weighted dual $\mathcal{D}_{\mathcal{V}}$, we can improve our output by being able to prioritise more highly connected vulnerabilities over less well connected ones. Although this does not change any of the proofs above, it is a significant improvement our analysis.

§ A WORKED EXAMPLE

We present the following vulnerability graph $\mathcal{V}$, with hosts $a$ to $g$, and vulnerabilities $1$ through $8$:

[nodeDecorate/.style=inner sep=2pt,draw,thick,minimum size=0.4cm,

1/0/0/blue/circle, 2/2/0/blue/circle, 3/4/0/blue/circle, 4/6/0/blue/circle, 5/8/0/blue/circle, 6/10/0/blue/circle, 7/12/0/blue/circle, 8/14/0/blue/circle,
a/1/2.5/red/rectangle, b/3/2.5/red/rectangle, c/5/2.5/red/rectangle, d/7/2.5/red/rectangle, e/9/2.5/red/rectangle, f/11/2.5/red/rectangle, g/13/2.5/red/rectangle

() at (,) [nodeDecorate] [label=left: ,fill=,̧shape=] ;

1/a, 1/b, 1/d, 1/f,
2/a, 2/b,
3/a, 3/d, 3/e,
4/b, 4/c, 4/f,
5/g, 6/f, 6/g,
7/b, 7/c, 7/f,
8/c, 8/d, 8/g

() edge[lineDecorate] node  ()

The host connections are given by the following vulnerability to host edges, summarised by:

* $1$ — $\{ a, b, d, f \}$

* $2$ — $\{ a, b \}$

* $3$ — $\{ a, d, e \}$

* $4$ — $\{ b, c, f \}$

* $5$ — $\{ g \}$

* $6$ — $\{ f, g \}$

* $7$ — $\{ b, c, f \}$

* $8$ — $\{ c, d, g \}$

From this graph we can get the following dual graph $\mathcal{D}_{\mathcal{V}}$:

[nodeDecorate/.style=inner sep=2pt,draw,thick,minimum size=0.5cm,shape=circle,

1/2/0, 2/4/0, 3/6/2, 4/6/4, 5/4/6, 6/2/6, 7/0/4, 8/0/2

() at (,) [nodeDecorate] $\name$;

() edge[lineDecorate] node  ()

If we compute the MVC we get the set $\{ 1, 2, 4, 8, 6 \}$. Removing these nodes from our original vulnerability graph $\mathcal{V}$ we get:

* $3$ — $\{ a, d, e \}$

* $5$ — $\{ g \}$

* $7$ — $\{ b, c, f \}$

Or as a diagram:

[nodeDecorate/.style=shape=circle,inner sep=2pt,draw,thick,minimum size=0.4cm,

3/8/0/blue/circle, 5/12/0/blue/circle,  7/4/0/blue/circle,
a/0/2.5/red/rectangle, b/2/2.5/red/rectangle, c/4/2.5/red/rectangle, d/6/2.5/red/rectangle, e/8/2.5/red/rectangle, f/10/2.5/red/rectangle, g/12/2.5/red/rectangle

() at (,) [nodeDecorate] [label=left: ,fill=,̧shape=] ;

3/a, 3/d, 3/e, 5/g, 7/b, 7/c, 7/f

() edge[lineDecorate] node  ()

This demonstrates that our algorithm and theorem is correct, as well as illustrating the density which can arise in vulnerability graphs and the duals we defined for them.

§ SOLVING WITH QUANTUM HARDWARE

We now demonstrate our approach using a DWave AQC system.

§.§ Implementation

To test the above process on quantum hardware, it was decided to make use of the built in functions within DWave's implementation of the  python library for graph programming. The back-end that was used was the `Advantage Solver 4.1' QPU system.

For benchmarking, the same algorithms were run through the accompanying  for QUBO's that solves using classical methods. We implemented our own solver, based on DWave's provided code, for finding independent sets with DWave's python interface. This gave us more control over the generation of solutions and checking. The code is available on our github: see [2].

To explore the hardware capabilities, python code was written to perform 3 solves on the duals automatically generated for increasingly larger random bipartite graphs with varying edge probabilities, timing the results on both classical and quantum solvers [2]. The given MVC candidates are then checked that all `vuln-host-vuln' paths have been eliminated.

2c|[HTML]EFEFEFc]@c@ Number of vulnerabilities
Disclosed to CVE Databases
1|l|[HTML]EFEFEFSoftware Packages     1r|[HTML]EFEFEF2020                                        All Time (as of 2021)
1|l|Windows 10                                             1r|807                      2990
1|l|Windows Server 2016                                    1r|794                      2764
1|l|Linux Kernel (all versions)                            1r|126                      3000
1|l|Debian Linux                                           1r|946                      7331
1|l|Ubuntu Linux                                           1r|483                      3680
1|l|Mac OSX (all versions)                                 1r|314                      3100
1|l|Android (all versions)                                 1r|859                      4707
1|l|iOS (for iPhone)                                       1r|322                      2820
1|l|Chrome Browser                                         1r|227                      2554
1|l|Firefox Browser                                        1r|141                      1993
1|l|Microsoft Office Suite                                 1r|71                       727

This table shows the number of disclosed vulnerabilities to the CVE database [14] for some select popular software packages, for 2020 specifically and `All Time'.

§.§ Limitations and Real World Patching Estimates

Owing to the lack of publicly available vulnerability data from scans that consisted of more than a handful of hosts, it was not possible to validate this approach using real world data.

However, one can estimate the possible occurrence of patches for a given host based on the number of vulnerabilities disclosed for major software operating systems, components, and suites. We use [14] as our primary resource for this data, a summary of some examples is in table <ref>.

For any vulnerability to enter the CVE database [14], it is usually accompanied by some kind of coordinated disclosure, for which a patch is generally also made available. As such, the number of publicly disclosed vulnerabilities should correlate reasonably accurately to the number of available patches for that software.

Considering our chart in table <ref> we find that anywhere from 7-27% of all vulnerabilities for a given software package were disclosed in 2020 alone. We can estimate the average annual patch inflation rate (that is, growth in the number of patches year on year) to be somewhere between 5-15%.

Although there is a lack of any hard evidence for how many clients and servers any given enterprise commissions, maintains, or decomissions in any year, we can estimate from our annual patch inflation rate that for a given 5 year period the number of patches will grow between 27% and 100% for our patch inflation range.

If we assume largely homogeneous networks - i.e. that most hosts are derived from a small number of `gold builds' to facilitate the recommendations in [19] - then the number of vulnerabilities  will be relatively small compared to the overall number of hosts.

Combining all of these estimates, we find that with an average (from [14]) of $\approx 3000$ vulnerabilities per major software package, an enterprise with 10 dominant software packages will accrue anywhere up to 8,000 new vulnerabilities in 5 years with a patch inflation rate of 5%, each package potentially adding $\approx 800$ vulnerabilities to the CVE database.

For an enterprise with 100,000 hosts (clients and servers) this would give a minimal guesstimate of $\approx 8\%$ probability that any vulnerability is connected to any given host, assuming an even distribution of software. This is the estimate that we will use later in our experimentation with a DWave AQC, but should be considered somewhat of an upper bound.

Clearly time really is of the essence when applying patches, and so any process to speed this activity up has significant potential utility given the potential growth rates involved.

§.§ Results and Analysis

The results of running the outlined procedure in section <ref> on three random graphs of varying edge probabilities, ranging in sizes 16 to 48 vertices (8 to 24 vulnerabilities). Edge probabilities determine the likelihood that any two vertices have an edge between them. Times for the quantum annealer are shown by the orange line, with the blue line indicating the time taken to find the MVC clasically, in microseconds (1e7 scale). The red dots indicate the number of invalid vector cover solutions found by the DWave annealer. The lower right graph shows the averaged difference in sizes between the classical method for finding MVCs (denoted cMVC), and the DWave AQC (deonted qMVC).

4c|50% Edge Probability
4c|33.34% Edge Probability
4c|10% Edge Probability
1|l|Num. Vulns     1l|c]@l@Mean Diff of
VC sizes     1l|cMVC ($\mu s$)     1l|qMVC ($\mu s$)     c]@l@Invalid
VCs     1l|c]@l@Mean Diff of
VC sizes     1l|cMVC ($\mu s$)     1l|qMVC ($\mu s$)     c]@l@Invalid
VCs     1l|c]@l@Mean Diff of
VC sizes     1l|cMVC ($\mu s$)     1l|qMVC ($\mu s$)     c]@l@Invalid
1|l|8              1l|-0.333          1l|2847               1l|16435              1      1l|0               1l|3677               1l|16435              0      1l|0               1l|1911               1l|16982              0
1|l|9              1l|0.333           1l|3450               1l|15579              0      1l|0               1l|2485               1l|17851              0      1l|0               1l|2734               1l|17395              0
1|l|10             1l|-0.667          1l|4876               1l|16451              3      1l|0.667           1l|3812               1l|17547              0      1l|0               1l|3385               1l|16150              0
1|l|11             1l|-0.667          1l|8156               1l|17095              2      1l|0               1l|7738               1l|17464              0      1l|0               1l|8022               1l|15646              0
1|l|12             1l|0.333           1l|12256              1l|16765              2      1l|-0.333          1l|11649              1l|16582              1      1l|0               1l|11038              1l|16214              0
1|l|13             1l|0.667           1l|22217              1l|16413              0      1l|1               1l|20335              1l|16505              0      1l|0               1l|18499              1l|16114              0
1|l|14             1l|0               1l|36199              1l|17154              2      1l|0               1l|33898              1l|36734              0      1l|0               1l|29909              1l|16129              0
1|l|15             1l|-0.333          1l|60240              1l|17238              2      1l|0.333           1l|57864              1l|17940              1      1l|0.333           1l|57489              1l|16195              0
1|l|16             1l|-0.667          1l|89672              1l|16380              2      1l|0.333           1l|89678              1l|17403              1      1l|0               1l|85638              1l|16473              0
1|l|17             1l|-0.667          1l|174828             1l|15852              3      1l|0.333           1l|137885             1l|17143              2      1l|1               1l|128818             1l|16615              0
1|l|18             1l|-0.333          1l|222740             1l|16653              2      1l|0.667           1l|225831             1l|16651              2      1l|0               1l|208862             1l|16617              0
1|l|19             1l|-1              1l|479130             1l|17296              3      1l|0.667           1l|472617             1l|17702              2      1l|0.667           1l|388965             1l|16940              0
1|l|20             1l|0               1l|950424             1l|16194              1      1l|0               1l|945096             1l|18166              2      1l|-0.333          1l|802647             1l|16615              1
1|l|21             1l|-0.667          1l|1928877            1l|17399              2      1l|-0.333          1l|1912908            1l|17339              3      1l|0.667           1l|1590244            1l|16001              0
1|l|22             1l|-0.667          1l|3949102            1l|18930              3      1l|0               1l|3950946            1l|18667              3      1l|0.333           1l|3219998            1l|16636              0
1|l|23             1l|-1.667          1l|8381824            1l|18085              3      1l|0.667           1l|8257182            1l|16495              3      1l|0.667           1l|6928612            1l|16599              0
1|l|24             1l|-1              1l|17086013           1l|16866              3      1l|-0.333          1l|18049568           1l|16306              3      1l|-0.333          1l|14240710           1l|17202              2

This table shows the results for comparing classical and quantum efforts undertaken to find the MVCs. Times are given in microseconds ($\mu s$). The mean difference in size between classical and quantum is gieven by $|qMVC| - |cMVC|$. Also provided are the number of vector covers generated by the DWave AQC that were found to be invalid - i.e. they did not provide valid vector covers that solve the problem.

We see in figure <ref> the results of our process. The benchmark testing was done on a Mac M1 with 8Gb RAM vs. a DWave QPU. The quantitative results are found in table <ref>.

What can be observed is that whilst the classical algorithm for finding vector covers increases in time taken proportional exponentially to the size of the graph, in line with what one would expect for an NP-hard problem, the quantum algorithm processing remains essentially constant, increasing only slightly in line with the number of variables presented in the QUBO - a number proportional to the number of edges in our dual graphs $\mathcal{D}_{\mathcal{V}}$, and so proportional secondarily to the size and density of the starting bipartite graphs.

What should be noted here is that as the randomly generated graphs move away from the worst-case edge density (that of around one half) the improvements in fidelity on the DWave system improves significantly. In fact, the DWave is producing comparable Vector Cover results for large graphs with 10% edge density, but in almost strictly static time. Indeed, the inflection point for when our problem test set takes longer to compute classically seems to be around 12 vertices, although this will likely change with larger and larger graphs.

For comparison, the MVC solution on the Dual graph for the random vulnerability graphs with edge probability 50% generated a QUBO problem with 78 variables when sent to the DWave system. For the 10% edge probability vulnerability graphs the QUBO problems had just 33 variables.

The difference in the classical exact and quantum vector cover solution set sizes tells us how close the DWave system came to the proper solution. When looking at the $\Delta$ for each of the values, we can observe; For 50% edge probability graphs, the DWave was missing vertices, giving a negative $\Delta$ value. Whilst for 33.34% and 10% edge probability graphs, the values were on average higher, meaning that the DWave solution had too many vertices.

Even though the DWave provided, on average, too many vertices in a vector cover for it to be minimal, this is in line with classical approximation methods [22]. It should be pointed out that in these results there only seems to be at most 1 extra vertex being provided, which is not a significant overhead.

As such, what we are seeing is that the quantum solver is at least bounded above by the classical solver's time, and would appear to be exponentially faster at solving this problem type at the larger node counts but with edge densities we expect to see in the real world, than its classical counterpart.

§ CONCLUSIONS

In conclusion, this work has found that there is a potential viable use for quantum computing in the field of vulnerability analysis, specifically the prioritisation of patches by this method. It shows some demonstrable promise of workable quantum advantage for the vulnerability graph densities that we expect to see in the real world.

The theoretical work presented here demonstrates how all possible kill-chains can be remediated through careful analysis of the connectivity of vulnerabilities. The experimental data matches the theoretical work, ultimately showing that these problems are made tractable by quantum hardware.

The proof for theorem <ref> shows that existing best practice controls, such as network segregation and zero trust [7], are valid and useful in this model. Indeed, this model validates that possible kill-chains are limited from a vulnerability chaining point of view by the connectivity reduction presented by these controls.

§.§ Further Work

To facilitate future development of these ideas, we have included all of the code for performing the experiments and data analysis in this paper on the Quantum Village Github [2].

There are likely other parallelization options for solving these problems with GPUs that were not explored here at all. There is work, e.g. in [25], that indicates this is a feasible route for future comparison.

The most obvious next step in this work is to remedy the lack of real data to operate on. However have estimated by considering the homogeneity of many modern enterprises [24] we might estimate that there are comparatively few vulnerabilities compared to number of hosts - that is, most vulnerabilities affect most hosts on a highly host-homogeneous network. How this affects the process presented here would be interesting to measure.

To develop this, it should be possible to add more layers to our bipartite graph, making it $n$-partite for $n$ data sources we wish to consider. Kill-chains here would be informed by vulnerability data, host data, host-to-host connectivity data (captured in, say, a subnet layer), threat intelligence data, and more. In a significantly more complex vulnerability graph like this, it might just be found that the quantum computing processing baseline witnessed here provides even more of an advantage.

There is also the question of prioritising the vulnerabilities identified through this MVC method - although, the weighted case and/or a subsequent analysis of centrality (e.g. degree centrality or betweenness centrality, see [4]) will likely provide these.

[1]
Devon Callahan, Timothy Curry, Hazel Davidson, Heytem Zitoun, Benjamin Fuller,
and Laurent Michel.
Fashion: Functional and attack graph secured hybrid optimization of
virtualized networks, 2019.

[2]
Mark Carney.
Quantum Vulnerability Analysis Notebooks.
<https://github.com/Quantum-Village/Quantum-Vuln-Analysis>, 11

[3]
Clarence Chio and David Freeman.
Machine Learning and Security.
O'Reilly Media, Sebastopol, CA, February 2018.

[4]
Reinhard Diestel.
Graph Theory.
Graduate texts in mathematics. Springer, Berlin, Germany, 4 edition,
July 2010.

[5]
Sadegh Farhang, Mehmet Bahadir Kirdan, Aron Laszka, and Jens Grossklags.
Hey google, what exactly do your security patches tell us? a
large-scale empirical study on android patched vulnerabilities, 2019.

[6]
Aguessy François-Xavier, Bettan Olivier, Blanc Grégory, Conan Vania, and
Debar Hervé.
Bayesian attack model for dynamic risk assessment, 2016.

[7]
Stephen Groat, Joseph Tront, and Randy Marchany.
Advancing the defense in depth model.
In 2012 7th International Conference on System of Systems
Engineering (SoSE), pages 285–290, 2012.

[8]
Richard M. Karp.
Reducibility among combinatorial problems.
In Complexity of Computer Computations, pages 85–103. Springer
US, 1972.

[9]
Gary Kochenberger, Jin-Kao Hao, Fred Glover, Mark Lewis, Zhipeng Lü, Haibo
Wang, and Yang Wang.
The unconstrained binary quadratic programming problem: a survey.
Journal of Combinatorial Optimization, 28(1):58–81, April

[10]
Andrew Lucas.
Ising formulations of many NP problems.
Frontiers in Physics, 2, 2014.

[11]
Isaac Matthews, John Mace, Sadegh Soudjani, and Aad van Moorsel.
Cyclic bayesian attack graphs: A systematic computational approach.
In 2020 IEEE 19th International Conference on Trust, Security
and Privacy in Computing and Communications (TrustCom). IEEE, December

[12]
Lily Hay Newman.
Equifax officially has no excuse, Sep 2017.

[13]
Michael A Nielsen and Isaac L Chuang.
Quantum Computation and Quantum Information.
Cambridge University Press, Cambridge, England, December 2010.

[14]
Serkan Özkan.
Top 50 products by total number of "distinct" vulnerabilities.
<https://www.cvedetails.com/top-50-products.php>, 2022.

[15]
Elijah Pelofske, Georg Hahn, and Hristo Djidjev.
Solving large minimum vertex cover problems on a quantum annealer.
In Proceedings of the 16th ACM International Conference on
Computing Frontiers. ACM, April 2019.

[16]
Winn Schwartau.
Analogue Network Security.
Schwartauhaus, 04 2018.

[17]
Jim Scott.
Cybersecurity is a data problem.
Sep 2021.

[18]
C. De Simone, M. Diehl, M. Jünger, P. Mutzel, G. Reinelt, and G. Rinaldi.
Exact ground states of ising spin glasses: New experimental results
with a branch-and-cut algorithm.
Journal of Statistical Physics, 80(1-2):487–496, July 1995.

[19]
Murugiah Souppaya.
Guide to enterprise patch management planning: Preventive maintenance
for technology.
Technical report, National Institute of Standards and Technology,

[20]
Giorgio Di Tizio, Michele Armellini, and Fabio Massacci.
Software updates strategies: a quantitative evaluation against
advanced persistent threats.
IEEE Transactions on Software Engineering, pages 1–1, 2022.

[21]
Xinda Wang, Shu Wang, Pengbin Feng, Kun Sun, Sushil Jajodia, Sanae Benchaaboun,
and Frank Geck.
PatchRNN: A deep learning-based system for security patch
In MILCOM 2021 - 2021 IEEE Military Communications
Conference (MILCOM). IEEE, November 2021.

[22]
David P. Williamson and David Bernard Shmoys.
The design of approximation algorithms.
Cambridge University Press, New York, 2011.
OCLC: ocn671709856.

[23]
Tarun Yadav and Arvind Mallari Rao.
Technical aspects of cyber kill chain.
In Communications in Computer and Information Science, pages
438–452. Springer International Publishing, 2015.

[24]
Qisheng Zhang, Abdullah Zubair Mohammed, Zelin Wan, Jin-Hee Cho, and
Terrence J. Moore.
Diversity-by-design for dependable and secure cyber-physical systems:
A survey.
IEEE Transactions on Network and Service Management,
19(1):706–728, 2022.

[25]
Wenyong Zhong, Jianhua Sun, Hao Chen, Jun Xiao, Zhiwen Chen, Chang Cheng, and
Xuanhua Shi.
Optimizing graph processing on gpus.
IEEE Transactions on Parallel and Distributed Systems,
28(4):1149–1162, 2017.