DOI HERE Advance Access Publication Date: Day Month Year Paper

Liuyi et al.

[$\ast$]Xuekui Zhang<EMAIL_ADDRESS>and Li<EMAIL_ADDRESS>
0Year 0Year 0Year

# BOSS - Biomarker Optimal Segmentation System

Liuyi Lan  Xuanjin Cheng  Li Xing 0000-0002-4186-7909  Xuekui Zhang
0000-0003-4728-2343  Department of Mathematics and Statistics, University of
Victoria, Victoria, BC, Canada  Department of Medical Genetics, The University
of British Columbia, Vancouver, BC, Canada  Department of Mathematics and
Statistics, University of Saskatchewan, Saskatoon, SK, Canada

(2023; 2023; Date; Date; Date)

###### Abstract

Motivation: Precision medicine is a major trend in the future of medicine. It
aims to provide tailored medical treatment and prevention strategies based on
an individual’s unique characteristics and needs. Biomarker is the primary
source of patients’ unique features used in precision medicine. We often need
to investigate many cutoff values of a continuous biomarker to find the
optimal one and test if it can help segment patients into two groups with
significantly different clinical outcomes. This requires multiple testing
adjustment on tests conducted on overlapped data. The permutation-based
approach is often a preferred solution, since it does not suffer the
limitations of state-of-art theoretical methods. However, permutation is
computationally expensive and limits its application scenarios, such as web
applications requiring a fast response or the analysis of genomic study
requiring to repeat analysis many times on tens of thousands of genes.
Results: We proposed a novel method BOSS, Biomarker Optimal Segmentation
System, to solve this problem. In simulation studies, we found BOSS’s
statistical power and type I error control are both non-inferior to the
permutation approach, and it is hundreds of times faster than permutation. To
illustrate our method, we applied BOSS to real data and revealed potentially
converging biomarkers that have referential importance in exploring synergy
and target-matched therapies in lung adenocarcinoma.
Availability: An R package boss is being developed and will be available on
https://cran.r-project.org/.

###### keywords:

multiple testing adjustment, patient segmentation, optimal cutoff, data
overlapping structures

##  1 Introduction

Precision medicine, the medical approach of employing patients’ specific
information to tailor medical treatments to satisfy their needs, has become
increasingly popular recently. It is partially due to the contribution of
technological advances, enabling individuals’ genomic data collection at a
fast speed with reduced cost (add a reference). And these genomic data serve
as the essential and major source for biomarkers in precision medicine
(Ashina, ; Vargas, ). One of the most common practices of implementing
precision medicine in health care and management is to use biomarkers to
stratify patients into subgroups according to their associated risks of a
disease or their likelihood of positive responses to treatments or
interventions. Different treatments are applied to the corresponding subgroups
based on such stratification to alleviate risks or improve positive responses.
In practice, we often seek a cutoff value of a continuous biomarker to segment
patients into two subgroups labelled as high-risk vs low-risk groups or
responder vs non-responder and assign different treatments to the two groups.
To achieve this goal, we must investigate multiple candidate cutoff values and
answer two closely related questions, (1) which cutoff is the optimal one, and
(2) whether the outcome variables differ significantly between these two
groups segmented using the biomarker. And the answers can be revealed based on
testing a series of candidate cutoffs to check their significance and compare
their importance based on the effects carried on the test statistics. In this
situation, we must apply multiple testing adjustments to correct the inflated
false positive probabilities due to conducting multiple testings.

The popular multiple testing adjustment methods control the false discovery
rate (FDR) or the familywise error rate (FWER). FDR is defined as the expected
proportion of false rejections among the total rejections, considering the
seriousness of the loss incurred by erroneous rejections to be inversely
related to the number of hypotheses rejected (Benjamini, ). FWER is defined as
the probability of making any false rejections. When a small number of false
positives do not affect the model’s overall effectiveness, controlling FDR is
often preferred since it is less stringent and leads to more statistical power
than FWER. So FDR is much more commonly used when multiple testing adjustments
are needed. However, the FDR adjustment method requires a large number of
hypothesis testing to estimate the FDR precisely. In the problem of finding
the optimal cutoff, the more candidate cutoff values investigated, the heavier
multiple testing penalties are introduced, leading to more severe lack-of-
power issues. Thus, we often cannot investigate enough candidate cutoff values
to estimate FDR reliably. Hence, we choose to control FWER in this problem.

Bonferroni correction (Bonferroni, ) is the most popular FWER control method.
To control FWER below a threshold $\alpha$, Bonferroni correction uses
$\alpha/n$ as the significance level for all tests, where $n$ is the number of
hypotheses to be tested. As it assumes all hypotheses are independent,
Bonferroni correction is often over-conservative and reduces statistical
power. For example, when we repeatedly conduct the ‘same’ test on the ‘same’
dataset twice, we will get the same p-value twice, and we should not apply any
multiple testing adjustments to this p-value even though two tests are carried
out. The rule can be further generalized. The penalty applied to multiple
testing should be relied on the correlations between these tests. Less penalty
should be considered under high correlations among the test statistics.
Different candidate cutoff converts the continuous biomarker into different
binary variables in the segmentation problem. Any two derived binary variables
overlap, and patients whose biomarker levels are above or under both cutoffs
have exactly the same 0/1 values in the derived variables. Consequently, the
corresponding test statistics obtained from these derived variables must be
correlated. The closer the values of the two cutoffs, the higher the
correlation between their test statistics. Therefore, Bonferroni correction is
unsuitable for our problem due to its independence assumption. Worsley
(Worsley, ) proposed an improved Bonferroni, which orders the datasets and
explicitly models correlations between test statistics obtained from every two
‘adjacent’ datasets. However, this method assumes that test statistics
obtained from non-adjacent datasets are independent, which is still
unrealistic in our problem.

Meanwhile, advanced methods are developed using the maximally selected
statistic and its asymptotic distribution to select the optimal cutoffs and
adjust for multiple tests. For binary outcomes, Miller and Siegmund MSC
proposed to use a maximally selected chi-squared statistic to seek optimal
cutoff and Boulesteix and StroblMSCS  proved that using standard chi-square
distribution can render more significance in coefficient estimation and
proposed to obtain true critical values by deriving the asymptotic Gaussian
distribution of maximally chi-square statistic. For continuous and time-to-
event outcomes, Lausen and SchumacherLausen1992  suggested the standardized
maximally selected rank statistic. However, these methods are only applicable
under a large number of candidate cutoffs due to the requirement of asymptotic
properties. As mentioned before, there cannot be many cutoffs in our case. To
adequately address these issues with fewer candidate cutoffs, Hothorn and
ZeileisSTA  proposed the permutation method, which randomly shuffles part of
the data to estimate the empirical distribution of the maximally selected
statistic test under the null hypothesis instead of employing their asymptotic
distribution. Therefore, their method is still valid for the nonlarge number
of candidate cutoffs. Hilsenbeck and ClarkHILSENBECK  compared multiple
testing adjustments based on the improved Bonferroni and the use of asymptotic
distribution and empirical permutation of the maximally selected statistic via
simulation studies and showed that the rejection rates from the permutation
method are most close to those of the independent validation sample.
Hilsenbeck et al. Susan  also noted that permutation and resample can better
identify prognostic factors than other methods regarding accuracy and
efficiency. Because of the advantages of the permutation methods, Cheng et al.
cSurv  applied the permutation-based multiple-testing adjustment and built a
web application providing a pipeline to segment cancer patients based on the
optimal cutoffs of the gene expression biomarkers. However, the website
suffers a low-speed issue due to the heavy computational load of permutation.
It responds around $20$ seconds to a response, which is far from the good load
time (e.g. $0-4$ seconds) recommended for the best website conversion rate
(portent, ). In reality, since we must generate and analyze a lot of copies of
permuted data to achieve a decent precision in the p-value, we cannot go
around the computational issue. And this drawback led to its inefficiency in
practical application.

The drawbacks of existing methods motivate us to develop a novel method, which
does not require investigating many candidate cutoffs to ensure the validity
of asymptotic distribution and is much less computationally expensive than
permutation approaches. Our novel approach borrow the idea from the literature
on Group Sequential Design (GSD) of clinical trials. In the GSD, a series of
interim analyses are conducted on datasets overlapped in a nested way, i.e. in
any two interim analyses, the early dataset is a subset of the latter. The GSD
textbook (group, ) (Page 2 of Chapter 11) introduced that the estimated
coefficients from the same regression model fitted to the nested data follow
the multivariate normal distribution, and their correlations can be explicitly
calculated. The discussion covers various outcome types under their associated
models, including linear and generalized linear regressions in Chapter 3 and
Cox regressions in Chapter 13. Moreover, Lausen and Schumacher Lausen1992
demonstrated that the log-rank statistics for survival outcome also follow the
multivariate normal distribution under random (or non-informative) censoring.
Note that these theoretical results in GSD are based on nested data. In the
patient segmentation problem, the generated binary variables are overlapped
but not in a nested way. This work will use theoretical derivations to extend
GSD results to solve this patient segmentation problem, and propose a novel
method for patient segmentation called BOSS (Biomarker Optimal Segmentation
System).

In BOSS, we model the joint distribution of the test statistics obtained from
all cutoff values as a multivariate normal distribution, utilize the specific
data overlapping structures of different cutoff values to calculate their
correlation structures and precisely calculate the FWER from its definition
without assuming any independence among the tests. Consequently, BOSS has two
major advantages. First, it directly calculates the FWER from a closed-form
formula, leading to a much faster speed than the permutation-based methods.
Therefore, it is suitable to be incorporated into the web system, such as
cSurvival (cSurv, ). Second, since it does not require many candidate cutoffs,
BOSS can ensure decent statistical power, especially when the sample size is
not big enough to investigate a lot of cutoffs.

The rest of this article is organized as follows. Section 2 sets up the
problem and presents our model in detail. Section 3 conducts simulation
studies to evaluate the performance of the BOSS method by comparing it with
the permutation test. Section 4 illustrates our method via its application to
real data. Section 5 includes conclusions and discussions.

##  2 Methods

###  2.1 Notation and Model

Let vector $\mathbf{Y}=(Y_{1},\ldots,Y_{n})$ denote the clinical outcomes of
$n$ patients and vector $\mathbf{B}=(B_{1},\ldots,B_{n})$ denote the
measurements of a quantitative biomarker, where $n$ is the sample size. Let
the $n$ by $(p+1)$ dimension matrix
$\mathbf{C}=(\mathbf{1},\mathbf{C}_{1},\ldots,\mathbf{C}_{p})$ denote the
design matrix consisting of a column of $1$’s and other columns of $p$
covariate values (e.g. age, gender, and BMI, etc). We aim to segment patients
into two groups based on a cutoff value for the biomarker. If the difference
between these two groups’ clinical outcomes is tested statistically
significant, we claim the biomarker can be used to segment these patients. We
investigate $k$ candidate cutoffs $(\tau_{1},\ldots,\tau_{k})$ to find the
best one and test if it can segment patients into two groups with
significantly different clinical outcomes after adjusting for the covariates.
The $i$-th cutoff $\tau_{i}$ converts continuous biomarker $\mathbf{B}$ into a
binary variable $\mathbf{X}_{i}=(X_{i1},\ldots,X_{in})$ defined as

$\displaystyle X_{is}=\begin{cases}1,&\mathrm{if}\;B_{s}>\tau_{i}\\\
0,&\mathrm{if}\;B_{s}\leq\tau_{i}.\end{cases}\;\;\mathrm{for}\;s=1,\ldots,n$
(1)

For each $i$, we consider a linear regression model as below.

$\displaystyle\mbox{E}(\mathbf{Y})=\beta_{i}\mathbf{X}_{i}+\gamma_{0i}+\Sigma_{l=1}^{p}\gamma_{li}\mathbf{C}_{l},$
(2)

where $\beta_{i}$ represents the difference of clinical outcomes between two
groups adjusted to other covariates, $\gamma_{0i}$ is the intercept, and
$\gamma_{li}$ is the regression coefficient for the $l$th covariate.

###  2.2 Hypothesis Testing

Based on the above model, we test the hypothesis of whether the biomarker can
segment patients. The null hypothesis is the biomarker cannot segment patients
at any given cutoff, i.e.

$H_{0}:\;\;\beta_{1}=\beta_{2}=\cdots=\beta_{k}=0.$ (3)

The alternative hypothesis, $H_{1}$, is that at least one of these
$\beta_{i}$s’ is not equal to $0$. We consider the Z statistics, which follows
standard normal distribution under $H_{0}$, i.e.

$\displaystyle Z_{i}=\hat{\beta}_{i}/sd(\hat{\beta}_{i})\sim
N(0,1)\;\;\;\mathrm{under}\;\;H_{0},$ (4)

where $\hat{\beta}_{i}$ is the least squares estimate of the coefficient and
$sd(\hat{\beta}_{i})$ is its standard deviation. If we derive $k$ p-values
from these standard normal distributions, they cannot be directly used to make
a decision, since we must correct the inflated false positive probability from
conducting multiple tests.

Following the literature about maximumly selected statistics, we use the most
extreme test statistics to make decisions, and show it is equivalent to FWER
of conducting $k$ tests in our situation. We denote $z_{i}$ as the observed
value of the Z-statistic, which can be obtained from the output of the fitted
regression models, and will be used to calculate FWER. Under $H_{0}$, the
joint distribution of $Z_{i}$’s is symmetric to zero in any dimension, so we
only need to consider the absolute value of $z_{i}$’s under $H_{0}$. That is,
when $z_{i}<0$, we have

$\displaystyle P(Z_{i}>z_{i})=P(Z_{i}<|z_{i}|)\;\;\;\textrm{when $z_{i}<0$ and
$H_{0i}$ is true.}$

So, to simplify the derivation of FWER, we assume $z_{i}\geq 0$ for all
$i=1,\ldots,k$; and we define the index of optimal cutoff $k^{*}$ as the one
corresponding to the largest observed test statistics,

$\displaystyle k^{*}=\mathrm{argmax}\\{z_{i}|i=1,\ldots,k\\}$ (5)

the type I error of maximumly selected test statistics or the FWER of $k$
tests can be defined as below,

$\displaystyle FWER$ $\displaystyle=1-P(\max_{i}Z_{i}\leq z_{k^{*}}|H_{0})$
$\displaystyle=1-P(Z_{1}\leq z_{k^{*}},\ldots,Z_{k}\leq z_{k^{*}}|H_{0}).$ (6)

We reject the null hypothesis $H_{0}$ if $\mathrm{FWER}<\alpha$, and fail to
reject if $\mathrm{FWER}\geq\alpha$, where $\alpha$ is the significance level
(usually $\alpha=0.05$). To calculate the FWER in the Formula (6), we need
these Z-statistics $Z_{i}$’s joint distribution under the null hypothesis
$H_{0}$, which follows a multivariate normal distribution(group, )

$\displaystyle(Z_{1},\ldots,Z_{k})|H_{0}\sim N(\mathbf{0},\mathbf{\Sigma}).$
(7)

This work’s most challenging part is calculating the elements in the
covariance matrix $\mathbf{\Sigma}$. We give its results below in formulas (8)
and (9), and provide their detailed derivation in the supplementary document.

The $(i,j)$ element of the covariance matrix $\mathbf{\Sigma}$ can be
calculated as

$\displaystyle
Cov(Z_{i},Z_{j})=\frac{\mathbf{X}_{i}^{T}(\mathbf{I}-\mathbf{H})\mathbf{X}_{j}}{\sqrt{\mathbf{X}_{i}^{T}(\mathbf{I}-\mathbf{H})\mathbf{X}_{i}}\sqrt{\mathbf{X}_{j}^{T}(\mathbf{I}-\mathbf{H})\mathbf{X}_{j}}}$
(8)

where $\mathbf{H}=\mathbf{C}(\mathbf{C}^{T}\mathbf{C})^{-1}\mathbf{C}^{T}$ is
hat matrix of covariates. If no covariates exist in the regression model, the
covariance in Formula (8) is reduced to a simpler form

$\displaystyle Cov(Z_{i},Z_{j})=\sqrt{(n-m_{j})m_{i}}/\sqrt{(n-m_{i})m_{j}},$
(9)

where $m_{i}\leq m_{j}$, $m_{i}$ and $m_{j}$ are the numbers of $1$s’ in
vectors $\mathbf{X}_{i}$ and $\mathbf{X}_{j}$ respectively, and $n$ is the
total number of samples. Formula (9) agrees with our intuition that without
covariates in the model, the correlation between the estimated regression
coefficients from models associated with two different cutoffs only depends on
the overlapping $1$’s proportion in these two derived binary variables.

Algorithm 1 Select optimal cutoff and test significance

Clinical Outcome $\mathbf{Y}$; Covariates $\mathbf{C}$; A Continuous Biomarker
$\mathbf{B}$; Candidate Cutoffs $(\tau_{1},\dots,\tau_{k})$; Significance
Level $\alpha=0.05$

1\. For $i$ in $1,\ldots,k$ do

(a) Derive $\mathbf{X}_{i}$ based on the $i$th cutoff $\tau_{i}$ and biomarker
measurements $\mathbf{B}$ using Formula (1).

(b) Regress $\mathbf{X}_{i}$ on $\mathbf{Y}$ with adjustment to $\mathbf{C}$
in a regression model to obtain the estimated regression coefficient
$\hat{\beta}_{i}$, its estimated standard deviation
$\hat{sd}(\hat{\beta}_{i})$, and the Z-statistics $z_{i}$.

end for

2\. Find the index of the optimal cutoff,
$k^{*}=\max_{k}(|z_{1}|,\ldots,|z_{k}|)$.

3\. Calculate the components of $Z_{i}$’s covariance $\mathbf{\Sigma}$, using
Formulas (8) or (9).

4\. Calculate FWER using Formula (6).

5\. If FWER $<\alpha$, reject $H_{0}$ and claim that using biomarker and its
cutoff $\tau_{k^{*}}$ can segment patients into significantly different
groups.

optimal cutoff $\tau_{k^{*}}$ and its FWER value

Finally, we summarize the analysis procedure of BOSS in Algorithm 1.

###  2.3 Extension to other types of regression models

In the GSD textbook (group, ), the correlation formula of Z-statistics among
the nested interim analyses is the same for linear, generalized, and Cox
regression, which shines the light on the belief that the same property should
hold for other types of regression. In this work, we skip the mathematical
derivation for other regression models and use simulation studies to
demonstrate that formula (8) is also valid in Cox regressions.

##  3 Simulation Studies

The main purpose of our simulation studies is to compare the performance of
our novel method, BOSS, with the gold standard permutation method. We simulate
data under various settings and apply both methods to the data to compare
their statistical powers, Type I errors, and computation times. The power
measures the ability of the tool to identify true signals, while the type I
error indicates false discoveries. An optimal tool is expected to have high
power and properly control type I error to a nominated threshold ($0.05$ being
the commonly used value ). The computation time is essential in many
applications, especially for website-based tools and when we must repeat the
analysis many times (e.g. for tens of thousands of genes).

###  3.1 Simulation Procedure

To ensure the simulated data closely mimic the real-world data, we employ real
genomic data to build a blueprint model and simulate outcome variables based
on such blueprint models while reserving genomic information and patients’
other characteristics in the real data. In particular, we use data from a
curated cancer outcome database (cSurv, ), containing expression levels of
$17,431$ genes of $500$ lung adenocarcinoma cases. The detailed background of
the data is illustrated in the next section of the data application. We carry
out simulations in the following two steps.

Step 1: Build Blueprint Model We apply BOSS to the $17,431$ genes separately
in the data mentioned above and obtain their FWERs defined in the Formula (6),
which results in $2,089$ identified significant genes via linear regression
model fitting and $3,502$ identified significant genes via Cox regression
model fitting based on the significance threshold FWER$<0.05$. From the
identified genes, we select the top $50$ genes with the smallest FWERs (i.e.
the top 50 strong-effect biomarkers associated with the outcome) and $50$
genes with the largest FWERs (i.e. the bottom 50 weak-effect biomarkers
associated with the outcome). In total, we obtain $100$ representative genes
for linear regression and Cox regression, respectively. To describe the
blueprint model, we denote the expression levels of these selected genes as
$\mathbf{B}_{g}$ and $\tau^{*}_{g}$ as the $g$-th gene’s optimal cutoff
selected by BOSS with $g=1,\ldots,100$. Let $\mathbf{Y}_{g}$ be the vector of
the outcome variable (quantitative or survival types). For quantitative
outcomes,

$E(\mathbf{Y}_{g})=\beta_{{}_{g}}I\\{\mathbf{B}_{g}>\tau^{*}_{g}\\}+\gamma_{0g},$
(10)

and for survival outcomes,

$h(\mathbf{Y}_{g})=h_{0}(\mathbf{Y}_{g})\mbox{exp}\left(\beta_{{}_{g}}I\\{\mathbf{B}_{g}>\tau^{*}_{g}\\}+\gamma_{0g}\right),$
(11)

where $\beta_{{}_{g}}$ is the association parameter between the $g$th gene
(biomarker) and the outcome variable, $\tau^{*}_{g}$ is the optimal threshold,
and $\gamma_{0g}$ is the corresponding linear component associated with the
covariates ($g=1,...,100$).

In summary, based on $100$ selected genes, we obtained $100$ blueprint models
for the linear regression models and $100$ for Cox regression models.

Step 2: Simulate Data All simulated data retain genomic information from
original real data, and simulate clinical outcomes using the blueprint models.
We simulate quantitative outcomes using the Gaussian random sampler and the
model (10), and simulate survival outcomes using the R package coxed (simsurv,
) and the model (11). This R package uses a random spline method to generate a
baseline hazard function via a cubic spline fitting at randomly drawn points.

For each blueprint model discussed above, we simulate $100$ copies of positive
data and $100$ copies of negative data. We denote the positive data as
$(\mathbf{Y}_{g},\mathbf{B}_{g})$. The outcomes of positive data are generated
using parameters of blueprint model learned from the original data, whose true
association parameters $\beta_{{}_{g}}$s’ are non-zero values. Hence, the
proportion of identified significant biomarkers in data analysis represents
the statistical power. We denote the negative data as
$(\mathbf{Y}^{0}_{g},\mathbf{B}_{g})$. The outcomes of negative data are
generated using zeroed-out association parameters (i.e. setting $\beta_{g}=0$)
in blue print models. Hence, the generated outcomes are unrelated to
biomarkers, and the proportion of identified significant biomarkers in data
analysis represents the type I error.

###  3.2 Analysis Scheme

We apply BOSS to analyze each copy of $(\mathbf{Y}_{g},\mathbf{B}_{g})$ and
$(\mathbf{Y}^{0}_{g},\mathbf{B}_{g})$ separately. We claim the $g$-th gene is
significant, if its BOSS FWER (6) is less than $0.05$. Note that the analyses
are repeated with various values of the number of candidate cutoffs,
$k=6,8,10,12,14$, and we choose the cutoff points $(\tau_{1},\ldots,\tau_{k})$
as equally spaced quantiles of observed biomarker values. In this way, we
obtain $100$ FWERs for every gene from each simulated data under various $k$
and two outcome types.

As for evaluation metrics, given a $g\in\\{1,2,\cdots,100\\}$, a
$k\in\\{6,8,10,12,14\\}$ and one outcome type (i.e. quantitative or survival
outcome), we use the proportion of FWER$<0.05$ from these $100$ copies
$(\mathbf{Y}_{g},\mathbf{B}_{g})$ to represent the statistical power of the
BOSS test and use the proportion of FWER$<0.05$ from these $100$ copies
$(\mathbf{Y}^{0}_{g},\mathbf{B}_{g})$ to represent the Type I error rate of
the BOSS test under such condition.

We repeat the above analyses using a permutation approach with $1000$
permutations to achieve the precision of $0.001$ in each estimated FWER, which
results in the corresponding pairs of estimated statistical power and Type I
error, based on the same rules for BOSS.

###  3.3 Simulation Results

####  3.3.1 Comparisons based on Evaluation Metrics

Figure 1: Boxplots of type I error rate and power from simulated data with
quantitative outcomes. In the top panel, the boxplots of type I error rates of
both methods under various $k$, the number of candidate cutoffs, are
presented. In the middle and bottom panels, the boxplots of the power of both
methods based on the $50$ strong biomarkers and $50$ weak biomarkers under
various $k$ are presented separately. Figure 2: Boxplots of type I error rate
and power from simulated data with survival outcomes. In the top panel, the
boxplots of type I error rates of both methods under various $k$, the number
of candidate cutoffs, are presented. In the middle and bottom panels, the
boxplots of the power of both methods based on the $50$ strong biomarkers and
$50$ weak biomarkers under various $k$ are presented separately.

We employ side-by-side boxplots of the resulting power and type I error to
demonstrate the performance of the BOSS and permutation methods in Figure 1
for quantitative outcomes and Figure 2 for survival outcomes separately.

First, in the top panel of each figure, boxplots of type I errors of the two
methods are presented under various $k$, the number of candidate cutoffs, with
a dashed horizon line indicating the nominated threshold $0.05$ for
controlling type I error rate. These boxplots indicate the type I errors of
both BOSS and the permutation tests are properly controlled around $0.05$
(i.e. centred at $0.05$ with most of IQRs being less than $0.04$). The one-
sample Wilcoxon tests are conducted to numbers in each box to compare these
resulting type I errors against the nominated level of $0.05$. We found none
of these boxes is significantly different from $0.05$ after the Bonferroni
correction. Therefore, both methods perform well in controlling false
positives at the desired level.

Second, the middle panels show the medians of the powers obtained from the
$50$ strong-effect genes of both BOSS and the permutation method are all above
$80$% for quantitative outcomes (Figure 1) and all close to $60$% for survival
outcomes (Figure 2) under various simulation settings, which indicates both
methods have similar decent powers under all the simulation settings for
biomarkers with strong effects. When the biomarkers carry weak signals, the
medians of the powers from both methods decrease to around $30$% for
quantitative outcomes (Figure 1) and $20$% for survival outcomes (Figure 2)
separately.

Based on the visual comparison from boxes of BOSS versus permutation approach,
we found these two methods perform very similarly. We confirmed the visual
pattern using the paired Wilcoxon tests, and found no significant p-value
after the Bonferroni correction. Therefore, we conclude BOSS performs equally
well as the permutation approach regarding statistical power and Type I error
control.

####  3.3.2 Comparisons based on Computation Time

Figure 3:  Summary of computation time of the BOSS and the permutation methods
in simulation studies. The mean and the corresponding one-standard-error bar
of computation time for each method are aggregated based on the pooled
computation time under each $k$, the number of candidate cutoffs, from the
simulation studies.

We summarize the computation time of both methods from simulation studies and
present them in Figure 3. It clearly shows that BOSS is much faster than a
permutation approach. In particular, BOSS decreases the computation time by
$175$, $146$, $90$, $64$, $60$ times when $k=6,8,10,12,$ and $14$
correspondingly. Furthermore, although both methods’ computation time
increases with the number of candidate cutoffs, the ratio of the two methods
is gradually levelled off.

Note that the permutation approach carried out in the simulation study only
consists of $1000$ permutations, providing the estimate of the FWER with the
precision of $0.001$ (i.e. three digits from the decimal point). The precision
of these estimates is often not enough in many applications. If we need to
increase the precision of the FWER by one more digit from the decimal point,
the number of required permutations increases to $10$ folds, leading to at
least a $10$ fold increase in the computing time of the permutation approach.
In contrast, BOSS estimates FWER theoretically and its precision does not
interfere with the computing time.

###  3.4 Sensitivity Analysis on Sample Size

As a sensitivity analysis, we explore if the sample size has an effect on the
performance of the two methods. The previous simulation studies are based on
the same sample size with the real data employed to generate blueprint models.
So in this analysis we generate data with different sample sizes. In
particular, we conduct random sampling with replacement to these simulated
datasets separately to generate one with half size and the other with double
size from each original dataset. After that we apply the analysis scheme to
these generated data and compare their results with the ones from the original
data to make conclusions.

Figure 4:  Boxplots of type I error rate and power from sensitivity analysis.
Based on survival outcome under $k=10$, boxplots from the half, original, and
doubled sized data are presented in the left, middle and right penals
seperately. And the type I error rates and two types of power from the top
$50$ and the bottom $50$ biomarkers are presented in the top, middle and
bottom panels separately. Boxplots from BOSS are filled with white colors,
while boxplots from the permutation method are filled with grey colors for
comparison.

The sensitivity analysis results confirm that under different sample sizes,
BOSS and permutation perform very similarly in statistical power and Type I
error. For example, Figure 4 shows the distributions of estimated power and
type I error rates from further investigations under a fixed $k=10$ and
survival outcomes with various-sized data. Furthermore, these results indicate
that increasing the sample size can largely improve the statistical power for
both BOSS and permutation. However, changes in sample sizes do not affect the
type I error rates.

In summary, simulation studies demonstrate that the novel BOSS method works as
good as the gold standard permutation method regarding the ability in
detecting positive signals and controlling false positive discoveries, which
are not affected by sample size. However, BOSS uses much less time than the
permutation method, satisifying the critieria of a good load time for a
webtool and is more useful when we need to repeat analysis many times (e.g.
analyze many genes in a genomic study).

##  4 Real Data Application

Lung cancer is the leading cause of cancer death worldwide (Siegel, ). Lung
adenocarcinoma (LUAD) is a type of lung cancer that develops in the epithelial
cells of the lungs. Drug and immune resistance has been extensively reported
in LUAD Spella ; Iglesias . Revealing the genetic mechanisms underlying LUAD
can facilitate understanding LUAD etiology and hence offers new entry points
for treatment and/or disease management. To demonstrate the novel method
provides new insights into the cancer prognosis, we employed the BOSS to
analyze the curated LUAD genomic and clinical data published in cSurv , which
contains the expressions of $17431$ genes collected from $500$ LUAD cases and
their corresponding clinical information including gender, overall survival
time, and censoring status. This dataset can be downloaded from the website
provided in the Data availability statement.

The results from BOSS provide the optimal cutoff and the FWER (i.e. the
adjusted p-value defined in Formula 6) for each gene. Note that the FWER makes
a correction in BOSS to avoid p-value inflation due to multiple testing of
different cutoff values for the same gene. In this analysis, we conducted
analysis for $17431$ genes, which require another multiple testing adjustment.
We further corrected these $17431$ BOSS FWERs to control the False Discovery
Rate (FDR) of testing many genes at level $0.05$. After the FDR adjustment, we
found expressions of $48$ genes having optimal cutoffs and showing significant
(FDR-adjusted $p<0.05$) correlation with overall survival time in LUAD
(Supplementary Table S1). We performed a literature review of these $48$ genes
and found that $20$ of them were specifically reported to be associated with
LUAD (Supplementary Table S1). For example, Xu  demonstrated that
overexpression of ANLN is associated with metastasis of LUAD complications.
Liu et al. Liu  showed that CXCL induced spinal metastasis of LUAD. Xu et al.
Wu  identified FAM117A as a novel prognostic marker for poor outcomes in lung
cancer patients. The supports from independent evidence from the literature
demonstrate that BOSS can find useful prognostic genes and generate biological
insights. There are $28$ genes discovered by our method but have no literature
support. We consider them as potential novel genes and suggest researchers
further investigate and validate them.

Figure 5: Enrichment network of overrepresented (pval $\textless$ 0.0001, padj
$<$ 0.075) GSs in genes correlated with overall survival in lung
adenocarcinoma. Node denotes GS; node size reflects the number of genes in
each GS. Edge reflects significant gene overlap as defined by a Jaccard
Coefficient larger than or equal to 0.25. Detailed statistics are provided in
Supplementary Table S2. ECM, extracellular matrix; GS, gene set; padj,
adjusted P-value; pval, P-value.

To understand the common messages in these $48$ genes, we performed functional
enrichment analysis with eVITTA (eVITTA, ) (Supplementary Table S2). We found
that these $48$ genes are overrepresented in cancer hallmark pathways (Han, ),
including cell-extracellular matrix (ECM) interaction, xenobiotic transport,
and hypoxia response (Figure 5). Notable genes include ITGB1, FLNC, SLC2A1,
and EGLN1. Among them, ITGB1, $\beta$1 integrin, is the predominantly
expressed integrin protein in normal and tumor cells playing a pivotal role in
ECM-related signaling; its expression is found to be prognostic in non-small-
cell lung cancer (Weiqi, ). FLNC, or filamin-C, as an important cytoskeleton
component, plays a key role in cell migration and signal transduction (Ai, ).
While its role in lung cancer is not extensively reported, it’s been shown to
be prognostic in other cancers such as prostate (Ai, ). SLC2A1 or glucose
transporter 1, and EGLN1 or hypoxia-inducible factor prolyl hydroxylase 2, are
two key mediators of hypoxia responses (Dengler, ). Hypoxia is an important
feature lung cancer microenvironment related to cancer progression,
metastasis, and metabolism (Iwona, ). In hypoxia responses, SLC2A1 helps
transition from oxidative towards glycolytic metabolism in cancer cells while
EGLN1 modifies HIFs Post-translationally (Dengler, ). This network highlights
potential converging mechanisms that might influence malignancy
synergistically, advancing LUAD and promoting resistance to treatment.

Together, our method was robust in detecting potentially prognostic biomarkers
by examining all expressed genes with reasonable computational resources and
time. The biomarkers identified with BOSS in this work have referential
importance in exploring synergy and target-matched therapies in LUAD.

##  5 Conclusion and Discussion

Research on selecting the optimal cutoffs is conducted in various fields, such
as machine learning, biology, and medicine. In machine learning, Shih ; Lausen
used the optimal cutoff selection in building decision and regression trees.
In biology, optimal cutoffs are utilized to distinguish species abundance
based on ecological factors (Miiller, ). In clinical trials, the optimal
cutoffs are used in different settings, such as assessing the association
between the characteristics of tumor cells and the prognosis of patients
(Galon, ), and obtaining the threshold for residual leukemic cells in acute
myeloid leukemia patients (Buccisano, ). However, these proposed methods
cannot be not directly applied to our case, i.e. featuring a special
overlapping data structure generated from several candidate cutoffs.
Therefore, we developed our novel statistical test method, BOSS, which has
three main advantages: fast calculating speed, high accuracy, and no
requirement for a large number of candidate cutoffs. First, BOSS has decent
speed, leading to great practical significance. weng  shows that every
0.4-second increase in Google page load time reduces at least 8 million daily
searches. And every 1-second rise in Amazon’s Web load time reduces at least
$1.6$ billion company’s annual revenue. The speed of BOSS satisfies the good
load time criteria. Therefore, it can greatly improve the user experience if
it replaces the current permutation method used in the web system for medical
analysis (cSurv, ). Second, BOSS is accurate. And its accuracy is consistent
and does not interfere with the speed. In contrast, the accuracy of the
permutation test is proportional to the number of permutations and inversely
proportional to the speed. High accuracy requirement in permutation leads to
severly prolonged running time. Third, BOSS does not require a large number of
candidate cutoffs. BOSS is designed based on controlling family wise error
rate, not requiring a large number of hypothesis testings. Furthermore, we
carry out theoretical derivation to calculate the adjusted p values from the
exact false positive probabilities based on specific data overlapping
structures, which improve the accuracy and remove the requirement of large
sample asymptotic properties.

In practice, the number of candidate cutoffs is related to the power in two
contradictory directions. On one side, investigating more cutoffs enables
finding an optimal value closer to the truth, which increases the power. On
the other hand, investigating more cutoffs requires more tests causing heavier
penalties in the multiple testing adjustment, which decreases the power. Our
simulation studies provide some evidence and insights into that paradox. For
quantitative outcomes in Figure 1, most powers of genes with strong signals
are above $80\%$, and the best median power is achieved at the maximum number
of cutoffs $k=14$, whereas most power of genes with weak signals are between
$30\%$ and $40\%$ and the best median power is achieved at $k=10$. For
survival outcomes, most powers are $60\%$ or lower, and the best median powers
are not achieved at the max value of $k$. This suggests exploring more
candidate cutoffs to find a more precise optimal cutoff that can better
segment patients, when the signal is strong enough to achieve a good power of
$80\%$; and suggests reducing the number of cutoffs to lower the penalty of
multiple testing adjustments, when the signal is not strong to achieve a good
power. In summary, when selecting the number of candidate cutoffs, we should
balance between two factors to achieve the best results: (1) find a more
precise optimal cutoff and (2) reduce multiple testing penalties.

In this article, we discuss the method to find the optimal cutoffs with one
biomarker, which is the most common case in practice, i.e., only using one to
segment patients into subpopulations. However, sometimes researchers need to
consider more than one biomarker. For example, researchers who are interested
in two biomarkers often explore differences between samples with double
positive or double negative biomarkers (i.e., subgroups where biomarkers are
both high or low). Patients with two biomarkers’ values in opposite directions
(one protective and one risk) are excluded from personalized medicine in this
situation due to the difficulty of defining their group membership. We can
convert two continuous biomarkers into a binary variable $X_{i}$ and apply
BOSS to compare double positives versus double negatives. The detailed process
of transforming the two-biomarker problem into a one-biomarker problem is
discussed in the Supplementary document with a numerical example.

In summary, BOSS is a powerful tool for selecting optimal cutoffs in various
research fields, especially in the scenario requiring fast computing speed to
develop a web application or needing to repeat analysis many times (e.g. on a
vast number of biomarkers in genomic studies). The derivation of BOSS inspires
researchers to develop multiple test adjustments based on specific data
overlapping structures.

## Supplementary document

The supplementary document includes 4 parts: (1) the derivation of the
covariance formula of Z-statistics in linear regression models, (2)scenarios
of finding optimal cutoffs for multiple variables, (3) the table of 48 genes
with significant optimal cutoffs, and (4) overrepresented genes in cancer
hallmark pathways.

## Data availability statement

The cancer patient survival data were obtained from the cSurvival program
(cSurv, ). Raw survival data can be downloaded from
https://tau.cmmt.ubc.ca/cSurvival/project_data/TCGA-LUAD/df_survival_o.csv.
Raw gene data can be downloaded from
https://tau.cmmt.ubc.ca/cSurvival/project_data/TCGA-LUAD/df_gene.csv.

## References

  * (1) Messoud Ashina, Gisela M Terwindt, et al.  Migraine: disease characterisation, biomarkers, and precision medicine.  The Lancet, 397(10283):1496–1504, 2021.
  * (2) Ashley J. Vargas and Curtis C. Harris.  Biomarker development in the precision medicine era: lung cancer as a case study.  Nature Reviews Cancer, 16(8):525–537, 2016.
  * (3) Yoav Benjamini and Yosef Hochberg.  Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing.  Journal of the Royal Statistical Society: Series B (Methodological), 57(1):289–300, 1995.
  * (4) Rupert G. Miller.  Simultaneous Statistical Inference.  Springer New York, NY, 1981.
  * (5) K J Worsley.  An improved Bonferroni inequality and applications.  Biometrika, 69(2):297–302, 1982.
  * (6) Rupert G. Miller and David O. Siegmund.  Maximally selected chi square statistics.  Biometrics, 38:1011, 1982.
  * (7) Anne-Laure Boulesteix and Carolin Strobl.  Maximally selected Chi-squared statistics and non-monotonic associations: An exact approach based on two cutpoints.  Computational Statistics & Data Analysis, 51(12):6295–6306, 2007\.
  * (8) Berthold Lausen and Martin Schumacher.  Maximally Selected Rank Statistics.  Biometrics, 48(1):73, 1992.
  * (9) Torsten Hothorn and Achim Zeileis.  Generalized Maximally Selected Statistics.  Biometrics, 64(4):1263–1269, 2008.
  * (10) Susan Galloway Hilsenbeck and Gary M. Clark.  Practical p-value adjustment for optimally selected cutpoints.  Statistics in Medicine, 15(1):103–112, 1996.
  * (11) Susan Galloway Hilsenbeck, Gary M Clark, and William L McGuire.  Why do so many prognostic factors fail to pan out?  Breast Cancer Research and Treatment, 22(3):197–206, 1992.
  * (12) Xuanjin Cheng, Yongxing Liu, Jiahe Wang, Yujie Chen, Andrew Gordon Robertson, Xuekui Zhang, Steven J M Jones, and Stefan Taubert.  cSurvival: a web resource for biomarker interactions in cancer outcomes and in cell lines.  Briefings in Bioinformatics, 23(3):bbac090, 2022.
  * (13) Michael Wiegand.  Site speed is (still) impacting your conversion rate, 2022.  https://www.portent.com/blog/analytics/research-site-speed-hurting-everyones-revenue.htm.
  * (14) David Harrington.  Group Sequential Methods with Applications to Clinical Trials.  Chapman and Hall/CRC, NY, 2001.
  * (15) Sam Brilleman.  simsurv: Simulate survival data.  Website, 2021.  https://cran.r-project.org/web/packages/simsurv/index.html.
  * (16) Rebecca L. Siegel, Kimberly D. Miller, and Ahmedin Jemal.  Cancer statistics, 2019.  CA: A Cancer Journal for Clinicians, 69(1):7–34, 2019.
  * (17) Magda Spella and Georgios T. Stathopoulos.  Immune Resistance in Lung Adenocarcinoma.  Cancers, 13(3):384, 2021.
  * (18) Venus Sosa Iglesias, Lorena Giuranno, Ludwig J. Dubois, Jan Theys, and Marc Vooijs.  Drug Resistance in Non-Small Cell Lung Cancer: A Potential for NOTCH Targeting?  Frontiers in Oncology, 8:267, 2018.
  * (19) Jie Xu, Hong Zheng, et al.  Overexpression of ANLN in lung adenocarcinoma is associated with metastasis.  Thoracic Cancer, 10(8):1702–1709, 2019.
  * (20) Wangmi Liu, Xiankuan Xie, and Jiayan Wu.  Mechanism of lung adenocarcinoma spine metastasis induced by CXCL17.  Cellular Oncology, 43(2):311–320, 2020.
  * (21) Chao Wu, Jiajin Zhang, et al.  FAM117A Is a New Prognostic Marker of Lung Adenocarcinoma and Predicts Sensitivity to PD0332991.  Evidence-Based Complementary and Alternative Medicine, 2022:3945446, 2022.
  * (22) Xuanjin Cheng, Junran Yan, Yongxing Liu, Jiahe Wang, and Stefan Taubert.  eVITTA: a web-based visualization and inference toolbox for transcriptome analysis.  Nucleic Acids Research, 49(W1):W207–W215, 05 2021.
  * (23) Douglas Hanahan and Robert A. Weinberg.  Hallmarks of Cancer: The Next Generation.  Cell, 144(5):646–674, 2011.
  * (24) Weiqi Zheng, Caihui Jiang, and Ruifeng Li.  Integrin and gene network analysis reveals that itga5 and itgb1 are prognostic in non-small-cell lung cancer.  OncoTargets and Therapy, 9:2317–2327, 2016.  PMID: 27143927.
  * (25) Jianzhong Ai, Tao Jin, et al.  Vinculin and filamin-C are two potential prognostic biomarkers and therapeutic targets for prostate cancer cell migration.  Oncotarget, 8(47):82430–82436, 2017.
  * (26) Veronica L. Dengler, Matthew D. Galbraith, and Joaquín M. Espinosa.  Transcriptional regulation by hypoxia inducible factors.  Critical Reviews in Biochemistry and Molecular Biology, 49(1):1–15, 2014.
  * (27) Iwona Ziółkowska-Suchanek.  Mimicking Tumor Hypoxia in Non-Small Cell Lung Cancer Employing Three-Dimensional In Vitro Models.  Cells, 10(1):141, 2021.
  * (28) Y.-S Shih.  A note on split selection bias in classification trees.  Computational Statistics & Data Analysis, 45(3):457–466, 2004\.
  * (29) Berthold Lausen, Torsten Hothorn, Frank Bretz, and Martin Schumacher.  Assessment of Optimal Selected Prognostic Factors.  Biometrical Journal, 46(3):364–374, 2004.
  * (30) Jörg Müller and Torsten Hothorn.  Maximally selected two-sample statistics as a new tool for the identification and assessment of habitat factors with an application to breeding-bird communities in oak forests.  European Journal of Forest Research, 123(3):219–228, 2004.
  * (31) Jérôme Galon, Anne Costes, et al.  Type, density, and location of immune cells within human colorectal tumors predict clinical outcome.  Science, 313(5795):1960–1964, 2006.
  * (32) F Buccisano, L Maurillo, V Gattei, et al.  The kinetics of reduction of minimal residual disease impacts on duration of response and survival of patients with acute myeloid leukemia.  Leukemia, 20(10):1783–1789, 2006.
  * (33) Zhuang Weng.  The Optimization and Performance Testing of Web Application Based on HTTP/2protocol.  15(1), 2017.


#  Enso
A general-purpose virtual machine ††thanks: This work is licensed under the
Creative Commons Attribution-ShareAlike 4.0 International License.

Bruno França

Version 4 – August 1, 2019  Trinkler Software
<EMAIL_ADDRESS>
###### Abstract

In this paper we introduce Enso, a virtual machine designed to be used as
general-purpose state transition function in blockchains. This design allows
the blockchain application logic to be coded into the state, instead of into
the state transition function, making it much more flexible and easier to
modify. A byproduct is reducing the frequency of forks, concerted or not.

##  I Introduction

Despite the diversity and complexity in the blockchain world today, all
blockchains are just distributed virtual machines.

Independently of the function for which it is used, a blockchain is just a
virtual machine running some code distributed across several computers.

The components of such a distributed virtual machine can be visualized as a
stack:

State
---
STF
Consensus
Networking

Starting at the bottom, we have the networking. This basically includes all
the tools that the computer needs to connect to the other nodes and to
receive, transmit and relay data across the network. It is a standard
component for any peer-to-peer network.

Then we have the consensus. This allows all nodes to reach an agreement on
what blocks are accepted. There is a huge variety of consensus protocols, but
they all have the same objective: making all honest nodes reach an agreement
on some input data.

Lastly, we have the state and the state transition function (STF). State
represents, of course, the current state of the virtual machine. The state,
for most users, is the only component that they interact with. The state
transition function is the function that takes as input the previous state and
a block, and produces a new state as output. So, it is what defines the rules
for modifying the state. Contrary to the networking and consensus, these two
components no longer interact with, or have any notion of, other nodes.
Together they just form a simple virtual machine.

The way in which each of these components is modified, or upgraded, varies
from blockchain to blockchain. The following image illustrates how legacy
chains, like Bitcoin [1] or Ethereum [2], are modified:

State | $\leftarrow$ | Extrinsics
---|---|---
STF | $\leftarrow$ | Forks
Consensus | $\leftarrow$ | Forks
Networking | $\leftarrow$ | Adoption

The state is always modified with extrinsics111We use the same definition for
transaction as the one taken by Parity. There are extrinsics, which are any
input to the state transition function, and both transactions and inherents
are mutually exclusive types of extrinsics. Transactions are extrinsics that
are propagated through the network and signed. Inherents are neither
propagated nor signed. An example of an inherent is a timestamp.. This is how
most people will interact with the blockchain and it is very straightforward.

Next we have the state transition function and the consensus. Both can be
modified with forks, either soft forks or hard forks. The problem with forking
is that forks are difficult to organize and can cause the chain to split into
two networks (ex: Bitcoin and Bitcoin Cash, or Ethereum and Ethereum Classic).

Finally we have the network. Changes in the network normally are not
contentious since they can happen with gradual user adoption. For example,
some nodes in the Bitcoin network can start accepting transactions via Tor and
relaying them to miners. There is no need for any type of consensus on whether
transactions should be accepted via Tor. Any node that wishes to do so can
make that change and make that service available to the network right away.
So, networking updates happen with organic user adoption.

Recently, blockchains appeared which claim to be upgradable by governance. We
are only interested in the technical aspect of upgrading a blockchain and not
in the social or political aspect (e.g. if voting should be proportional to
stake, etc). The most famous examples in this category are probably Tezos [3]
and Parity’s Substrate (Substrate is not actually a blockchain but rather a
blockchain development kit). With this type of chain we have the following
stack:

State | $\leftarrow$ | Extrinsics
---|---|---
STF | $\leftarrow$ | Swaps
Consensus | $\leftarrow$ | Swaps
Networking | $\leftarrow$ | Adoption

Note that the state and networking are still modified in the same fashion as
legacy chains, the difference is that the state transition function (STF) and
the consensus use module swaps instead of forks.

In these blockchains, there is a shell that contains only the networking and
some other auxiliary code. Then there are two slots (one for the STF, another
for the consensus) that accept arbitrary code blocks, called modules.

The idea is that a blockchain can be modified by simply swapping one module
for another. If the blockchain also has governance functionality coded into
it, it is possible for the nodes to reach an agreement on whether a module
should be swapped or not.

Governance and module swaps are different components. Governance only allows a
network to reach an agreement on whether the blockchain is to be updated or
not. The actual method of doing that update is the module swapping.
Objectively, the main difference between forks and swaps is that swaps are
concerted, while forks are not. Otherwise the process is similar. The miners
(or validator, or authorities) all agree to, at a specific time, stop using
one module and start using the next one.

The breakthrough on the part of Tezos was realizing that legacy blockchains
have their governance off-chain and that instead the governance should be
placed on-chain into the STF.

Another way of analyzing this stack is seeing what functionality is coded into
each component. As we have said, a blockchain is a distributed virtual
machine, but of course, each blockchain is running a specific application
(also commonly called the runtime). For example, Bitcoin is a distributed
virtual machine running a ledger application. We can visualize where each part
of the application is located in the stack:

App data | $\equiv$ | State
---|---|---
App logic | $\equiv$ | STF
|  | Consensus
|  | Networking

The application logic and data are separated, with the logic being hardcoded
into the STF and the data being stored into the state.

This separation of code and data is how early computers worked. Programs were
hardcoded into the hardware using plugboards and then data would be fed into
the computer using some storage medium (for example, punch cards). This was
time-consuming, requiring the plugboard to be rewired whenever a new program
was to be run in the computer. Then along came the stored-program computer,
which treats code and data equally. The paradigm change was thinking of the
hardware as a general-purpose shell without any specific functionality
predefined into it. Now programs could be fed into the computer just like
regular data.

Enso applies the stored-program computer concept to blockchains. It adopts a
very simple general-purpose virtual machine as the state transition function
so that the application logic can be put into the state. This allows for both
the application logic and data to be modified with extrinsics, while the
consensus and the virtual machine itself are still modified with forks or
module swaps.

App logic and data | $\equiv$ | State | $\leftarrow$ | Extrinsics
---|---|---|---|---
Virtual machine | $\equiv$ | STF | $\leftarrow$ | Forks/Swaps
|  | Consensus | $\leftarrow$ | Forks/Swaps
|  | Networking | $\leftarrow$ | Adoption

By treating the logic and the data equally, blockchains can be built with a
greater degree of flexibility. And extrinsics allow for simpler and more
granular modifications to the application logic. Changes to the consensus or
to the virtual machine still need to be done using the module swaps, but we
expect that consensus and virtual machine updates will be fairly rare, at
least when compared with changes to the application logic.

##  II Actor model

The actor model is a mathematical model of computation, especially useful when
concurrency is desired, that was first published in 1973 and is inspired in
Physics. It treats all computation as the result of actors exchanging messages
between themselves.

An actor is an object with an address, code and storage. They communicate with
each other by sending messages. For this, each of them has a mailbox, which is
basically just a queue for the incoming messages. Upon receiving a message, an
actor will run it’s code and can alter it’s own storage.

An actor is thus capable of the following actions:

  * •

Pulling messages from its mailbox queue.

  * •

Deciding what to do with a message, including ignoring it.

  * •

Modifying its internal state.

  * •

Sending messages to other actors at addresses that it knows about.

  * •

Creating new actors.

There is no assumed sequence to the above actions and they can be performed in
parallel. Messages are asynchronous and thus may arrive in any order.

##  III Specification

In this section we will delve into the specification of the Enso virtual
machine. Enso is heavily modeled after the actor model, with the main
difference being that Enso does not support concurrency, although concurrency
could in theory be added.

Just like in the actor model, in Enso everything is either an actor or a
message.

An actor is an entity composed of the following components:

  * •

ID: A unique identifier of the actor.

  * •

Code: A block of code containing functions that can be called by other actors.

  * •

Storage: A data structure containing arbitrary information and that can be
read and modified by the code.

Using Enso, the state of the blockchain will be simply the set of all actors.
So, all of the application logic and data will be coded into the actors.

Note that actors in Enso do not have a mailbox. That is because, since no
concurrency is needed, there is only a single global message queue that orders
all messages sent by all actors.

Messages, put simply, are asynchronous function calls from one actor to
another. By asynchronous we mean that when an actor sends a message, the actor
will not be blocked while it waits for the message to be processed. Instead it
will keep working and accepting other messages. This asynchronicity is
necessary because we want the virtual machine to keep running even if a
function call fails.

A message takes the form (ID_to, function_call, parameters), where:

  * •

ID_to: The ID of the actor that will receive the message.

  * •

function_call: The name of the function that will be called.

  * •

parameters: A list of parameters that will be passed to the function.

Now that we know that the state is composed of actors, and that actors
communicate between themselves with messages, we can state which actions are
available to an actor in Enso:

  * •

Receive messages from the global message queue.

  * •

Processing messages by calling the corresponding internal functions.

  * •

Modifying its own ID, code or storage.

  * •

Send messages to other actors at addresses that it knows about.

  * •

Read storage of other actors at addresses that it knows about.

  * •

Create new actors.

Contrary to the actor model, actors in Enso are able to read each others
storage. The blockchain state is public and accessible to anyone, so it only
makes sense to allow actors to also read the entire state.

The only component of Enso that remains to be described is the global message
queue, which is just a simple first in first out queue for messages.

When an message is processed, it is taken out of the queue and the function
function_call of the actor ID_to is called with the input parameters
parameters. If, for some reason, the function cannot be called, for example if
the actor doesn’t exist or there is no function with that name, the message is
simply ignored.

##  IV Implementation

The objective of the last section was to describe the components of the Enso
virtual machine and its operation in an abstract form. In this section we want
to give some details on how Enso can be implemented in practice.

First, even though both the application data and logic will coexist in the
state, in most cases there will be a clear separation between both. By this we
mean that, when developing a new blockchain with Enso, it is useful to divide
the actors into two categories:

  * •

Kernel actors: All actors created at the genesis block and during updates to
the blockchain. These are actors that define the rules for how the runtime
works and are unique actors, meaning that there is only one instance of each
actor type.

  * •

User actors: All actors created by the users. There is a predetermined list of
actor templates, from which users can create new user actors.

The idea is that the entire application logic will be coded into the kernel
actors, for example we may have one actor to handle the incoming extrinsics,
another that implements some governance method, etc. While the user actors
will be regular actors like accounts or smart contracts, effectively
representing the application data.

Second, an extrinsic is just a message that is sent by an user, or by some
other external source, instead of by an actor. Otherwise, an extrinsic should
be able to be processed by an actor just like a message.

Third, a block is just an ordered list of extrinsics. When the virtual machine
receives a new block, it just needs to process those extrinsics in order and
in the end it will have reached a new state.

So, to process a new block, a node follows these rules:

  1. 1.

Receive an ordered list of extrinsics.

  2. 2.

Add all the extrinsics, in order, to the global message queue.

  3. 3.

Take the top message and send it to the corresponding actor.

  4. 4.

The receiving actor then may create more messages, which will also be added to
the queue.

  5. 5.

Keep repeating steps 3 and 4 until there are no more messages left

  6. 6.

When there are no more messages to be processed, the state transition has
finished.

##  V Conclusion

This paper analyzes the organization of a blockchain’s code by thinking of it
as a distributed virtual machine. We then argue that a distributed virtual
machine needs to act analogously to a stored-program computer, which implies
that the application logic needs to be coded into the state and that the state
transition function needs to be a general-purpose virtual machine. Using this
method it is possible to have blockchains that are simpler and faster to
develop and operate, and that require far fewer forks.

## References

  * [1] S. Nakamoto _et al._ , “Bitcoin: A peer-to-peer electronic cash system,” 2008\. [Online]. Available: https://bitcoin.org/bitcoin.pdf
  * [2] V. Buterin _et al._ , “A next-generation smart contract and decentralized application platform,” _white paper_ , 2014. [Online]. Available: https://cryptorating.eu/whitepapers/Ethereum/Ethereum_white_paper.pdf
  * [3] L. Goodman, “Tezos: a self-amending crypto-ledger,” _URL: https://www. tezos. com/static/papers/white_paper. pdf_ , 2014. [Online]. Available: https://everipedia-storage.s3.amazonaws.com/NewlinkFiles/16739988/8d4d1f1b-8/white_paper.pdf