# Statistical monitoring of European cross-border physical electricity flows
using novel temporal edge network processes

Anna Malinovskaya
Leibniz University Hannover, Germany
Rebecca Killick
Lancaster University, United Kingdom
Kathryn Leeming
Leicestershire, United Kingdom
Philipp Otto
University of Glasgow, United Kingdom

###### Abstract

Conventional modelling of networks evolving in time focuses on capturing
variations in the network structure. However, the network might be static from
the origin or experience only deterministic, regulated changes in its
structure, providing either a physical infrastructure or a specified
connection arrangement for some other processes. Thus, to detect change in its
exploitation, we need to focus on the processes happening on the network. In
this work, we present the concept of monitoring random Temporal Edge Network
(TEN) processes that take place on the edges of a graph having a fixed
structure. Our framework is based on the Generalized Network Autoregressive
statistical models with time-dependent exogenous variables (GNARX models) and
Cumulative Sum (CUSUM) control charts. To demonstrate its effective detection
of various types of change, we conduct a simulation study and monitor the
real-world data of cross-border physical electricity flows in Europe.

Keywords: CUSUM Control Charts; GNARX Model; Network Modelling; Network
Monitoring.

##  1 Introduction

Network analysis is an important statistical inference tool in various
disciplines including economics, social sciences, and chemistry (De Masi and
Gallegati, 2007; Khrabrov and Cybenko, 2010; McDermott et al., 2021); It
provides insight into the network structure and beyond. Usually, the
underlying network structure is considered to be random, meaning that both
nodes and edges may appear and disappear over time. However, processes also
originate in network structures with limited or no temporal dynamics. In both
cases, to understand whether a change has occurred we need to inspect the data
over time. This procedure is also known as network monitoring. Hence, to
distinguish the case when the network is considered to be a random observation
from the case when the network structure is deterministic but the process on
it is random, we utilise the terminology random network monitoring and fixed
network monitoring introduced in Stevens et al. (2021b).

Usually, the modelling and monitoring of networks are focused on detecting
changes and anomalies reflected in the geometric properties of a graph,
falling under the type of random network monitoring. Illustrative studies
include the monitoring of e-mail communication and the daily flights within a
country for detecting unusual communication patterns (cf. Perry, 2020;
Malinovskaya and Otto, 2021). However, for some networks, the development of
connections or inclusion of new nodes is either no longer possible or would
not be of considerable monitoring interest.

In general, the integration of nodal or edge attributes for improving the
modelling of the underlying network formation mechanism and in turn the
network monitoring is not a new concept (Azarnoush et al., 2016; Shaghaghi and
Saghaei, 2020). However, in existing research, the attributes are considered
to regulate the presence or absence of an edge. In other words, the contextual
information available either on vertices or edges is viewed as an extra
dimension to the graph helping to explain the likelihood of the observed links
(Miller et al., 2013). In contrast, we consider the process happening on edges
as the primary or only source of information about the network state.

It is beneficial to distinguish general random processes we regard on the
fixed network structures from specific areas of analysis of processes on
networks. The first well-established perspective is the analysis of spatial
point patterns on networks (cf. Baddeley et al., 2021), where the accurate
location of the object on the physical network is of concern. The second area
concerns the analysis of network flows (cf. Ahuja et al., 1993; Kolaczyk and
Csárdi, 2020), where the physical constraints of a network structure and the
flow itself play a vital role. In this work, however, we could think of the
analysed process as being a flow or traffic with no associated constraints
unless they are imposed explicitly.

Starting with the description of TEN processes in Section 2, we explain our
monitoring framework in Section 3. Afterwards, we perform a simulation study
by testing the proposed methodology under different anomalous scenarios in
Section 4. To illustrate the importance and the idea of how to monitor the
real-world TEN processes, we monitor cross-border physical electricity flows
in Europe in Section 5. In the last section, we discuss further research
direction and summarise the perspectives and limitations of the proposed
approach.

##  2 TEN Process

Consider a network $G=(V,E)$, where the elements of $V$ represent vertices (or
nodes) and $E$ – edges (or links). Further, we assume a fixed structure of $G$
over time described by an adjacency matrix
$\bm{Y}\coloneqq(Y_{ij})_{i,j=1,\dots,|V|}$, with $|V|$ being the number of
nodes. When two vertices $i,j\in V$ are connected by an edge $e\in E$, they
are called adjacent. Usually, the entries of $\bm{Y}$ are binary, i.e.
$Y_{ij}=1$ if $(i,j)\in E,\,i\neq j$, and $0$ otherwise. The graph can be
directed or undirected, in the latter case is $\bm{Y}$ symmetric. To each
existing edge $e\in E$ we relate time series $\\{x_{e,t}\\}$, where
$t=1,\dots,T$, being the attributed process of $G$. The complete
representation $\bm{X}=(X_{e,t})_{e=1,\ldots,|E|,t=1,\ldots,T}$, where $|E|$
denotes the number of edges, together with the graph $G$ forms a dynamic
Temporal Edge Network (TEN) process illustrated in Figure 1 (top) for two time
stamps. The following section proposes one approach to modelling and
monitoring TEN processes.

Figure 1: Visual illustration of a TEN process (top) for two time points with
an adapted representation (below) described in Section 3.1.2.

##  3 Network Monitoring Framework

Usually, statistical network monitoring is subdivided into two parts: 1)
Network modelling and 2) Process monitoring (Stevens et al., 2021a). Instead
of formal network modelling, it is possible to collect network features and
monitor those. For example, Flossdorf and Jentsch (2021) comprehensively
discuss which network characteristics are suitable for which types of change.
In our work, prior to the monitoring procedure in Section 3.2, we first define
a model that can be suitable for TEN processes.

###  3.1 Modelling of TEN Processes

The majority of research on network models considers the network (or graph)
structure to be a random variable. Available models include Markov process
models (cf. Snijders, 2005), (temporal) exponential random graph models (cf.
Hanneke et al., 2010) and latent process-based models, e.g. (dynamic)
stochastic block models (cf. Matias and Miele, 2017). In contrast, TEN
processes assume a static structure of a network over time. One could adapt
the aforementioned random network models to this context but they involve
computationally intensive estimation that relies on numeric approximation.
Hence, we consider alternative model forms designed for fixed network
structures.

The modelling of multivariate time series with an underlying network
dependency structure is gaining popularity. Knight et al. (2020) present the
Generalized Network Autoregressive (GNAR) framework, where a network is
associated with multivariate time series and modelled as one. The GNAR models
time series that occur at nodes of the network, e.g. growth rates of gross
domestic product. The GNAR model was extended to incorporate time-dependent
exogenous variables (GNARX) in Nason and Wei (2022). After the introduction of
the original GNARX model in Section 3.1.1, we present its extension from nodal
time series to time series on network edges in Section 3.1.2.

####  3.1.1 GNARX Model

In this section, consider $\\{x_{i,t}\\}$ to be a time series related to each
node $i$. The GNARX model $(p,\bm{s},\bm{q})$ with $H$ exogenous regressors
$\\{z_{h,i,t:h=1,\dots,H,\,i\in V,\,t=1,\dots,T}\\}$ and the autoregressive
order $p$, where
$(p,\bm{s},\bm{q})\in\mathbb{N}\times\mathbb{N}^{p}_{0}\times\mathbb{N}^{H}_{0}$
holding for all vertices $i\in V$, is specified as

$x_{i,t}=\sum_{l=1}^{p}\bigg{(}\alpha_{i,l}x_{i,t-l}+\sum_{r=1}^{s_{l}}\beta_{l,r}\sum_{j\in\mathscr{N}^{(r)}(i)}\omega_{i,j}x_{j,t-l}\bigg{)}+\sum_{h=1}^{H}\sum_{q=0}^{q_{h}}\gamma_{h,q}z_{h,i,t-q}+\epsilon_{i,t}.$
(1)

The order $p$ also determines the maximum order of neighbour time lags, i.e.
$\bm{s}=(s_{1},\dots,s_{p})$ with $s_{l}$ being the maximum stage of neighbour
dependence for time lag $l$. For example, $s_{1}=2$ means that nodes depend on
their first and second-stage neighbours in $G$ in the first time lag.
Similarly, the maximum time lag of the $h$-th exogenous regressor
$\\{z_{h,i,t}\\}$ is defined as $q_{h}$ and collectively these are
$\bm{q}=(q_{1},\dots,q_{H})$. In case $q_{h}=0$, the current value of the
exogenous variable is considered at time point $t$.

The noise is denoted by $\epsilon_{i,t}$ and is assumed to be independent and
identically distributed at each vertex $i$ with mean zero and variance
$\sigma^{2}_{i}$. The parameters $\alpha_{i,l},\ \beta_{l,r},\
\gamma_{h,q}\in\mathbb{R}$ define autoregressive influence, neighbouring
influence and external influence from regressors, respectively. It is possible
to estimate a global-$\alpha$ model, where $\alpha_{i,l}=\alpha_{l}$, assuming
the same autoregressive process for all nodes.

The set $\mathscr{N}^{(r)}(i)$ denotes the $r$-th stage neighbourhood set of
node $i\in G$. For instance, the stage-1 neighbours of a node $i\in V$ are the
adjacent nodes $j\in V$, connected by an edge. Further, the stage-2
neighbourhood set of a node $i\in V$ is the stage-1 neighbours of the adjacent
nodes $j\in V$ as can be seen in Figure 2. Moreover, there are weights
$\omega\in[0,1]$ associated with every pair of nodes that, in our case, depend
on the size of the neighbour set as explained in Knight et al. (2020).

Figure 2: Stage-1 (blue) and stage-2 (green) neighbourhood of the orange node.

By fitting the GNARX model, we obtain the estimates of the parameters
$\alpha_{i,l},\ \beta_{l,r},\ \gamma_{h,q}$ that can be used for predicting
the $\hat{x}_{i,t+1}$ value from $\\{x_{i,t}\\}$. The one-step-ahead forecast
errors are then determined as

$u_{i,t+1}=x_{i,t+1}-\hat{x}_{i,t+1}$ (2)

which can be utilised in our monitoring framework. However, the GNARX model is
defined for time series related to nodes. Thus, we need to extend the
representation of TEN processes to enable a correct application of the GNARX
model for performing the monitoring.

####  3.1.2 Extension to TEN Processes

There are different ways of thinking about why we require an alternative
representation. First, our main focus lies in discovering changes happening
between and within the streams, i.e. a process captured on the edges. Thus, we
actually indirectly consider it to be a network not of single nodes anymore
but of pairs of nodes. Second, the influence from node to node and the
influence from one stream to its neighbouring stream portray two distinct
standpoints. In the case of the dependency between two traffics, we have more
than two nodes involved and other underlying mechanisms that determine the
exchange and its strength. Therefore, we need a different adjacency matrix to
describe this novel dependency structure.

For distinguishing between $G$ and $\bm{Y}$, and the novel representation of
the TEN process as $G^{\prime}$, we introduce notation for the new adjacency
matrix as $\bm{Y}^{\prime}$. Now, as displayed in Figure 1 (below), the edges
$e\in E$ from the original graph $G$ become vertices $\iota\in V^{\prime}$
with the number of nodes being $|V^{\prime}|=|E|$ so that the time series
$\\{x_{e,t}\\}$ placed on edges $e\in E$ are now placed on the nodes, becoming
$\\{x_{\iota,t}\\}$. Consequently, we would also have new edges $\xi\in
E^{\prime}$. Considering the connectivity structure captured by the adjacency
matrix, without any expert knowledge, the natural choice for creating
$\bm{Y}^{\prime}$ would be to use a model for generating random graphs. A
potential candidate is Erdős-Rényi model, where all graphs on a fixed vertex
set with a fixed number of edges are equally likely. Another possibility is to
consider sampling from a stochastic block model which produces graphs with a
community structure. For example, links within a community may be more
probable than between communities. Both of those options are tested in Section
4.

There is, however, an additional approach to choosing a new connectivity
structure by referring to other types of graphs from graph theory. In this
case, we use a deterministic approach to construct $\bm{Y}^{\prime}$. One
relevant representation technique is the line (also known as edge-to-vertex
dual) graph $L(G)$ of $G$. Here, $L(G)$ is constructed on $E$, where $e\in E$
are adjacent as nodes $\iota\in V^{\prime}$ if and only if they share a common
node as edges in $G$ (Diestel, 2017). In Figure 1 (below) this connectivity
structure is adapted to the original TEN process shown above. Such formation
can suit well real-world applications as it offers a more logical structure in
case of limited knowledge about any existing communities or insufficiency of a
connectivity structure generated from a random graph model.

After redefining the representation of TEN processes, we can proceed with
their modelling by applying the GNARX model and their monitoring using a
residual-based control chart introduced in the next section.

###  3.2 Monitoring of TEN Processes

Recall that for network monitoring we need a model and a monitoring procedure.
For monitoring, it is vital to determine the framework in which we are aiming
to perform it. One possibility is to monitor the estimates of the model
parameters although this is computationally intensive (cf. Malinovskaya and
Otto, 2021). A more efficient alternative is to perform monitoring based on
residuals. For instance, Miller et al. (2013) compute graph residuals as the
difference between the observed graph and its expected value. In modelling the
TEN with the GNARX model (see Equation 2), it is uncomplicated to calculate
the residuals as the deviations from the current observation and the GNARX
prediction. As soon as a change occurs, we would expect the residuals to
increase indicating a lack of fit. According to Alexopoulos et al. (2004), who
offers a detailed introduction to Statistical Process Monitoring (SPM) and
especially to the forecast-based monitoring methods, when the model or
prediction is accurate, the prediction errors are uncorrelated. Thus, we can
apply traditional SPM techniques such as control charts to these forecast
errors. A technique that is particularly suitable for testing whether the
current residual-based statistic is anomalous is a control chart.

Belonging to the methods of SPM (cf. Montgomery, 2012), control charts are
both straightforward to implement and interpret. Grundy (2021) proposes a
framework for monitoring forecast models, showing that changes in complex data
are reflected in the forecast errors. They utilise the Cumulative Sum (CUSUM)
control chart, and we extend their approach to monitoring residuals of the
GNARX model for TENs as shown subsequently. A comprehensive overview of
residual control charts and their comparison to other types of control charts
is provided in Knoth and Schmid (2004) and Jensen et al. (2006).

Applying control charts in an online manner, i.e. performing a sequential
change point detection, at each time point $t$ we test whether the process
remains in control. Hence, the specific hypothesis as well as the definition
of a change point are clarified in Section 3.2.1 before coming to the
description of the implemented control chart in Section 3.2.2.

####  3.2.1 Hypothesis Formulation

As we base our monitoring part on the residuals obtained by finding the
difference between the actual process and the predictor at each flow (Equation
2), it is vital to obtain an accurate model before monitoring starts. Our
objective is then to test

$H_{0,t}:\,\text{The observed TEN process coincides with the fitted GNARX
model}$

against the alternative

$\hskip 36.98866ptH_{1,t}:\,\text{The observed TEN process does not coincide
with}\text{the fitted GNARX model}$

for each $t$. Now the question arises of which variable to use for the test
statistic that could provide us with the information about when $H_{0}$ is
violated. The change point $\tau$ is defined as

$\bm{x}_{t}\quad\sim\quad\left\\{\begin{array}[]{cc}F(\bm{\mu},\,\Sigma)&\text{
if }t<\tau\\\ F_{\tau}(\bm{\mu}_{\tau},\,\Sigma_{\tau})&\text{ if
}t\geq\tau\\\ \end{array}\right.\,,$

where $\bm{x}_{t}=(x_{1,t},\dots,x_{|V^{\prime}|,t})^{\prime}$, $\bm{\mu}$ and
$\Sigma$ define the mean vector and variance-covariance matrix of the network
flow distribution $F$ in our case, respectively. However, the change in the
mean or/and in the variance of the raw data would lead to the respective
changes in the forecast errors (Grundy, 2021). Hence, by constructing a test
statistic based on $\bm{u}_{t}=(u_{1,t},\dots,u_{|V^{\prime}|,t})^{\prime}$,
we can accordingly determine the time point $\tau$ when the change has
occurred. Following, we introduce the change point detection framework
specified by Grundy (2021) and discuss the respective test statistic for
applying residual-based CUSUM control charts.

####  3.2.2 Residual-Based CUSUM Control Chart

Before describing the technicalities of the considered control chart, it is
worth beginning with an explanation of its key elements and framework for
application in practice. In general, a control chart, being a graphical tool
for detecting unusual variations of the process, consists of three components.
These are a central line ($CL$) for plotting the average of the process, an
upper control line for the upper control limit and a lower control line for
the lower control limit. When the process considerably deviates from its
expected state, the control statistic starts crossing the limits, indicating a
substantial change beyond random variation.

Usually, the are two phases involved in the implementation of control charts.
Phase I corresponds to the exploration and calibration period, where we
estimate the target process and calibrate the control chart, e.g. computing
control limits. It is based on the assumption that the process during this
phase is in control, i.e. stable, predictable and repeatable (Vining, 2009).
In Phase II the actual application of the control chart begins which could be
viewed as a sequential implementation of the hypothesis test introduced in
Section 3.2.1. The more precise the estimation of the parameters in Phase I,
the more reliable the performance of the control chart in Phase II. The
desired operation of control charts during Phase II is a quick detection of
the change point in a process when it experiences an out-of-control state,
i.e. unusual variation of a process is present. At the same time, we strive
for a long-running scheme without false alarms occurring, meaning the case
when the process actually remains in control but the scheme signals a change.
This translates to a balance between small control limits for fast detection
and large control limits for small false alarm rates.

In real-world applications, it is usually unknown what properties, e.g. mean
or variance, of a process may change. Thus, an important criterion for
selecting a suitable control chart is its ability to track many types of
change simultaneously. Another point considers the decision about what exactly
to monitor – the data itself or some process related to it. As discussed by
Grundy (2021), the monitoring of the data directly might deteriorate in case
of temporal dependency, complex trends or seasonality effects within its
structure. However, monitoring forecast residuals of the process omits these
issues and is capable of reflecting changes in either the process mean or the
process variance.

As under the real settings, it is usually unknown whether the change occurs in
a mean and/or variance of the process, we utilise a more general type of
Page’s CUSUM detector, which is based on the (centred) squared data and is
able to detect a combination of the mean and/or variance change in the
original data. For monitoring forecast errors $u_{\iota,t}$ obtained for each
flow $\iota$, Grundy (2021) adapts Page’s CUSUM test statistic (cf. Page 1954)
to the centered squared forecast errors as follows

$Q_{\iota}(m,k)=\sum_{t=m+1}^{m+k}(u_{\iota,t}-\hat{b})^{2}-\frac{k}{m}\sum_{t=1}^{m}(u_{\iota,t}-\hat{b}),$
(3)

where $m$ corresponds to the length of Phase I, $k$ is the current time point
in Phase II and $\hat{b}$ is the mean estimate of the forecast errors computed
from Phase I. From that, we compute the control chart statistic as

$D_{\iota}(m,k)=\max_{0\leq a\leq k}|Q(m,k)-Q(m,a)|,$ (4)

and the corresponding upper control limit is given by

$UCL=\hat{\sigma}_{\iota}\zeta_{\alpha}g(m,k,\nu),$ (5)

with $\zeta_{\alpha}$ being the critical value of the distribution with the
significance level $\alpha$ and $\hat{\sigma}_{\iota}$ the estimate of the
standard deviation of the centred squared forecast errors belonging to the
flow $\iota$. The part $g(m,k,\nu)$ defines the weight function with the
tuning parameter $\nu$, both are described in detail in Horváth et al. (2004)
and Grundy (2021). In this work, we consider $\nu=0$. The current time point
$k=\tau$, i.e. the control chart detects a change point if $D_{\iota}(m,k)\leq
UCL$. It is worth noting that no lower control limit is defined as the control
statistic obtains non-negative values only and, therefore, has its natural
lower control limit which is zero.

####  3.2.3 Monitoring Scope and Performance Function

The considered CUSUM chart belongs to the univariate control charts, meaning
that only one process variable is monitored by the scheme. In our case, it
means that we simultaneously implement $n=|V^{\prime}|$ control charts to
monitor the TEN process completely. Consequently, it is also possible to
perform local monitoring and control the behaviour of only one essential flow
from a complete TEN process.

In this work, however, we are interested in proposing a monitoring procedure
suitable for tracking the complete TEN process. Thus, as a performance metric,
we introduce a cumulative intensity function of a change that presents a
cumulative sum of flows that triggered a change at a time point $t$

$I_{\bm{X}}(T)=\sum_{t=1}^{T}\frac{\sum_{\iota=1}^{n}{\mathbbm{1}_{[UCL,\infty)}[D_{\iota}(m,k)]}}{n}.$
(6)

If in one flow $\iota$ a change was detected, the monitoring of this flow
stops, resulting in $\max(I)=1$. Thinking about how a signal is produced by
computing $I_{\bm{X}}(T)$, we need to define a suitable threshold $W$ based on
expert knowledge or system requirements.

In the subsequent section, we illustrate a simulation study that handles
different cases of changes as well as compares monitoring results with respect
to the selected model for constructing the adjacency matrix $\bm{Y}^{\prime}$.

##  4 Simulation Study

The reason to conduct a simulation study is twofold: First, we aim to quantify
the detection speed of changes in the effects captured by the model parameters
$\alpha_{i,l}$, $\beta_{l,r}$ or $\gamma_{h,q}$. Second, as mentioned in
Section 2, there are different possibilities to construct the adjacency matrix
$\bm{Y}^{\prime}$ as introduced in Section 3.1.2. Here, we examine two
distinct random graph models that cover both the availability and absence of
expert knowledge about suitable sampling models. If a TEN process is well
understood, the application of a Stochastic Block Model (SBM), where the flows
should be first subdivided into clusters, is a good choice. Alternatively, it
is possible to run a clustering algorithm. If the cluster structure is not
suitable or no specific knowledge about the TEN process is available, the
structure of $G^{\prime}$ can be sampled by applying Erdős-Rényi model. An
example of both structure types is presented in Figure 3.

Figure 3: Representation of two different connectivity types: A graph
structure generated by an Erdős-Rényi model (left) vs. by a SBM (right).

In the next part, we explain the design of the conducted simulation study,
presenting the results as well as comparing them with regard to the generated
graph structure.

###  4.1 Experimental Setting

To perform a simulation study, we have to decide on the design of a TEN
process, the initial parameters and the test cases. To keep the study simple
for reproducing it but comprehensive for designing different anomaly types as
well as testing two distinct $\bm{Y}^{\prime}$ structures, we create a TEN
process with 10 flows and medium connectivity. That means, for sampling
$\bm{Y}^{\prime}$ from the Erdős-Rényi model, we set the number of nodes to be
10 and the number of connections to be 30, considering that everything else is
random. For sampling from an SBM, we separate the flows into two clusters
$c_{1}$ and $c_{2}$ of equal size. Additionally, we need to define the
probabilities $P$ that the flows are connected within a cluster or between the
clusters.

Regarding the settings for sampling from a GNARX model, we design a process
with a global autoregressive effect of order $p=1$. Also, the stage-1
neighbourhood is considered together with two exogenous regressors and
corresponding time lags $\bm{q}=(0,0)$. To be precise, we generate separate
data for estimating the GNARX model (600 observations) before designing Phase
I with 200 observations and Phase II with 100 observations which is subdivided
into 50 in-control and 50 out-of-control samples. The detailed description of
each simulation step is presented in Procedure 1. In total, we conduct 500
iterations for each test case, presenting the outcomes subsequently.

Procedure 1

  1. 1.

Setting: Choose the lengths of the burn-in period (300), model fitting period
(600), Phase I (200) and Phase II (100) subdivided into in-control (50) and
out-of-control (50) parts. Set the data-generating parameters of the GNARX
model: $\alpha_{\iota,l}=\alpha=0.2$, $\beta_{l,r}=\beta=0.3$ with $p=1$ and
$\bm{s}=(1)$, $\gamma_{h,q}=\gamma_{1}=2$ and $\gamma_{2}=3$ with
$\bm{q}=(0,0)$. Select the change: $\beta_{\tau}=\beta+0.3$.

  2. 2.

TEN Generation:

    1. (a)

Generate the adjacency matrix
$\bm{Y}^{\prime}\coloneqq(Y^{\prime}_{\iota\upsilon})_{\iota,\upsilon=1,\dots,|V^{\prime}|}$
either by sampling from the Erdős-Rényi model with $|V^{\prime}|=10$ and
$|E^{\prime}|=30$ or the SBM with $P_{c_{1}c_{2}}=P_{c_{2}c_{1}}=0.2$ and
$P_{c_{1}c_{1}}=P_{c_{2}c_{2}}=0.8$.

    2. (b)

For each $t=1,\dots,1200$ generate the GNARX process:

      1. i.

Sample error values $\epsilon_{\iota,t}$, where
$\epsilon_{\iota,t}\sim\mathcal{N}(0,1)$.

      2. ii.

Sample covariate values $z_{1,\iota,t}\sim\mathcal{N}(0,1)$ and
$z_{2,\iota,t}\sim\mathcal{N}(0,1)$.

      3. iii.

For $\iota=1$ to $\iota=|V^{\prime}|$

        1. A.

If $t\leq 1150$, generate the data with $\alpha,\beta,\gamma_{1}$ and
$\gamma_{2}$ parameters:

$x_{\iota,t}=\alpha
x_{\iota,t-1}+\beta\sum_{\upsilon\in\mathscr{N}^{(1)}(\iota)}\omega_{\iota,\upsilon}x_{\upsilon,t-1}+\gamma_{1}z_{1,\iota,t-1}+\gamma_{2}z_{2,\iota,t-1}+\epsilon_{\iota,t}.$

        2. B.

If $t\geq 1151$, generate the data with
$\alpha,\beta_{\tau},\gamma_{1},\gamma_{2}$ parameters having an effect either
on all flows or only on flows from $c_{1}$, adjusting the equation in A.

  3. 3.

TEN Modelling:

    1. (a)

Fit the GNARX model with $p=r=1$ using the first 600 generated observations
after discarding the burn-in period $t=1,\dots,300$.

    2. (b)

Save the estimated parameters
$\tilde{\alpha},\tilde{\beta},\tilde{\gamma}_{1}$ and $\tilde{\gamma}_{2}$.

  4. 4.

TEN Monitoring:

    1. (a)

Phase I: Estimate the target process by using
$\tilde{\alpha},\tilde{\beta},\tilde{\gamma}_{1}$ and $\tilde{\gamma}_{2}$
parameters to compute forecast residuals $u_{\iota,t}$ and calibrate the CUSUM
control charts described in Section 3.2.2.

    2. (b)

Phase II: Run monitoring using the calibrated CUSUM control charts and compute
the cumulative intensity function provided in Equation 6. Save the time stamp
$t$ when a change is detected.

Procedure 1 Simulation steps for one iteration for the test case $\beta+0.3$

###  4.2 Results and Comparison of Connectivity Structures

As explained in Section 3.2, we aim to detect changes in the whole TEN
process, i.e. the anomalous network states. Therefore, we compute the
cumulative change intensity given in Equation 6. For summarising the
performance in the simulation study, for each time point in Phase II we
average the values of the cumulative change intensity function, obtaining a
mean cumulative change intensity. It is worth noting that we do not define a
threshold $W$ here, focusing on the general detection capability of the
approach as well as fluctuation in the performance. It is important to note
that as soon as the control chart signals for a particular flow, the
monitoring of that flow stops. It means, that no multiple change points for
the same flow are possible.

For each parameter $\tilde{\alpha},\tilde{\beta}$ and $\tilde{\gamma}_{1}$ we
have three test cases of changes whose values gradually increase. Comparing
the different types of sampling the adjacency matrix $\bm{Y}^{\prime}$, we
exchange the initially estimated parameters with anomalous either in all flows
or only in flows being part of $c_{1}$.

Change in $\alpha+0.2$ (Original $\alpha=0.2$)

Change in $\alpha+0.3$ (Original $\alpha=0.2$)

Change in $\alpha+0.5$ (Original $\alpha=0.2$)

Figure 4: Simulation study: Visualisation of the mean of cumulative change
intensity functions $I_{\bm{X}}(T)$ over 500 iterations with differences in
$\alpha$.

Change in $\beta+0.3$ (Original $\beta=0.3$)

Change in $\beta+0.4$ (Original $\beta=0.3$)

Change in $\beta+0.5$ (Original $\beta=0.3$)

Figure 5: Simulation study: Visualisation of the mean of cumulative change
intensity functions $I_{\bm{X}}(T)$ over 500 iterations with differences in
$\beta$.

Change in $\gamma_{1}+1.0$ (Original $\gamma_{1}=2.0$)

Change in $\gamma_{1}+1.5$ (Original $\gamma_{1}=2.0$)

Change in $\gamma_{1}+2.0$ (Original $\gamma_{1}=2.0$)

Figure 6: Simulation study: Visualisation of the mean of cumulative change
intensity functions $I_{\bm{X}}(T)$ over 500 iterations with differences in
$\gamma_{1}$.

In the following Figures 4, 5 and 6 we display the monitoring results in Phase
II (PII). The blue area defines the in-control and the red the out-of-control
part. First, in all plots, we notice the same pattern: The curve that shows
the simulation with $\bm{Y}^{\prime}$ sampled from the SBM stays below another
curve that corresponds to the setting with an adjacency matrix generated from
the Erdős-Rényi model. That corresponds to the way the change in the effects
was implemented: While in the case of the SBM structure, only one community
(half of the flows) was affected by anomalous parameter sizes, the simulation
that involves the Erdős-Rényi structure experiences the change in the complete
network, i.e. all flows were affected. Thus, the cumulative intensity change
reaches a value of only 0.5 which is the normalised size of a community with 5
flows that presents anomalous behaviour when using the SBM. Second, we clearly
recognise that the most difficult type of anomaly for the proposed approach is
the change in the neighbourhood effects captured by the parameter $\beta$. It
can be seen from the wider uncertainty span and longer run length to detect
the insertion of a change in $\beta$. Overall, the simulation study reveals
that the monitoring approach is highly effective, meaning that no signals
occur when no actual change has been introduced but is only efficient in
detecting quickly the time stamp after the change if we are aware of the
underlying structure of $G^{\prime}$ or aim to detect medium or large process
innovations.

In case one specific flow reflects the state of the whole system or is
particularly relevant for the flawless functionality of a network, we could
focus solely on its monitoring. For that, we would apply the CUSUM control
chart and directly determine the change point. Visually, a possible outcome
could look as displayed in Figure 7.

Change in $\alpha+0.3$ (Original $\alpha=0.2$)

Figure 7: Simulation study: Example of monitoring one specific flow from a TEN
process.

In the following section, as a real-world illustration, we present the
monitoring of cross-border physical electricity flows across Europe.

##  5 Monitoring of Cross-Border Physical Electricity Flows

Cross-country analysis of electricity trade, especially in the context of
renewable energy, is a relevant field of study that i.a. contributes to the
evaluation of transmission infrastructure policies (cf. Abrell and Rausch,
2016). It is worth noting that the cross-border physical flow network
indicates the actual flow of electricity that corresponds to the laws of
physics and could differ from the scheduled commercial exchanges which reflect
the economic relations between the market parties and, therefore, be of
particular technical importance, allowing for the flawless electricity trade
across Europe. To the best knowledge of the authors, there are no studies
related to the change point detection in European Cross-Border Physical Flows
(CBPFs), so we illustrate the monitoring of these processes below.

###  5.1 Data Description

European Network of Transmission System Operators for Electricity (ENTSO-E)
provides data about the cross-border physical flow and scheduled commercial
electricity exchanges of more than 40 countries111Central collection and
publication of electricity generation, transportation and consumption data and
information for the pan-European market. ENTSO-E Transparency platform (2023).
https://transparency.entsoe.eu. It offers the data for each hour for years
starting from 2014. However, after careful investigation of the amount of
missing data, we decided to concentrate on the period from January 1, 2018, to
November 27, 2022, and use a weekly aggregation of the electricity values
given in megawatts.

The graph in Figure 8 (a) represents the aggregated quantity of electricity
which was transmitted in the year 2019 so that the wideness of edges reflects
the total amount of power which was passed between the two countries. The
higher this value is, the wider the edge between two nodes. However, what we
can also notice is the existence of parallel edges, i.e. one flow $f_{1}$
coming from France (FR) to Spain (ES) and another flow $f_{2}$ back. If we
decide to directly apply a new representation as introduced in Section 3.1.2,
we would obtain a network that would be considerably bigger than the original
graph. Hence, we need to aggregate both flows $f_{1}$ and $f_{2}$ and take one
vertex to represent a country pair. There are different possibilities for
doing it, we decide on three of them and discuss each of them subsequently.

(a) The edge thickness implicates the strength of the flow.

(b) The size of the nodes implicates the strength of the electricity exchange
and the colour its proportion of electricity generated with renewable energy
sources (light green corresponds to the higher proportion, blue corresponds to
the missing information).

Figure 8: Conventional representation of the cumulated cross-border physical
electricity flow in the year 2019 (a) and new representation (b) following the
description in Section 3.1.2.

###  5.2 Phase I Modelling

To avoid a considerable expansion in the new representation of the CBPFs (see
Section 3.1.2), we introduce three different aggregation strategies of both
flows $f_{1}$ and $f_{2}$, fitting and comparing three models, respectively.
The first statistic is the Box-Cox transformation of the sum of both flows
$\mathscr{M}_{1}=ln(f_{1}+f_{2}+1)$ with $\lambda_{1}=0$ and $\lambda_{2}=1$
(cf. Box and Cox, 1964) that reflects the overall strength of the exchange.
The second statistic measures the asymmetry in both flows as the difference of
the Box-Cox transformed variables $f_{1}$ and $f_{2}$, i.e.
$\mathscr{M}_{2}=\ln(f_{1}+1)-\ln(f_{2}+1)$. It implies whether it is more
common for one country to import or export electricity. The third statistic
captures proportionally differences in both flows, being
$\mathscr{M}_{3}=\frac{f_{1}-f_{2}}{f_{1}+f_{2}}$. No additional
transformation is required in this case as
$\mathscr{M}_{3}\in[-1,1]$\$\\{0\\}$. Figure 8 (b) illustrates a new
representation of CBPFs, where the new adjacency structure is constructed
using the rules of a line graph described in Section 3.1.2.

As no exceptional events are known in year 2018 and 2019 that could
considerably disturb the CBPF network, we decide to choose these two years as
Phase I. Testing different model settings, the final set-up involves a global
autoregressive parameter $\tilde{\alpha}$ with $p=1$, 1-stage neighbourhood
$\tilde{\beta}$ and two covariates that correspond to the Box-Cox transformed
total amount of electricity generated by renewable energy sources from both
countries with the effects $\tilde{\gamma_{1}}$ and $\tilde{\gamma_{2}}$,
having $\bm{q}=(0,0)$.

Table 1 presents the results of estimating the parameters for each choice of
the aggregation statistic $\mathscr{M}$ using 7828 observations
($76\times(52\times 2-p)$), where 76 defines the number of bilateral
exchanges. As we can observe, the GNARX model fits the data in Phase I well
and confirms the relevance of accounting for the network structure as well as
external effects when modelling TEN processes. Using the displayed
coefficients, we proceed with the monitoring of Phase II which consists of the
years 2020–2022.

| Model $\mathscr{M}_{1}$ | Model $\mathscr{M}_{2}$ | Model $\mathscr{M}_{3}$
---|---|---|---
$\tilde{\alpha}$ | $0.896^{***}$ | $0.904^{***}$ | $0.910^{***}$
| $(0.005)$ | $(0.005)$ | $(0.005)$
$\tilde{\beta}$ | $0.094^{***}$ | $0.019^{\cdot}$ | $0.020^{\cdot}$
| $(0.006)$ | $(0.011)$ | $(0.011)$
$\tilde{\gamma}_{1}$ | $0.003^{\cdot}$ | $-0.034^{***}$ | $-0.003^{***}$
| $(0.002)$ | $(0.006)$ | $(0.001)$
$\tilde{\gamma}_{2}$ | $0.005^{\cdot}$ | $0.032^{***}$ | $0.003^{***}$
| $(0.002)$ | $(0.006)$ | $(0.001)$
$\bar{\mathscr{M}}$ | 10.84 | 0.14 | 0.03
$\tilde{\sigma}_{\mathscr{M}}$ | 2.00 | 6.81 | 0.84
RMSE | 0.89 | 2.75 | 0.34
${}^{***}p<0.001$, ${}^{**}p<0.01$, ${}^{*}p<0.05$, ${}^{\cdot}p\leq 0.1$
Table 1: Phase I modelling.

###  5.3 Phase II Monitoring

Figures 9, 10 and 11 display the monitoring results during Phase II. Out of 76
country pairs, between 27 and 39 pairs have had a change point according to
the implemented CUSUM control charts based on either $\mathscr{M}_{1}$,
$\mathscr{M}_{2}$ or $\mathscr{M}_{3}$ statistics. We highlight two periods
during Phase II: The period related to severe situations due to the COVID-19
pandemic and the period starting in February 2022 where several crises are
expected or happened due to the Russian-Ukrainian war. Considering the
threshold $W$, we define it to be $W=0.2$, i.e. when in at least 16 bilateral
exchanges a change point is detected.

Comparing three different aggregation methods, we notice similarities in
detecting changes during the pandemic. However, starting from June 2021 the
behaviour of the control charts significantly differs. In the case of the
results with $\mathscr{M}_{1}$, the time between two defined periods remains
relatively stable compared to two other charts with $\mathscr{M}_{2}$ and
$\mathscr{M}_{3}$. Then, the first detection on March 6, 2022 (aggregating the
days February 28 – March 6) considers the beginning of the war between Russia
and Ukraine. Surprising could be the fact that another country pair
experiences a change in the same week, namely Belarus and Ukraine. Later,
another change detection occurs at the end of May in the pair Russia and
Estonia, reflecting the political agreement of the Baltic states to stop any
energy trade with Russia. The same date and the reason are also given in the
control chart with $\mathscr{M}_{3}$ in the pair Russia and
Lithuania222https://enmin.lrv.lt/en/news/no-more-russian-oil-gas-and-
electricity-imports-in-lithuania-from-sunday. In terms of the signal in the
pair Ukraine and Moldova (control chart with $\mathscr{M}_{2}$), it could be
related to the launch of the planned commercial exchanges announced at the end
of June 2022 and further being increased in the beginning of autumn
333https://www.entsoe.eu/news/2022/09/04/transmission-system-operators-of-
continental-europe-decide-to-further-increase-trade-capacity-with-the-ukraine-
moldova-power-system/.

Figure 9: Results of monitoring $\mathscr{M}_{1}$, obtaining in total 27
change points. Figure 10: Results of monitoring $\mathscr{M}_{2}$, obtaining
in total 33 change points. Figure 11: Results of monitoring $\mathscr{M}_{3}$,
obtaining in total 39 change points.

Overall, we can notice that without expert knowledge of the CBPF network, it
is challenging to reflect a particular reason for the change point, however,
we also can see a reliable performance of the proposed framework in detecting
major events such as the pandemic and the energy crisis related to the begin
of the Russian-Ukrainian war. The caveat to be aware of is the interpretation
of the chosen aggregation statistic
$\mathscr{M}_{1},\mathscr{M}_{2},\mathscr{M}_{3}$ as in some of the cases the
signalled anomalies considerably differ.

##  6 Conclusion

The network with a given structure but a random process on its edges that we
define in this work as a Temporal Edge Network (TEN) can be of particular
interest for guaranteeing the safety of the infrastructure but also for
foreseeing possible accidents. In this manuscript, we present the monitoring
framework to detect anomalies in TENs by combining the GNARX model and the
CUSUM control chart based on residuals.

There might still exist the question of why monitoring networks with a fixed
structure needs special treatment. For example, why not try to “randomise” the
connections? The explanation is straightforward: It could be possible to
select a flow threshold for creating dynamics in the graph by deleting or
adding the links, however, the anomalies which would be detected in this case
are similar to those which are the focus of the random network monitoring. In
our case, we are rather interested in anomalies occurring in the edge process,
e.g. changes in the flow’s strength or some other temporal deviations. Hence,
we introduce a new way of representing TEN processes where edges become the
main focus of both modelling and monitoring parts.

The proposed change point detection framework is applied only to the TENs
observed at discrete times. It is an open research field on how and when to
extend the monitoring to continuous times. Equally important is to research
when different adjacency matrices for different time points should be
introduced and whether it will benefit the overall monitoring procedure. Also,
for monitoring a TEN process at once, we have introduced a cumulative change
intensity function. However, an appealing way to perform such monitoring would
be a suitable multivariate control chart or another statistical monitoring
tool that allows for medium-sized multivariate processes.

As an empirical illustration, we have monitored a network of bilateral
electricity flows across Europe. Beyond that, our framework could potentially
be applied for monitoring traffic flows on roads or rail systems to optimise
transportation infrastructure as well as within the computer communication
domain, aiding in cybersecurity. In an environmental context, TENs may be
useful in evaluating water transfers between areas or companies, contributing
to sustainable water resource management.

This work is devoted to networks with static structures that fit well with the
considered process on edges being cross-border physical electricity flow.
However, thinking of potential processes with a changing underlying structure
over time, it is important to determine whether the monitoring framework would
experience substantial changes in terms of re-estimating the parameters in
Phase I or challenges in performing multivariate monitoring when some of nodes
or edges disappear with time. Moreover, the effectiveness of the proposed
monitoring procedure strongly relies on the assumption of the GNARX model
being suitable to the considered data. Thus, in case the data cannot be well
represented by the GNARX model, i.e. the count data are given, an alternative
modelling approach or an extension to the currently existing model is
required.

## References

  * (1)
  * Abrell and Rausch (2016) Abrell, J. and Rausch, S. (2016).  Cross-country electricity trade, renewable energy and european transmission infrastructure policy, Journal of Environmental Economics and Management 79: 87–113.
  * Ahuja et al. (1993) Ahuja, R. K., Magnanti, T. L. and Orlin, J. B. (1993).  Network Flows: Theory, Algorithms, and Applications, Prentice Hall.
  * Alexopoulos et al. (2004) Alexopoulos, C., Goldsman, D., Tsui, K.-L. and Jiang, W. (2004).  SPC monitoring and variance estimation, Frontiers in Statistical Quality Control 7: 194–209.
  * Azarnoush et al. (2016) Azarnoush, B., Paynabar, K., Bekki, J. and Runger, G. (2016).  Monitoring temporal homogeneity in attributed network streams, Journal of Quality Technology 48(1): 28–43.
  * Baddeley et al. (2021) Baddeley, A., Nair, G., Rakshit, S., McSwiggan, G. and Davies, T. M. (2021).  Analysing point patterns on networks—a review, Spatial Statistics 42: 100435.
  * Box and Cox (1964) Box, G. E. and Cox, D. R. (1964).  An analysis of transformations, Journal of the Royal Statistical Society Series B: Statistical Methodology 26(2): 211–243.
  * De Masi and Gallegati (2007) De Masi, G. and Gallegati, M. (2007).  Debt-credit economic networks of banks and firms: the Italian case, Econophysics of Markets and Business Networks: Proceedings of the Econophys-Kolkata III, Springer, pp. 159–171.
  * Diestel (2017) Diestel, R. (2017).  Graph Theory, Springer.
  * Flossdorf and Jentsch (2021) Flossdorf, J. and Jentsch, C. (2021).  Change detection in dynamic networks using network characteristics, IEEE Transactions on Signal and Information Processing over Networks 7: 451–464.
  * Grundy (2021) Grundy, T. D. (2021).  On Aspects of Changepoint Analysis Motivated by Industrial Applications, Lancaster University (United Kingdom).
  * Hanneke et al. (2010) Hanneke, S., Fu, W., Xing, E. P. et al. (2010).  Discrete temporal models of social networks, Electronic Journal of Statistics 4: 585–605.
  * Horváth et al. (2004) Horváth, L., Hušková, M., Kokoszka, P. and Steinebach, J. (2004).  Monitoring changes in linear models, Journal of statistical Planning and Inference 126(1): 225–251.
  * Jensen et al. (2006) Jensen, W. A., Jones-Farmer, L. A., Champ, C. W. and Woodall, W. H. (2006).  Effects of parameter estimation on control chart properties: a literature review, Journal of Quality technology 38(4): 349–364.
  * Jeske et al. (2018) Jeske, D. R., Stevens, N. T., Tartakovsky, A. G. and Wilson, J. D. (2018).  Statistical methods for network surveillance, Applied Stochastic Models in Business and Industry 34(4): 425–445.
  * Khrabrov and Cybenko (2010) Khrabrov, A. and Cybenko, G. (2010).  Discovering influence in communication networks using dynamic graph analysis, 2010 IEEE Second International Conference on Social Computing, IEEE, pp. 288–294.
  * Knight et al. (2020) Knight, M., Leeming, K., Nason, G. and Nunes, M. (2020).  Generalized network autoregressive processes and the GNAR package, Journal of Statistical Software 96(5): 1–36.
  * Knoth and Schmid (2004) Knoth, S. and Schmid, W. (2004).  Control charts for time series: A review, Frontiers in statistical quality control 7: 210–236.
  * Kolaczyk and Csárdi (2020) Kolaczyk, E. D. and Csárdi, G. (2020).  Analysis of Network Flow Data, Springer International Publishing, Cham, pp. 169–186.
  * Malinovskaya and Otto (2021) Malinovskaya, A. and Otto, P. (2021).  Online network monitoring, Statistical Methods & Applications 30(5): 1337–1364.
  * Matias and Miele (2017) Matias, C. and Miele, V. (2017).  Statistical clustering of temporal networks through a dynamic stochastic block model, Journal of the Royal Statistical Society Series B: Statistical Methodology 79(4): 1119–1141.
  * McDermott et al. (2021) McDermott, M. J., Dwaraknath, S. S. and Persson, K. A. (2021).  A graph-based network for predicting chemical reaction pathways in solid-state materials synthesis, Nature communications 12(1): 3097.
  * Miller et al. (2013) Miller, B. A., Arcolano, N. and Bliss, N. T. (2013).  Efficient anomaly detection in dynamic, attributed graphs: Emerging phenomena and big data, 2013 IEEE International Conference on Intelligence and Security Informatics, IEEE, pp. 179–184.
  * Montgomery (2012) Montgomery, D. C. (2012).  Statistical Quality Control, Wiley Global Education.
  * Nason and Wei (2022) Nason, G. P. and Wei, J. L. (2022).  Quantifying the economic response to COVID-19 mitigations and death rates via forecasting purchasing managers’ indices using generalised network autoregressive models with exogenous variables, Journal of the Royal Statistical Society Series A: Statistics in Society 185(4): 1778–1792.
  * Page (1954) Page, E. S. (1954).  Continuous inspection schemes, Biometrika 41(1/2): 100–115.
  * Perry (2020) Perry, M. B. (2020).  An EWMA control chart for categorical processes with applications to social network monitoring, Journal of Quality Technology 52(2): 182–197.
  * Shaghaghi and Saghaei (2020) Shaghaghi, M. and Saghaei, A. (2020).  PCA likelihood ratio test approach for attributed social networks monitoring, Communications in Statistics-Theory and Methods 49(12): 2869–2886.
  * Snijders (2005) Snijders, T. A. (2005).  Models for longitudinal network data, Models and methods in social network analysis 1: 215–247.
  * Stevens et al. (2021a) Stevens, N. T., Wilson, J. D., Driscoll, A. R., McCulloh, I., Michailidis, G., Paris, C., Parker, P., Paynabar, K., Perry, M. B., Reisi-Gahrooei, M. and Sparks, R. (2021a).  Research in network monitoring: Connections with SPM and new directions, Quality Engineering 33(4): 736–748.
  * Stevens et al. (2021b) Stevens, N. T., Wilson, J. D., Driscoll, A. R., McCulloh, I., Michailidis, G., Paris, C., Paynabar, K., Perry, M. B., Reisi-Gahrooei, M., Sengupta, S. and Sparks, R. (2021b).  Foundations of network monitoring: Definitions and applications, Quality Engineering 33(4): 719–730.
  * Vining (2009) Vining, G. (2009).  Technical advice: Phase I and phase II control charts, Quality Engineering 21(4): 478–479.


2021

[1]Mamata Das

[1]Selvakumar K These authors contributed equally to this work.

These authors contributed equally to this work.

1]Department of Computer Applications, NIT Trichy, NH 83, Trichy, 620015,
Tamil Nadu, India

# Analyzing and Comparing Omicron Lineage Variants Protein-Protein Interaction
Network using Centrality Measure

<EMAIL_ADDRESS><EMAIL_ADDRESS>P.J.A. Alphonse
<EMAIL_ADDRESS>[

###### Abstract

The Worldwide spread of the Omicron lineage variants has now been confirmed.
It is crucial to understand the process of cellular life and to discover new
drugs need to identify the important proteins in a protein interaction network
(PPIN). PPINs are often represented by graphs in bioinformatics, which
describe cell processes. There are some proteins that have significant
influences on these tissues, and which play a crucial role in regulating them.
The discovery of new drugs is aided by the study of significant proteins.
These significant proteins can be found by reducing the graph and using graph
analysis. Studies examining protein interactions in the Omicron lineage
(B.1.1.529) and its variants (BA.5, BA.4, BA.3, BA.2, BA.1.1, BA.1) are not
yet available. Studying Omicron has been intended to find a significant
protein. $68$ nodes represent $68$ proteins and $52$ edges represent the
relationship among the protein in the network. A few centrality measures are
computed namely page rank centrality (PRC), degree centrality (DC),closeness
centrality (CC), and betweenness centrality (BC) together with node degree and
Local Clustering Co-efficient(LCC). We also discover $18$ network clusters
using Markov clustering. $8$ significant proteins (candidate gene of Omicron
lineage variants) were detected among the $68$ proteins, including AHSG,
KCNK1, KCNQ1, MAPT, NR1H4, PSMC2, PTPN11 and, UBE21 which scored the highest
among the Omicron proteins. It is found that in the variant of Omicron
protein-protein interaction networks, the MAPT protein’s impact is the most
significant.

###### keywords:

Omicron lineage variant, Centrality, PPIN, MCL

##  1 Introduction

An impartial panel of scientists known as TAG-VE (Technical Advisory Group on
SARS-CoV-2 Virus Evolution) regularly observes and examines the appraise of
the SARS-CoV-2 virus to determine if specific mutations or combinations of
mutations have an impact on the behaviour of the virus. The B.1.1.529 variants
of SARS-CoV-2 was the subject of an evaluation by the TAG-VE on November 26,
2021. South Africa disclosed the B.1.1.529 version on November 24, 2021 to
World Health Organization (WHO) 10665-62429 for the first time. The WHO has
classified B.1.1.529 as a VOC under the name Omicron considering the data that
a negative shift in COVID-19 epidemiology has occurred. Similar to other SARS-
CoV-2 variations, there are numerous lineages and sublineages in the Omicron
variation. Omicron presently has $3$ main lineages: BA.5, BA.4, and BA.2.
Omicron Pango Lineage currently has six different variations or sublineages
(BA.5, BA.4, BA.3, BA.2, BA.1.1, BA.1). Although these lineages are frequently
extremely similar to one another, there may be variations between lineages
that influence how the virus behaves. In our research, we have created seven
PPI networks of Omicron Pango Lineage including all the variants. The network
has been created on STRING, analyzed the network and find the most influential
proteins from the network. The networks that describe the interactions between
the parts of such complex systems are easier to analyse than it is to
investigate each component separately. The placement of some significant or
influential elements in most networks such as crucial proteins in PPI networks
is a well-known fact in the analysis of biological and social networks. These
locations, or vertices, have some unique structural characteristics. Such
facts are quantified using various centrality metrics. The vertices and edges
of a graph can be ranked from several perspectives based on centrality
measurements. To pinpoint "central" nodes in extensive networks, numerous
centrality measures (CM) have been developed. The user can choose whatever
metric best fits the study of a certain network because there are several
options available for ranking influential nodes. The effect of the network
architecture on how influential nodes are ranked by centrality metrics further
complicates the selection of an appropriate measure. In order to find the
centrality metric that is most successful at predicting influential proteins,
we looked at the centrality profiles of the nodes of Omicron PPINs. We looked
at how a broad range of widely used centrality measurements reflects various
topological network properties. This study demonstrates the state-of-the-art
in biological network centrality estimations. In order to identify the most
significant protein in the network, this research presents $4$ centrality
metrics (page rank centrality (PRC), degree centrality (DC), closeness
centrality (CC), and betweenness centrality (BC)) that are added with some
significant scores (node degree and Local Clustering Co-efficient(CCo), and
p-value) on Omicron variant’s PPI networks.

##  2 Related Work

Graph structures known as biological networks and social networks can be used
to describe a variety of complex systems, including biological and social
systems ghasemi2014centrality . For determining significant functional
characteristics of a network das2021markov das2022analytical , selecting an
appropriate set of centrality measurements is essential. landherr2010critical
the paper has been considered in relation to a critical analysis of centrality
measures in social networks. Three straightforward conditions for the
behaviour of centrality measures were used to analyse certain centrality
measures (BC, CC, DC, and eigenvector centrality). The author has been
analysis of PPI using Skyline Query on Parkinson’s disease
diansyah2019analysis . One of the disorders with the highest rate of global
growth, Parkinson’s disease, was shown to have 12 important proteins. The PPI
network features have been represented by attributes based on centrality
measures. The target genes for cancer illnesses were discovered by the author
using protein-protein interaction networks amala2019identification . Hubs and
centrality measurements were used to examine the possible genes. They extract
the genes with the highest scores in both mutation rates and graph centrality
in order to identify the target genes. The author compared 27 popular
centrality measurements using yeast PPINs ashtiani2018systematic . The
measurements classify and arrange the networks’ influential nodes. They have
also used hierarchical clustering and principal component analysis (PCA), and
they discovered that the topology of the network affects which metrics are the
most useful. The author has provided both historical and contemporary research
on social network centrality measures in das2018study survey paper. They
discussed created centrality measurements and mathematical definitions.
Additionally, they demonstrate various centrality measure uses in the fields
of education research grunspan2014understanding , biology
ghasemi2014centrality , traffic jayaweera2017centrality , transportation
wang2017novel , and security sparrow1991application carrington2011crime .
There are so many applications of centrality measure in different field
network brohee2006evaluation like psychological networks
bringmann2019centrality khojasteh2022comparing , brain networks joyce2010new ,
differential privacy models laeuchli2022analysis , etc.

##  3 Methods

This study used Omicron lineage variants data. The research has been completed
in different steps like, data collection, data cleaning, data validation,
creation of PPIN data, centrality measure and finally clustered the whole
network in different clusters. The clustering is done by the MCL (Markov
clustering algorithm) satuluri2009scalable . The objective of this research
work is to get the significant protein or prioritize the protein. For this we
have focuses on the centrality measure of the network. Figure 1 illustrates
the research workflow.

###  3.1 Data Collection

We have taken the real dataset of Omicron from Universal Protein
Resource/SwissProt (UniProt/ SwissProt) uniprot2021uniprot database which is
reviewed and found in the human body. In addition to storing experimental
results, computational features, and scientific conclusions, Swiss-Prot is a
highly annotated, non-redundant protein sequence database. Currently, the
UniProt Knowledgebase is comprised of UniProtKB/Swiss-Prot, which has been
reviewed. It provides accurate, consistent, and rich annotations for
functional information about proteins. Initially, we have taken a total of 228
proteins: B.1.1.529 (27), BA.5 (30), BA.4 (31), BA.3 (34), BA.2 (38), BA.1.1
(34), BA.1 (34) and analyzed individual Omicron lineage PPIN. The PPIN of
Omicron Lineage Variants are shown in Figure 7 to Figure 7. Then we sum up the
data and cleaned the data by removing duplicate data entries to create the
Omicron PPIN. The data validation and PPIN data creation in all the cases are
done by STRING szklarczyk2015string . There are several sources of information
within the STRING database, including computational prediction methods,
experimental data, and public text collections. A regular update keeps it up-
to-date and it is free to access. Additionally, it generates network images
using a spring model. In this model, nodes are considered masses, and edges
are considered springs. After cleaning the data we gate unique 68 proteins
which create the Omicron PPIN.

###  3.2 Centrality Measure

Here we will discuss very interasting aspect of network measure called
centrality. Centrality is basically widely used measure of how central a
particular node is with respect to the network.The network that results from
the PPI data is thought to be an undirected graph. Each node’s weight in the
graph is determined by the centrality approach. The BC, CC, DC, and PRC are a
few centrality techniques that can be applied to undirected graphs. Figure 7
to 7 and 15 depicts a protein network as an example of an undirected graph.
The variant BA.1 and BA.1.1 has the same PPIN only the difference in mutation.
The edges of the graph reflect the functional interaction or relationship that
takes place between proteins, whereas the nodes in the graph demonstrate the
proteins that affect Omicron’s activity.

####  3.2.1 Degree Centrality

The first basic centrality measure is the degree centrality (DC)
freeman1977set . We know that the degree is basically the number of edges
which are adjacent on a particular node. The DC is esentially is a degree of a
node but it is normalized.

The DC of a node $\mathit{v}$ is a degree of the node $\mathit{v}$ and divided
by the maximum degree of a node present in the graph. A node’s degree
centrality $C_{d}\mathit{(v)}$ in a network G(V, E) is denoted mathematically
as follows:

$C_{d}\mathit{(}v)=\frac{deg\mathit{(}v)}{\mbox{max
deg}_{u\in\mathit{v}}\mathit{(}u)}.$ (1)

It basically ranges between $0$ to $1$ and more the degree centrality mean
higher the likelyhood that the node has maximum degree. The
$C_{d}\mathit{(}v)$ can use to identify the more prominent or influential node
from a network.

####  3.2.2 Closensee Centrality

The Closensee Centrality ($CC$) sabidussi1966centrality indicates how close a
node from the rest of the network. A approach to identify nodes that can
efficiently spread information throughout a graph is through their CC. Average
distance between a node and all other nodes is measured by its proximity
centrality. The distances between nodes that have a high proximity score are
the shortest. A node’s Closeness centrality $C_{c}(\mathit{v})$ in a graph
G(V, E) is denoted mathematically as follows:

$C_{c}\mathit{(}v)=\frac{\lvert V\rvert-1}{\sum_{u\in
V-{\mathit{\\{v\\}}}}d\mathit{(}u,v)}.$ (2)

where, number of nodes is given by $\lvert V\rvert$ and the distance between
two nodes $\mathit{u}$ and $\mathit{v}$ is represented as $\mathit{d(u,v)}$.
Higher the value of $CC$ better would be the quality of the particular node.
The measure is useful in examining or restricting the spreds of disease in
epidemic modelling.

Start Data Collection BA.2 BA.3BA.4BA.5BA.1.1BA.1BA.1.259 Candidate Gene
(Omicron) Data Cleaning (68 unique protein) Merge all protein (242 protein)
PPIN Construction using STRING Validation of Protein in STRING Network
Analysis Node Centrality Closeness Centrality Betweenness Centrality Page Rank
Centrality $C_{S}\geq Th_{vl}$ Weak Interaction Strong Interaction Significant
Protein Endyesno Figure 1: Flow chart for numerical solution procedure

####  3.2.3 Betweenness Centrality

Betweenness (BC) freeman1977set is the measure to compute how central a node
is in between paths of the network or we can say to compute how many
paths(shortest) of the network passes through the node. A node’s Betweenness
centrality $C_{b}(\mathit{v})$ in a network G(V, E) is denoted mathematically
as follows:

$C_{b}\mathit{(}v)=\sum_{xy\in
V-\mathit{\\{v\\}}}\frac{\sigma_{xy}(\mathit{v})}{\sigma_{xy}}$ (3)

where the frequency of shortest paths in the network between nodes $x$ and $y$
is indicated by $\sigma_{xy}$ and $\sigma_{xy}(\mathit{v})$ denotes the same
passing through $\mathit{v}$. If $x=1$, then $\sigma_{xy}=1$. The BC is useful
in identifying the super spreaders in analyzing disease spreading in
epidemiology.

####  3.2.4 Page Rank

PageRank centrality ivan2011web is an adaptation of Eigen centrality that
ranks web content by using the value of linkages between sites. Any type of
network, including protein interaction networks, can be used with it.
Mathematically, the pagerank Centrality $C_{P}R(\mathit{v})$ in a network
$G(V,E)$ of a node $\mathit{v_{i}}$ is defined as:

$C_{PR}\mathit{(}v_{i})=\frac{1-d}{|V|}+d\sum_{\mathit{(}v_{t})\in
Inneighbor\mathit{(}v_{i})}\frac{C_{PR}(\mathit{v_{t}})}{outdeg(\mathit{v_{t}})}$
(4)

Where $d$ is constant and called damoing factor, usually the constant value is
considered as $0.85$.

###  3.3 Markov Clustering

At the Centre for Mathematics and Computer Science in the Netherlands, Stijn
van Dongen created the Markov Cluster Algorithm, MCL algorithm
vandongen2000cluster . It is an unsupervised cluster approach for networks
that is extremely quick and scalable and is based on the simulation of graph
flow. It is employed in bioinformatics and other fields. The distance matrix
derived from the STRING global scores in our study serves as the input to MCL.
Higher global scores for these interacting proteins increase the likelihood
that they will cluster together. The MCL das2021markov das2022analytical
operates primarily in two ways: Expansion-the operation corresponds to the
multiplication of standard matrices and simulates how a flow spreads and
becomes more homogeneous. The next is inflation which is described logically
as a diagonal scaling proceeded by a Hadamard power. Flow is compressed by
inflation by thickening only in areas where current density is high and
thinning only in areas where current density is low. There is no way to know
how many clusters there are. With the help of the inflation parameter, it is
implicitly managed. Higher inflation results in more clusters being obtained,
which is indirectly connected to the clustering’s precision. Here, the
inflation value has been set at $2$.

##  4 Results and Discussion

The global properties of Omicron base lineage variants are shown in 2. All the
seven network except $BA.1.1.259$ has an average node degree greater than $1$.
The $3$ base lineage ($BA.1.1.259,BA.1,BA.1.1$) has same density $0.0284$. The
highest density is $0.06719$ (BA.4) and the lowest density belongs to BA.3
($0.00416$). The average LCC is pretty good (highest $0.771$). The best
network is the BA.2 with the smallest p-value($0.00038$). The Table 2 shows
the global features of the Omicron PPIN. Node degree is $1.53$ on average and
the density is $0.0228$. The information in Table 4 contains the centrality
scores of $68$ proteins, which allow us to identify the protein’s relevance.
The network has a maximum degree of $7$ with an average Local Clustering
Coefficient (LCC) of $0.385$. The LCC range from $0$ to $1$, and they
represent the density of connections among neighbours. Nodes that have higher
values belong to densely connected clusters. The node is considered a part of
the clique if it has a value of $1$. The proteins GRB7, KCNK17, NDUFB5,
NDUFV1, RPSA, SNRPB, SNRPD1, and SNRPE in Table 4 are containing CCo value as
$1$ as they are part of the clique. Figure 15 is showing the PPI network of
Omicron and the score of the CM are visualizing in Figure 15 to 15. We have
calculated the maximum value of each centrality measure and divided it by two
to get each category’s threshold value. The threshold value will help us to
signify the important protein in the network. We have highlighted the
significant protein by getting the intersection of all the important proteins
of each category (CC, DC, and, PCR). A total of $8$ significant proteins were
detected from $68$ unique proteins. In our research work, we have extracted
the $18$ network cluster from the Omicron main network with the help of the
Markov clustering algorithm shown in Figure 15 and 15. In Table 3 we can see,
cluster $\mathcal{C}_{1},\mathcal{C}_{2}$ and $\mathcal{C}_{3}$ has $4$
protein in each, $\mathcal{C}_{4}$ to $\mathcal{C}_{8}$ has $4$ protein in
each and rest of the clusters are containing $2$ protein in each.

Table 1: Global properties of Omicron lineage variant’s network

Variants | #node | #edge | avg. node degree | max degree | density | avg. LCC | P-value
---|---|---|---|---|---|---|---
B.1.1.529 | 27 | 10 | 0.741 | 2 | 0.0284 | 0.407 | 0.3370
BA.1 | 34 | 18 | 1.06 | 3 | 0.0284 | 0.627 | 0.00365
BA.1.1 | 34 | 18 | 1.06 | 3 | 0.0284 | 0.627 | 0.00365
BA.2 | 38 | 25 | 1.32 | 4 | 0.03556 | 0.535 | 0.00038
BA.3 | 34 | 19 | 1.12 | 3 | 0.00416 | 0.657 | 0.00104
BA.4 | 31 | 17 | 1.1 | 3 | 0.06719 | 0.624 | 0.00271
BA.5 | 30 | 17 | 1.13 | 3 | 0.05666 | 0.711 | 0.00236
|  |  |  |  |  |  |

00footnotetext: avg. LCC: Average Local clustering coefficient, P-value: PPI
enrichment P-value

Table 2: Global properties of Omicron network

#Node | #Edges | Max Degree | Avg. node degree | Density | Avg. LCC | p-value
---|---|---|---|---|---|---
68 | 52 | 7 | 1.53 | 0.0228 | 0.385 | 0.0963
|  |  |  |  |  |

11footnotetext: Avg LCC: Avg. Local Clustering Coefficient, p-value: PPI
enrichment p-value

Table 3: Generated $18$ clusters from MCL algorithm cluster | gene count | protein names
---|---|---
$\mathcal{C}_{1}$ | 4 | KCNK1,KCNK16,KCNK17,KCNQ1
$\mathcal{C}_{2}$ | 4 | ERBB2,GRB7,HLA-DRB1,PTPN11
$\mathcal{C}_{3}$ | 4 | FABP6,NR1H4,SLC27A5,UGT1A3
$\mathcal{C}_{4}$ | 3 | PDPK1,RPS6KA3,YWHAH
$\mathcal{C}_{5}$ | 3 | HSF1,MX1,UBE2I
$\mathcal{C}_{6}$ | 3 | AHSG,BMPR2,ITGB7
$\mathcal{C}_{7}$ | 3 | SNRPB,SNRPD1,SNRPE
$\mathcal{C}_{8}$ | 3 | CHMP1B,IST1,PHF1
$\mathcal{C}_{9}$ | 2 | RPS26,RPSA
$\mathcal{C}_{10}$ | 2 | CLASP2,MAPT
$\mathcal{C}_{11}$ | 2 | COPB1,RAB2A
$\mathcal{C}_{12}$ | 2 | PSMC2,PSMD13
$\mathcal{C}_{13}$ | 2 | GTF2B,RARA
$\mathcal{C}_{14}$ | 2 | NDUFB5,NDUFV1
$\mathcal{C}_{15}$ | 2 | GPC1,IHH
$\mathcal{C}_{16}$ | 2 | CNNM2,NIPA1
$\mathcal{C}_{17}$ | 2 | AKR1B10,ALDH7A1
$\mathcal{C}_{18}$ | 2 | CFB,CFH
|  |

Figure 2: PPIN of B.1.1.529

Figure 3: PPIN of BA.1

Figure 4: PPIN of BA.2

Figure 5: PPIN of BA.3

Figure 6: PPIN of BA.4

Figure 7: PPIN of BA.5

Table 4: Centrality measure and some important score of 68 Omicron protein

S.N. | Protein Name | Node Degree | CCo | DC | CC | BC | PRC
---|---|---|---|---|---|---|---
1 | AHSG | 4 | 0 | 0.06 | 0.123 | 0 | 0.036
2 | AKR1B10 | 1 | 0 | 0.015 | 0.015 | 0 | 0.02
3 | ALDH7A1 | 1 | 0 | 0.015 | 0.015 | 0 | 0.02
4 | ASIC1 | 0 | 0 | 0 | 0 | 0 | 0.003
5 | BMPR2 | 2 | 0 | 0.03 | 0.104 | 0 | 0.02
6 | CFB | 2 | 0 | 0.03 | 0.115 | 0 | 0.018
7 | CFH | 2 | 0 | 0.03 | 0.121 | 0 | 0.018
8 | CHMP1B | 2 | 0 | 0.03 | 0.094 | 0 | 0.021
9 | CLASP2 | 1 | 0 | 0.015 | 0.108 | 0 | 0.01
10 | CNNM2 | 1 | 0 | 0.015 | 0.077 | 0 | 0.012
11 | COPB1 | 2 | 0 | 0.03 | 0.121 | 0 | 0.019
12 | DDC | 0 | 0 | 0 | 0 | 0 | 0.003
13 | EIF2B1 | 0 | 0 | 0 | 0 | 0 | 0.003
14 | EPHB4 | 0 | 0 | 0 | 0 | 0 | 0.003
15 | ERAP1 | 0 | 0 | 0 | 0 | 0 | 0.003
16 | ERBB2 | 3 | 0.333 | 0.045 | 0.137 | 0.333 | 0.023
17 | FABP6 | 1 | 0 | 0.015 | 0.114 | 0 | 0.01
18 | FAM20A | 0 | 0 | 0 | 0 | 0 | 0.003
19 | GNPNAT1 | 0 | 0 | 0 | 0 | 0 | 0.003
20 | GP1BA | 0 | 0 | 0 | 0 | 0 | 0.003
21 | GPC1 | 1 | 0 | 0.015 | 0.015 | 0 | 0.02
22 | GRB7 | 2 | 1 | 0.03 | 0.126 | 1 | 0.016
23 | GTF2B | 1 | 0 | 0.015 | 0.108 | 0 | 0.01
24 | HEXA | 0 | 0 | 0 | 0 | 0 | 0.003
25 | HLA-DRB1 | 2 | 0 | 0.03 | 0.132 | 0 | 0.017
26 | HSF1 | 1 | 0 | 0.015 | 0.124 | 0 | 0.01
27 | IHH | 1 | 0 | 0.015 | 0.015 | 0 | 0.02
28 | IST1 | 1 | 0 | 0.015 | 0.081 | 0 | 0.012
29 | ITGB7 | 1 | 0 | 0.015 | 0.102 | 0 | 0.011
30 | KCNK1 | 4 | 0.333 | 0.06 | 0.163 | 0.333 | 0.028
31 | KCNK16 | 3 | 0.667 | 0.045 | 0.147 | 0.667 | 0.021
32 | KCNK17 | 2 | 1 | 0.03 | 0.129 | 1 | 0.015
33 | KCNQ1 | 4 | 0.167 | 0.06 | 0.171 | 0.167 | 0.027
34 | KPNA2 | 0 | 0 | 0 | 0 | 0 | 0.003
35 | MAPT | 7 | 0.048 | 0.104 | 0.132 | 0.048 | 0.054
36 | MX1 | 1 | 0 | 0.015 | 0.124 | 0 | 0.01
37 | NDUFB5 | 2 | 1 | 0.03 | 0.109 | 1 | 0.017
38 | NDUFV1 | 2 | 1 | 0.03 | 0.109 | 1 | 0.017
39 | NIPA1 | 2 | 0 | 0.03 | 0.089 | 0 | 0.022
40 | NIPAL1 | 0 | 0 | 0 | 0 | 0 | 0.003
41 | NIPAL4 | 0 | 0 | 0 | 0 | 0 | 0.003
42 | NR1H4 | 5 | 0 | 0.075 | 0.14 | 0 | 0.044
43 | OBP2B | 0 | 0 | 0 | 0 | 0 | 0.003
44 | PDPK1 | 2 | 0 | 0.03 | 0.093 | 0 | 0.017
45 | PHF1 | 2 | 0 | 0.03 | 0.11 | 0 | 0.018
46 | PPP3CA | 0 | 0 | 0 | 0 | 0 | 0.003
47 | PSMC2 | 5 | 0.1 | 0.075 | 0.149 | 0.1 | 0.039
48 | PSMD13 | 1 | 0 | 0.015 | 0.119 | 0 | 0.01
49 | PTPN11 | 4 | 0.167 | 0.06 | 0.149 | 0.167 | 0.029
50 | RAB2A | 1 | 0 | 0.015 | 0.1 | 0 | 0.011
51 | RARA | 3 | 0 | 0.045 | 0.132 | 0 | 0.025
52 | RBM6 | 0 | 0 | 0 | 0 | 0 | 0.003
53 | RPS26 | 3 | 0.333 | 0.045 | 0.161 | 0.333 | 0.022
54 | RPS6KA3 | 2 | 0 | 0.03 | 0.109 | 0 | 0.017
55 | RPSA | 2 | 1 | 0.03 | 0.138 | 1 | 0.016
56 | SLC27A5 | 1 | 0 | 0.015 | 0.114 | 0 | 0.01
57 | SLC4A4 | 0 | 0 | 0 | 0 | 0 | 0.003
58 | SNRPB | 2 | 1 | 0.03 | 0.03 | 1 | 0.02
59 | SNRPD1 | 2 | 1 | 0.03 | 0.03 | 1 | 0.02
60 | SNRPE | 2 | 1 | 0.03 | 0.03 | 1 | 0.02
61 | TERF2IP | 0 | 0 | 0 | 0 | 0 | 0.003
62 | TNFRSF13C | 0 | 0 | 0 | 0 | 0 | 0.003
63 | TRPM1 | 0 | 0 | 0 | 0 | 0 | 0.003
64 | UBE2I | 5 | 0 | 0.075 | 0.157 | 0 | 0.04
65 | UGT1A3 | 1 | 0 | 0.015 | 0.114 | 0 | 0.01
66 | UGT2B7 | 0 | 0 | 0 | 0 | 0 | 0.003
67 | YWHAH | 2 | 0 | 0.03 | 0.109 | 0 | 0.017
68 | ZRANB3 | 0 | 0 | 0 | 0 | 0 | 0.003
|  |  |  |  |  |  |

11footnotetext:  Significant proteins:- AHSG: Alpha 2-HS Glycoprotein, KCNK1:
Potassium channel subfamily K member 1, KCNQ1: Potassium Voltage-Gated Channel
Subfamily Q Member 1, MAPT: Microtubule Associated Protein Tau, NR1H4: Nuclear
Receptor Subfamily 1 Group H Member 4, PSMC2: Proteasome 26S Subunit, ATPase
2, PTPN11: Tyrosine-protein phosphatase non-receptor type 11, UBE21: UBE2I
ubiquitin conjugating enzyme E2 I

Figure 8: protein vs. node degree

Figure 9: protein vs. CCo

Figure 10: protein vs. DC

Figure 11: protein vs. CC

Figure 12: protein vs. BC

Figure 13: protein vs. PR

Figure 14: Omnicron PPI network

Figure 15: PPIN after clustering

##  5 Conclusion

Centrality analysis are very useful for analyzing large biological networks.
Using a candidate gene network of Omicron as a case study, we investigated and
compared different centrality measures. According to the findings, it is
beneficial to explore candidate gene networks using methods from other fields
of science such as social network analysis. On the $7$ base lineage of Omicron
variations, including the $68$ unique protein encoded by the Omicron candidate
gene, graph analysis is done. From the Omicron main network, we extracted the
Markov clustering algorithm’s findings i.e., $18$ network clusters. The
primary Omicron network has $68$ nodes, each of which represents a protein. Of
the $68$ proteins, $8$ were found to be significant, including AHSG, KCNK1,
KCNQ1, MAPT, NR1H4, PSMC2, PTPN11, and UBE21, with the MAPT protein receiving
the highest score. The MAPT protein has the most dominating influence on the
protein-protein interaction network of the Omicron candidate gene, according
to the centrality score.Medical researchers as well as the general public will
benefit from this work, as it will allow them to to consider biological
knowledge in network analysis of the Omicron virus.

Analysis of networks can benefit greatly from centrality measures. However,
they are also required to be properly informed, selected, and applied. As part
of our main research work, we present information about the four major
centrality measures that have been found to be relevant for finding the most
significant proteins in the Omicron Lineage Variants PPIN. A wide range of new
and large networks are being created and developed due to different
applications and different centrality measures. The majority of studies have
tried to demonstrate the uniqueness and superiority of their centrality
measures. We still have a lot to learn about making a difference and applying
them properly. This is how we presented it.

  * •

Authors and Affiliations:

    * –

Mamata Das, Department of Computer Applications, NIT Trichy, India

    * –

Selvakumar K., Department of Computer Applications, NIT Trichy, India

    * –

P.J.A. Alphonse, Department of Computer Applications, NIT Trichy, India

    * –

Corresponding author: Mamata Das and Selvakumar K.

  * •

Compliance with Ethical Standards

    * –

Funding: This is the work of the first author under her doctoral. This
research received no external funding.

    * –

Disclosure of potential conflicts of interest: On behalf of all authors, the
corresponding author states that there is no conflict of interest.

    * –

Research involving human participants and/or animals: This article does not
contain any studies with human participants or animals performed by any of the
authors.

## References

  * [1] Arumugam Amala and Isaac Arnold Emerson.  Identification of target genes in cancer diseases using protein–protein interaction networks.  Network Modeling Analysis in Health Informatics and Bioinformatics, 8(1):1–13, 2019.
  * [2] Minoo Ashtiani and Ali Salehzadeh-Yazdi.  A systematic survey of centrality measures for protein-protein interaction networks.  BMC systems biology, 12(1):1–17, 2018.
  * [3] Laura F. Bringmann and Timon Elmer.  What do centrality measures measure in psychological networks?  Journal of abnormal psychology, 128(8):892, 2019.
  * [4] Sylvain Brohée and Jacques van Helden.  Evaluation of clustering algorithms for protein-protein interaction networks.  BMC bioinformatics, 7(1):1–19, 2006.
  * [5] Peter J. Carrington.  Crime and social network analysis.  The SAGE handbook of social network analysis, pages 236–255, 2011\.
  * [6] Kousik Das, Sovan Samanta, and Madhumangal Pal.  Study on centrality measures in social networks: a survey.  Social network analysis and mining, 8(1):1–11, 2018.
  * [7] Mamata Das, Pja Alphonse, and Selvakumar Kamalanathan.  Markov Clustering Algorithms and Their Application in Analysis of PPI Network of Malaria Genes.  In IDAACS, volume 2, pages 855–860. IEEE, 2021.
  * [8] Mamata Das, Pja Alphonse, and Selvakumar Kamalanathan.  An Analytical Study of COVID-19 Dataset Using Graph-Based Clustering Algorithms.  In Smart Intelligent Computing and Applications, volume 1, pages 1–15. Springer, 2022.
  * [9] Mohammad R. Diansyah and Wisnu A. Kusuma.  Analysis of protein-protein interaction using skyline query on parkinson disease.  In ICACSIS, pages 175–180. IEEE, 2019.
  * [10] Linton C. Freeman.  A set of measures of centrality based on betweenness.  Sociometry, pages 35–41, 1977.
  * [11] Mahdieh Ghasemi, Hossein Seidkhani, Faezeh Tamimi, and Rahgozar.  Centrality measures in biological networks.  Current Bioinformatics, 9(4):426–441, 2014.
  * [12] Daniel Z. Grunspan, Benjamin L. Wiggins, and Steven M. Goodreau.  Understanding classrooms through social network analysis: A primer for social network analysis in education research.  CBE—Life Sciences Education, 13(2):167–178, 2014.
  * [13] Gábor Iván and Vince Grolmusz.  When the Web meets the cell: using personalized PageRank for analyzing protein interaction networks.  Bioinformatics, 27(3):405–407, 2011.
  * [14] I. Jayaweera, K. Perera, and Jayantha Munasinghe.  Centrality measures to identify traffic congestion on road networks: A case study of sri lanka.  IOSR Journal of Mathematics (IOSRJM), 2017.
  * [15] Karen E. Joyce, Paul J. Laurienti, and Jonathan H. Burdette.  A new measure of centrality for brain networks.  PloS one, 5(8):e12200, 2010.
  * [16] Hakimeh Khojasteh, Alireza Khanteymoori, and Mohammad Hossein Olyaee.  Comparing protein–protein interaction networks of SARS-CoV-2 and (H1N1) influenza using topological features.  Scientific reports, 12(1):1–11, 2022.
  * [17] Jesse Laeuchlia and Yunior Ramírez-Cruzb.  Analysis of centrality measures under differential privacy models.  Applied Mathematics and Computation, 412:126546, 2022.
  * [18] Andrea Landherr, Bettina Friedl, and Julia Heidemann.  A critical review of centrality measures in social networks.  Business & Information Systems Engineering, 2(6):371–385, 2010\.
  * [19] World Health Organization. Office of Library and Health Literature Services.  Styles for bibliographic citations : guidelines for WHO-produced bibliographies.  ONLINE, 1988.
  * [20] Gert Sabidussi.  The centrality index of a graph.  Psychometrika, 31(4):581–603, 1966.
  * [21] Venu Satuluri and Srinivasan Parthasarathy.  Scalable graph clustering using stochastic flows: applications to community discovery.  In Proceedings of the 15th ACM SIGKDD international conference on Knowledge discovery and data mining, pages 737–746, 2009.
  * [22] Malcolm K. Sparrow.  The application of network analysis to criminal intelligence: An assessment of the prospects.  Social networks, 13(3):251–274, 1991.
  * [23] Damian Szklarczyk, Andrea Franceschini, Stefan Wyder, et al.  STRING v10: protein–protein interaction networks, integrated over the tree of life.  Nucleic acids research, 43(D1):D447–D452, 2015.
  * [24] The UniProt Consortium.  UniProt: the universal protein knowledgebase in 2021.  Nucleic acids research, 49(D1):D480–D489, 2021.
  * [25] Stijn vanDongen.  A cluster algorithm for graphs.  Information Systems, 2000.
  * [26] Junyi Wang, Xiaoni Hou, Kezan Li, and Yong Dinga.  A novel weight neighborhood centrality algorithm for identifying influential spreaders in complex networks.  Physica A: Statistical Mechanics and its Applications, 475:88–105, 2017.