# Detection of Global Anomalies on Distributed IoT Edges
with Device-to-Device Communication

Hideya Ochiai  The University of TokyoJapan<EMAIL_ADDRESS>0000-0002-4568-6726 ,  Riku Nishihata  The University of TokyoJapan nishihata-
<EMAIL_ADDRESS>,  Eisuke Tomiyama  The University of TokyoJapan
<EMAIL_ADDRESS>,  Yuwei Sun  The University of TokyoJapan
<EMAIL_ADDRESS>and  Hiroshi Esaki  The University of TokyoJapan
<EMAIL_ADDRESS>0000-0001-5657-9216

(2023)

###### Abstract.

Anomaly detection is an important function in IoT applications for finding
outliers caused by abnormal events. Anomaly detection sometimes comes with
high-frequency data sampling which should be carried out at Edge devices
rather than Cloud. In this paper, we consider the case that multiple IoT
devices are installed in a single remote site and that they collaboratively
detect anomalies from the observations with device-to-device communications.
For this, we propose a fully distributed collaborative scheme for training
distributed anomaly detectors with Wireless Ad Hoc Federated Learning, namely
“WAFL-Autoencoder”. We introduce the concept of Global Anomaly which sample is
not only rare to the local device but rare to all the devices in the target
domain. We also propose a distributed threshold-finding algorithm for Global
Anomaly detection. With our standard benchmark-based evaluation, we have
confirmed that our scheme trained anomaly detectors perfectly across the
devices. We have also confirmed that the devices collaboratively found
thresholds for Global Anomaly detection with low false positive rates while
achieving high true positive rates with few exceptions.

Anomaly Detection, Collaborative Learning, Device-to-Device Communication, The
Internet of Things

††journalyear: 2023††conference: The Twenty-fourth International Symposium on
Theory, Algorithmic Foundations, and Protocol Design for Mobile Networks and
Mobile Computing; October 23–26, 2023; Washington, DC, USA (Preprint
version)††booktitle: The Twenty-fourth International Symposium on Theory,
Algorithmic Foundations, and Protocol Design for Mobile Networks and Mobile
Computing (MobiHoc ’23), October 23–26, 2023, Washington, DC, USA††copyright:
acmcopyright††copyright: none††booktitle:  ††price:  ††doi:  ††isbn:  ††ccs:
Computing methodologies Anomaly detection††ccs: Networks Ad hoc networks††ccs:
Hardware Wireless integrated network sensors

##  1\. Introduction

Anomaly detection is an important function in Internet of Things (IoT)
applications for finding outliers caused by abnormal events such as machine
faults, electric surges, or security incidents. Sensors such as cameras,
accelerometers, or electric meters can be used for detecting anomalies in
remote sites. Here, we must remember that the data generated by these sensors
are relatively large because of their sampling frequencies necessary for
anomaly detection.

Uploading all the data to the Cloud is cost-ineffective because most of the
data are normal and the normal data consume both Cloud storage and
communication bandwidth, increasing the cost of their business. Edge IoT
computing is a promising approach for finding anomalies at remote sites from
highly-frequent sampling data and for picking up anomaly cases only when they
were detected.

In this paper, we consider the case that multiple IoT devices are installed in
a single remote site, such as in a building, and that they collaboratively
detect anomalies from the observations. They are expected (1) to learn the
normal features from daily observations and (2) to detect anomalies when they
occurred.

In such a case, because devices are located physically nearby, we can utilize
device-to-device communications for training anomaly detectors. Conventional
federated learning (Konečnỳ et al., 2016; Li et al., 2020b) can be also used
but all the IoT devices need to exchange model parameters with the Cloud
repeatedly, which multiplies the uplink (i.e., cellular) traffic.

We propose WAFL-Autoencoder – a fully distributed collaborative learning
scheme for autoencoder utilizing device-to-device communications. Here, WAFL
stands for wireless ad hoc federated learning (Ochiai et al., 2022). Even
though the name has federated learning, the architecture is fully distributed
and different from conventional federated learning. In our previous study, we
confirmed WAFL works efficiently for supervised cases. In this paper, we focus
on an unsupervised case and study how WAFL allows training for distributed
autoencoders.

We also try to use the trained autoencoders for distributed anomaly detection
applications. Here, we introduce the concept of Global Anomaly and Local
Anomaly as Fig. 1. A local anomaly is something rare to the device, but
familiar to others. For example, as Fig. 1, in device A, the major data are
handwritten character 0, and the handwritten character 2 is a local anomaly
because the 2 is rare to device A, but a major sample at another node. The
local anomaly can be easily detected without collaborating with other devices.

Please note that the shirt image in device D is an abnormal image for all the
devices. This is what we call Global Anomaly. Detection of a global anomaly is
challenging because the ML model needs to know that the event is also rare to
others without sharing the data itself. We focus on the detection of global
anomalies as an application of WAFL-Autoencoder. This raises another technical
challenge of threshold finding.

A threshold for anomaly detection is usually calculated from the validation
data that the device has. However, the calculated threshold in this way is
usually skewed for the device’s local data. For this problem, we also propose
a threshold-finding algorithm in this paper.

Figure 1. Local anomaly and global anomaly examples. Local anomaly is
something rare to the device, but well-known to others. Global anomaly is rare
to every device, which is difficult to detect with low false-positive rates.

Because this is a work-in-progress paper, we focus on understanding the novel
characteristics of WAFL-Autoencoder with standard benchmark datasets. The
evaluation includes the capability of reconstruction of input samples and the
performances of global anomaly detection with pure and dirty cases.

In summary,

  * •

We propose a WAFL-Autoencoder that utilizes device-to-device communication for
learning the features of distributed data samples at IoT edge devices.

  * •

We introduce the concept of Global Anomaly and Local Anomaly for distributed
cases, and focus on Global Anomaly detection, which is a challenging topic.

  * •

We propose a distributed threshold-finding algorithm for finding Global
Anomaly, which should be carried out along with the training of WAFL-
Autoencoder.

  * •

This paper focuses on benchmark-based evaluation with standard MNIST families
for understanding the novel characteristics of WAFL-Autoencoder and Global
Anomaly detection.

This paper is organized as follows. In section 2, we describe our related
work. In section 3, we propose WAFL-Autoencoder and threshold finding
algorithm, introducing the concept of global anomaly. In section 4, we provide
our evaluation. We conclude this paper in Section 5.

##  2\. Related Work

There are two types of federated learning: i.e., centralized federated
learning (CFL) and decentralized federated learning (DFL). CFL is a well-known
architecture, which has a parameter server for model sharing and aggregation.
DFL is a new architecture that does not have any central services.

In terms of anomaly detection, we can find some CFL-based studies in computer
networks (Pei et al., 2022; Zhao et al., 2019), IoT applications (Sater and
Hamza, 2021; Weinger et al., 2022; Cui et al., 2021), and so forth. Within
this architecture, variational autoencoder (Zhang et al., 2021; Palato, 2021),
GAN (Tabassum et al., 2022; Akcay et al., 2018), One-Class SVM (Anaissi et
al., 2022) were applied and studied. These systems would be useful if we can
set up a centralized server by ourselves, or if we can trust the third-party
service provider. However, they still consume the global communication traffic
for repeated model exchanges.

The architecture of DFL is studied with several peer-to-peer system
architectures. Blockchain-based systems (Li et al., 2020a; Bao et al., 2019;
Pokhrel and Choi, 2020; ur Rehman et al., 2020; Warnat-Herresthal et al.,
2021) were proposed for updating global models without any centralized
mechanisms. (Hu et al., 2019; Pappas et al., 2021; Roy et al., 2019) were
proposed with other forms of peer-to-peer systems. However, those studies
assume the Internet connection, where all the devices need to know remote
global addresses with each other. They do not investigate global anomaly
detection, which is an important topic of this paper. Besides, usually peer-
to-peer systems increase global network traffic.

Wireless ad hoc federated learning (WAFL) focuses more on opportunistic and ad
hoc scenarios where the communication among the devices happens only when they
are within the radio range(Ochiai et al., 2022). The idea is based on
communication without infrastructure, which is named “wireless ad hoc network
(WANET)” (Frodigh et al., 2000). WANET is one of the famous domains in
computer network research. The family of WANET include mobile ad hoc network
(Abolhasan et al., 2004), vehicular ad hoc network (Hartenstein and
Laberteaux, 2008), and delay/disruption tolerant network (Fall, 2003). They
allow communications among today’s smart devices with device-to-device
communication, i.e., without third-party intervention. The legacy ad hoc
networks was for communication not for collaborative machine learning as WAFL
proposed.

In this paper, for finding global anomalies in a distributed scenario, we have
taken an approach of distributed autoencoders that learn the features
collaboratively from distributed sources, which is new in these research
fields. We have previously studied the basic learning paradigm in (Ochiai et
al., 2022), and published the application to distributed generative
adversarial network(Tomiyama et al., 2023). This paper focus on autoencoder
for learning the features and finding anomaly in such a distributed scenario.

##  3\. WAFL-Autoencoder for Anomaly Detection

To learn the features of data in a distributed environment, we take the
approach of distributed autoencoders with distributed anomaly thresholds.

###  3.1. Definition of Terms

We assume a set of mobile devices $N$ that collaboratively train models based
on their local data. Let $W^{n}$ be the model parameters of device $n\in N$,
and $X^{n}$ be a set of data samples in $n$. With function $AE$, the output of
autoencoder $\hat{x}$ to an input $x\in X^{n}$ at device $n$, can be
formulated as:

(1) $\hat{x}=AE(x,W^{n})$

With appropriate loss function, such as mean squared error, i.e.,
$MSE(\hat{x},x)$, $W^{n}$ can be trained to fit specifically for $X^{n}$.

###  3.2. Global Anomaly

In the case of distributed environment, we have at least two types of
anomalies as we discussed: Local Anomaly and Global Anomaly. A local anomaly
is a rare data sample to the device. It might be important to the device, but
it is still “legitimate” in general. A global anomaly is rare to everyone,
which is more important to be detected. Thus, we focus on the detection of
global anomalies in this paper.

####  3.2.1. Implicit Data Category and Class

To formulate a global anomaly, we associate data $x$ with implicit data
category and class. We consider that data sample $x$ is implicitly associated
with a class $y$. In other words, $x$ could be paired with $y$ as ($x$,$y$).
Here, class $y$ belongs to a category denoted by $Y$: $y\in Y$. In this paper,
we consider category $Y_{0}$, e.g., MNIST, as a world of the legitimate
dataset. Other categories such as $Y_{1}$ and $Y_{2}$, e.g., Noisy-MNIST and
Occluded-MNIST respectively, are the global anomaly domains.

Please note that $y$ and $Y$ are given implicitly and not used for training
autoencoder, but we consider there can be such data structure behind $X$. Only
the collection of $x$ is given as a dataset for training, and they are not
explicitly associated with classes or categories in practice.

####  3.2.2. Distributions of Train, Validation, and Test Data

It is common to assume three types of data: i.e., (1) train, (2) validation,
and (3) test data. In this study, we use train data for training the model
parameters of the autoencoder. We use validation data for finding a threshold
for anomaly detection. We use test data for evaluating the performances of the
developed models and anomaly detection. Thus, train and validation data are
possessed in the device’s local storage, but test data is not stored as an
assumption.

We consider two cases for local data: i.e., pure and dirty cases. In the pure
case, we consider that all the train and validation data are legitimate, i.e.,
implicitly associated with category $Y_{0}$ only. In the dirty case, we
consider that they also contain a small number of global anomaly samples from
$Y_{i}(i\neq 0$).

Test data should be derived from both legitimate category $Y_{0}$ and anomaly
categories $Y_{i}(i\neq 0)$.

###  3.3. Aggregation of Distributed Models with Device-to-Device
Communication

Let $nbr(n)$ be a set of neighbors of device $n$ at a certain time. Here,
$nbr(n)$ does not include itself, i.e., $n\notin nbr(n)$. $nbr(n)$ may
dynamically change based on the physical mobility of the devices. Device $n$
can directly communicate with all of $nbr(n)$ over wireless channels such as
Wi-Fi (ad hoc mode) or Bluetooth.

Device $n$ is expected to receive all the model parameters from $nbr(n)$.
Then, device $n$ aggregates the models with the local model $W^{n}$ by the
following formula,

(2) $W^{n}\leftarrow W^{n}+\lambda\frac{\sum_{k\in
nbr(n)}{(W^{k}-W^{n})}}{|nbr(n)|+1}.$

Here, $\lambda$ is a coefficient parameter for adding the differences between
the neighbors. The aggregated model needs additional mini-batch-based training
with its local data.

After interacting with many devices, $W^{n}$ becomes generalized and it will
be able to correctly predict the data samples that device $n$ does not have
but others have. Please note that during this process, WAFL-AE achieves model
training without relying on any centralized mechanisms or exchanging the local
data itself.

###  3.4. Finding Anomaly Threshold with Device-to-Device Communication

Because of the Non-IID characteristics of distributed environment, a
calculated anomaly threshold from local validation data is also bound to the
distributions of the local dataset. The thresholds of anomaly scores
calculated in this way can be different by devices. For this issue, we propose
aggregation of those calculated thresholds to find the proper threshold. This
improves the performance of anomaly detection in a distributed environment.

Let $s(x,W^{n})$ be an anomaly score of sample $x$ with model parameter
$W^{n}$ at device $n$. According to the work (Martin, 2021), there are many
design choices for anomaly score $s$. As an example, we use the following
anomaly score $s$ throughout the paper:

(3) $s(x,W^{n})=\frac{MSE(AE(x,W^{n}),x)}{density(x)},$ (4)
$density(x)=\frac{\sum_{i\in x}p_{i}}{|x|}.$

Here, $i$ is the index of a feature of $x$. $p_{i}\in[0,1]$ is the normalized
value of the feature. $|x|$ is the number of features. This anomaly score is
based on the common MSE-based method with normalization of the input sample.

With a threshold of anomaly detection $\alpha^{n}$, device $n$ identifies that
$x$ is anomaly if $s(x,W^{n})>\alpha^{n}$, otherwise $x$ is not anomaly. If
$W^{n}$ is generalized enough for legitimate samples over the distributed
data, the detected anomaly is a global anomaly.

To find proper $\alpha^{n}$, we propose to exchange and aggregate $\alpha^{n}$
among the encounter devices during the model training phase. This can be
formulated as follows:

(5) $\alpha^{n}\leftarrow\frac{\sum_{k\in
nbr(n)}{\alpha^{k}+\alpha^{n}+\gamma\beta^{n}}}{|nbr(n)|+1+\gamma}.$

This calculates the average of $\alpha$ among the devices encountered.
$\beta^{n}$ – a locally calculated threshold with local validation data at
device $n$ is also introduced for averaging (with coefficient $\gamma$).

We calculate $\beta^{n}$ at every time after the model update in the following
manner. Given a threshold rate $\delta\in[0,1]$, for example, $\delta=0.999$,
the device finds such $s(x,W^{n})$ which cumulative distribution function
(CDF) reaches $\delta$ as $\beta^{n}$, i.e.,
$\beta^{n}=CDF^{-1}\left(\delta;s(X^{n}_{\text{val}},W^{n})\right)$. Here,
$X^{n}_{\text{val}}$ is the validation dataset of device $n$.

##  4\. Evaluation

To understand the basic characteristics of the WAFL-Autoencoder, i.e., the
ability to learn the features, reconstruct legitimate data, and detect global
anomalies, we conducted a benchmark-based evaluation as follows.

###  4.1. Experiment Setting

For training, we have configured 10 devices ($n=0,1,\ldots,9$) to have 99.95%
Non-IID MNIST dataset for train and validation. MNIST data samples were
scattered so that each label $n$ sample was assigned to device $n$ with a
probability of 99.95% and other devices uniformly for the rest of 0.05%. The
ratio between train and validation was 4:1.

For testing, we have used the original MNIST test dataset, which is IID, as
legitimate test data. To evaluate the detection of global anomalies, we have
prepared four types of test datasets: (1) Noisy-MNIST, (2) Occluded-MNIST, (3)
Fashion-MNIST, and (4) Kuzushiji-MNIST (Clanuwat et al., 2018). We generated
Noisy-MNIST by adding periodic (every 10 pixels) white spots to the original
MNIST test images, and Occluded-MNIST by putting a 5-pixel diameter dark
circle at position (14,20) of the original MNIST test images.

For the evaluation of the dirty case discussed in Section 3.2.2, we have mixed
global anomaly samples in the training phase by swapping 50 MNIST samples with
50 Fashion-MNIST samples. In this configuration, the ratio of global anomaly
samples in the training dataset is about 1%.

Table 1. Configuration of the autoencoder used in the evaluation. $(i,o,k,p,op)$ stands for ($input$, $output$, $kernel\\_size$, $padding$, $output\\_padding$). We have used $stride=2$ in Conv2D and TransConv2D. Comp. | Layer (Parameters)
---|---
Encoder | Conv2D ($i$=1, $o$=8, $k$=3, $p$=1) - ReLU -
| Conv2D ($i$=8, $o$=16, $k$=3, $p$=1) -
| BatchNorm2d (16) - ReLU -
| Conv2D ($i$=16, $o$=32, $k$=3, $p$=0) - Flatten -
| ReLU - Linear ($i$=288, $o$=128) -
| ReLU - Linear ($i$=128, $o$=64)
Decoder | Linear ($i$=64, $o$=128) - ReLU -
| Linear ($i$=128, $o$=288) - ReLU - Unflatten -
| TransConv2D ($i$=32, $o$=16, $k$=3, $p$=0, $op$=0) -
| BatchNorm2d (16) - ReLU -
| TransConv2D ($i$=16, $o$=8, $k$=3, $p$=1, $op$=1) -
| BatchNorm2d (8) - ReLU -
| TransConv2D ($i$=8, $o$=1, $k$=3, $p$=1, $op$=1) -
| Sigmoid

We have used random waypoint mobility to generate the contacts of devices,
which is very common in the studies of WANET (Bettstetter et al., 2004). We
have used RWP0500 – the most representative mobility pattern in (Ochiai et
al., 2022). According to the work, there are no huge differences in the
performance among other mobility patterns. Of course, there can be a
discussion of some special cases that a device is isolated from all other
devices all the time. But, such physical failure cases are not our target.

(a) Legitimate Inputs

(b) Global Anomaly Inputs

Figure 2. Reconstructed images of (a) legitimate inputs and (b) global anomaly
inputs at epoch = 0, 100, 500, and 5000 (at device 3). The model could
reconstruct minor legitimate class images successfully while keeping global
anomaly outputs unconstructed.

Table 1 shows the detailed configuration of the autoencoder used in our
experiment. For aggregation, model parameters except for batch normalization
layers were exchanged among the encountered devices.

We have used an SGD optimizer with momentum $=$ 0.9. Other hyper-parameters
were: learning rate $=$ 0.001, batch size $=$ 32, coefficient of WAFL-
aggregation $=$ 0.1, coefficient of threshold aggregation $\gamma=0.01$, and
threshold ratio $\delta=0.999$.

###  4.2. Reconstructed Images

Fig. 2 shows the reconstructed examples of (a) MNIST and (b) Fashion-MNIST
test samples. At the early stage, e.g., at epochs 0 and 100, it could not
reconstruct the images well except for the major classes of training data,
which is ’3’ as observed.

After several interactions among the devices, i.e., after epoch 500, they
could reconstruct the images very well for all the MNIST classes in Fig. 2(a).
They did not reconstruct the images in the case of Fashion-MNIST in Fig. 2(b),
which is successful for global anomaly detection.

This result indicates that general legitimate features were trained
successfully across the partitioned data even with the severe 99.95% Non-IID
case whereas global anomalies can be detected by comparing the inputs and the
reconstructions.

###  4.3. Detection of Global Anomalies

Fig. 3 shows the trend of positive rates for the detection by class or by
category of the test data from epoch 0 to epoch 5000 at device 0 in the pure
case. The number after MNIST (e.g., MNIST-0 $\ldots$ MNIST-9) indicates the
class. MNIST(ALL) means all the classes are mixed together. As MNIST is
legitimate in the global anomaly context, the rates of MNIST are false
positive rates (FPR). For Noisy-, Occluded-, Fashion-, and Kuzhushiji-MNIST
are the intentional global anomaly samples, the rates of them are true
positive rates (TPR). In this experiment, we have used three instances of RWP
mobilities with the same hyper-parameters but with different random seeds. The
FPRs and TPRs are the averages of them.

Figure 3. Positive rates of anomaly detection to the test data at device 0 for the models from epoch 0 to 5000. MNIST are false positive rates (FPR). Noisy-, Occluded-, Fashion-, and Kuzhushiji-MNIST are true positive rates (TPR). The number $y$ in MNIST-$y$ indicates FPR for class $y$ test samples.  Table 2. True/False positive rates by the developed models and thresholds at epoch 5000 per node. N, O, F, and K indicate Noisy-, Occluded-, Fashion-, and Kuzushiji-MNIST. Higher TPR with lower FPR is better. Self-train fails in achieving low FPRs. Device | WAFL-AE (train without anomaly) | WAFL-AE (train with 1% anomaly) | Self-train (without anomaly)
---|---|---|---
TPR | FPR | TPR | FPR | TPR | FPR
N | O | F | K | N | O | F | K | N | O | F | K
0 | 100% | 16.2% | 99.5% | 99.5% | 0.68% | 100% | 15.9% | 92.3% | 99.1% | 0.54% | 100% | 84.3% | 98.1% | 99.6% | 42.2%
1 | 100% | 5.75% | 100% | 100% | 3.66% | 100% | 4.53% | 100% | 100% | 3.27% | 100% | 80.7% | 98.7% | 99.6% | 63.5%
2 | 100% | 5.21% | 99.5% | 99.6% | 0.33% | 100% | 5.16% | 97.3% | 99.7% | 0.32% | 100% | 25.3% | 97.1% | 97.7% | 3.06%
3 | 100% | 5.60% | 99.6% | 99.5% | 0.40% | 100% | 5.94% | 97.9% | 99.7% | 0.43% | 100% | 18.3% | 97.8% | 96.4% | 2.50%
4 | 100% | 3.09% | 99.8% | 99.8% | 0.39% | 100% | 3.00% | 99.6% | 99.9% | 0.37% | 100% | 39.6% | 99.6% | 99.7% | 22.1%
5 | 100% | 2.82% | 99.7% | 99.7% | 0.33% | 100% | 2.74% | 99.6% | 99.9% | 0.32% | 100% | 20.9% | 99.4% | 99.1% | 3.54%
6 | 100% | 9.97% | 99.8% | 99.8% | 1.22% | 100% | 10.8% | 98.7% | 99.9% | 0.97% | 99.8% | 52.9% | 94.1% | 98.8% | 24.4%
7 | 100% | 2.43% | 99.9% | 99.9% | 0.48% | 100% | 2.36% | 99.3% | 99.9% | 0.46% | 100% | 41.8% | 96.7% | 98.4% | 16.2%
8 | 100% | 8.97% | 99.7% | 99.6% | 0.38% | 100% | 7.75% | 98.0% | 99.8% | 0.38% | 99.9% | 9.01% | 90.2% | 94.4% | 1.72%
9 | 100% | 3.69% | 99.9% | 99.9% | 0.40% | 100% | 3.27% | 99.4% | 99.9% | 0.36% | 100% | 24.9% | 99.3% | 98.8% | 8.58%
Avg. | 100% | 6.37% | 99.7% | 99.7% | 0.83% | 100% | 6.12% | 98.2% | 99.8% | 0.74% | 100% | 40.0% | 97.1% | 98.2% | 18.8%

Please note that the major training data at device 0 are MNIST class 0, and it
achieved good (i.e., low) FPR for MNIST-0 all the time even from epoch 0.
Other MNIST classes had much higher FPR at the beginning, but it almost
reached zero from epoch 200. This indicates that the WAFL-Autoencoder has
worked effectively for learning other legitimate data samples beyond the major
class of the device.

As for global anomalies, device 0 started to detect from epoch 500, finally
reaching the TPRs of 100%, 16%, 99%, and 100% for Noisy-, Occluded-, Fashion-,
and Kuzushiji-MNIST respectively at 5000.

We have observed two steps of stabilization for anomaly detection. The first
step was model stabilization. It took the first 1000 epochs for generalizing
the autoencoder model itself. The next step was threshold stabilization. After
stabilizing the model, the threshold could be stabilized.

Table 2 shows the summary of True / False positive rates by the developed
models and thresholds per device at epoch 5000 under (1) train without
anomaly, (2) train with 1% anomaly samples, and (3) self-training without
anomaly as a baseline. Compared to the self-training, WAFL-Autoencoder
achieved good performance on average showing lower false positive rates and
higher positive rates except for Occluded-MNIST in both pure and dirty cases.

##  5\. Conclusion

We have proposed the WAFL-Autoencoder that trains distributed autoencoders
with device-to-device communications assuming the application of multiple IoT
devices physically located nearby. We have also proposed a distributed
threshold-finding algorithm for finding Global anomalies with WAFL-
Autoencoder.

With our standard benchmark-based evaluation, we have confirmed that WAFL-
Autoencoder trains autoencoders perfectly for reconstructing the inputs of
legitimate data across the devices. We have also confirmed that the devices
collaboratively found thresholds for global anomaly detection with low false
positive rates while achieving high true positive rates for global anomalies
with few exceptions.

Future work may include the application of more realistic data such as
electric monitoring sequences, accelerometer signals, and photo images.

## Acknowledgement

This paper is based on the results obtained from a project commissioned by the
New Energy and Industrial Technology Development Organization (NEDO), Japan.
This work was also supported by JSPS KAKENHI Grant Number JP 22H03572.

## References

  * (1)
  * Abolhasan et al. (2004) Mehran Abolhasan, Tadeusz Wysocki, and Eryk Dutkiewicz. 2004\.  A review of routing protocols for mobile ad hoc networks.  _Ad hoc networks_ 2, 1 (2004), 1–22.
  * Akcay et al. (2018) Samet Akcay, Amir Atapour-Abarghouei, and Toby P Breckon. 2018\.  Ganomaly: Semi-supervised anomaly detection via adversarial training. In _Asian conference on computer vision_. Springer, 622–637.
  * Anaissi et al. (2022) Ali Anaissi, Basem Suleiman, and Widad Alyassine. 2022\.  A Personalized Federated Learning Algorithm for One-Class Support Vector Machine: An Application in Anomaly Detection. In _International Conference on Computational Science_. Springer, 373–379.
  * Bao et al. (2019) Xianglin Bao, Cheng Su, Yan Xiong, Wenchao Huang, and Yifei Hu. 2019.  FLchain: A blockchain for auditable federated learning with trust and incentive. In _2019 5th International Conference on Big Data Computing and Communications (BIGCOM)_. IEEE, 151–159.
  * Bettstetter et al. (2004) Christian Bettstetter, Hannes Hartenstein, and Xavier Pérez-Costa. 2004.  Stochastic properties of the random waypoint mobility model.  _Wireless networks_ 10, 5 (2004), 555–567.
  * Clanuwat et al. (2018) Tarin Clanuwat, Mikel Bober-Irizar, Asanobu Kitamoto, Alex Lamb, Kazuaki Yamamoto, and David Ha. 2018\.  Deep learning for classical japanese literature.  _arXiv preprint arXiv:1812.01718_ (2018).
  * Cui et al. (2021) Lei Cui, Youyang Qu, Gang Xie, Deze Zeng, Ruidong Li, Shigen Shen, and Shui Yu. 2021.  Security and privacy-enhanced federated learning for anomaly detection in iot infrastructures.  _IEEE Transactions on Industrial Informatics_ 18, 5 (2021), 3492–3500.
  * Fall (2003) Kevin Fall. 2003\.  A delay-tolerant network architecture for challenged Internets. In _Proceedings of the 2003 conference on Applications, technologies, architectures, and protocols for computer communications_. 27–34.
  * Frodigh et al. (2000) Magnus Frodigh, Per Johansson, and Peter Larsson. 2000\.  Wireless ad hoc networking: the art of networking without a network.  _Ericsson review_ 4, 4 (2000), 249.
  * Hartenstein and Laberteaux (2008) Hannes Hartenstein and LP Laberteaux. 2008.  A tutorial survey on vehicular ad hoc networks.  _IEEE Communications magazine_ 46, 6 (2008), 164–171.
  * Hu et al. (2019) Chenghao Hu, Jingyan Jiang, and Zhi Wang. 2019.  Decentralized federated learning: A segmented gossip approach.  _arXiv preprint arXiv:1908.07782_ (2019).
  * Konečnỳ et al. (2016) Jakub Konečnỳ, H Brendan McMahan, Felix X Yu, Peter Richtárik, Ananda Theertha Suresh, and Dave Bacon. 2016\.  Federated learning: Strategies for improving communication efficiency.  _arXiv preprint arXiv:1610.05492_ (2016).
  * Li et al. (2020b) Tian Li, Anit Kumar Sahu, Ameet Talwalkar, and Virginia Smith. 2020b.  Federated learning: Challenges, methods, and future directions.  _IEEE Signal Processing Magazine_ 37, 3 (2020), 50–60.
  * Li et al. (2020a) Yuzheng Li, Chuan Chen, Nan Liu, Huawei Huang, Zibin Zheng, and Qiang Yan. 2020a.  A blockchain-based decentralized federated learning framework with committee consensus.  _IEEE Network_ 35, 1 (2020), 234–241.
  * Martin (2021) Damien W Martin. 2021\.  _Fault detection in manufacturing equipment using unsupervised deep learning_.  Ph. D. Dissertation. Massachusetts Institute of Technology.
  * Ochiai et al. (2022) Hideya Ochiai, Yuwei Sun, Qingzhe Jin, Nattanon Wongwiwatchai, and Hiroshi Esaki. 2022.  Wireless Ad Hoc Federated Learning: A Fully Distributed Cooperative Machine Learning.  _arXiv preprint arXiv:2205.11779_ (2022).
  * Palato (2021) Mirko Palato. 2021\.  Federated variational autoencoder for collaborative filtering. In _2021 International Joint Conference on Neural Networks (IJCNN)_. IEEE, 1–8.
  * Pappas et al. (2021) Christodoulos Pappas, Dimitris Chatzopoulos, Spyros Lalis, and Manolis Vavalis. 2021.  IPLS: A framework for decentralized federated learning. In _2021 IFIP Networking Conference_. IEEE, 1–6.
  * Pei et al. (2022) Jiaming Pei, Kaiyang Zhong, Mian Ahmad Jan, and Jinhai Li. 2022\.  Personalized federated learning framework for network traffic anomaly detection.  _Computer Networks_ 209 (2022), 108906.
  * Pokhrel and Choi (2020) Shiva Raj Pokhrel and Jinho Choi. 2020.  Federated learning with blockchain for autonomous vehicles: Analysis and design challenges.  _IEEE Transactions on Communications_ 68, 8 (2020), 4734–4746.
  * Roy et al. (2019) Abhijit Guha Roy, Shayan Siddiqui, Sebastian Pölsterl, Nassir Navab, and Christian Wachinger. 2019\.  BrainTorrent: A peer-to-peer environment for decentralized federated learning.  _arXiv preprint arXiv:1905.06731_ (2019).
  * Sater and Hamza (2021) Raed Abdel Sater and A Ben Hamza. 2021.  A federated learning approach to anomaly detection in smart buildings.  _ACM Transactions on Internet of Things_ 2, 4 (2021), 1–23.
  * Tabassum et al. (2022) Aliya Tabassum, Aiman Erbad, Wadha Lebda, Amr Mohamed, and Mohsen Guizani. 2022.  FEDGAN-IDS: Privacy-preserving IDS using GAN and Federated Learning.  _Computer Communications_ (2022).
  * Tomiyama et al. (2023) Eisuke Tomiyama, Hiroshi Esaki, and Hideya Ochiai. 2023\.  WAFL-GAN: Wireless Ad Hoc Federated Learning for Distributed Generative Adversarial Networks. In _IEEE International Conference on Knowledge and Smart Technology_.
  * ur Rehman et al. (2020) Muhammad Habib ur Rehman, Khaled Salah, Ernesto Damiani, and Davor Svetinovic. 2020.  Towards blockchain-based reputation-aware federated learning. In _IEEE Conference on Computer Communications Workshops (INFOCOM WKSHPS)_. IEEE, 183–188.
  * Warnat-Herresthal et al. (2021) Stefanie Warnat-Herresthal, Hartmut Schultze, Krishnaprasad Lingadahalli Shastry, Sathyanarayanan Manamohan, et al. 2021\.  Swarm learning for decentralized and confidential clinical machine learning.  _Nature_ 594, 7862 (2021), 265–270.
  * Weinger et al. (2022) Brett Weinger, Jinoh Kim, Alex Sim, Makiya Nakashima, Nour Moustafa, and K John Wu. 2022\.  Enhancing IoT anomaly detection performance for federated learning.  _Digital Communications and Networks_ (2022).
  * Zhang et al. (2021) Kai Zhang, Yushan Jiang, Lee Seversky, Chengtao Xu, Dahai Liu, and Houbing Song. 2021\.  Federated variational learning for anomaly detection in multivariate time series. In _2021 IEEE International Performance, Computing, and Communications Conference (IPCCC)_. IEEE, 1–9.
  * Zhao et al. (2019) Ying Zhao, Junjun Chen, Di Wu, Jian Teng, and Shui Yu. 2019.  Multi-task network anomaly detection using federated learning. In _Proceedings of the tenth international symposium on information and communication technology_. 273–279.


# Inferring probabilistic Boolean networks from steady-state gene data samples

Vytenis Šliogeris, Leandros Maglaras, Sotiris Moschoyiannis Vytenis Šliogeris
and Sotiris Moschoyiannis are with the School of Computer Science & Electronic
Engineering, University of Surrey, GU2 7XH, UK (e-mail:
<EMAIL_ADDRESS>v.sliogeris@surrey.ac.uk). Sotiris Moschoyiannis
and Vytenis Sliogeris have been partly funded by UKRI Innovate UK, grant
77032. Leandros Maglaras is with the School of Computer Science and
Informatics, De Montfort University, Leicester, UK (e-mail:
leandros.maglaras@dmu.ac.uk).

###### Abstract

Probabilistic Boolean Networks have been proposed for estimating the behaviour
of dynamical systems as they combine rule-based modelling with uncertainty
principles. Inferring PBNs directly from gene data is challenging however,
especially when data is costly to collect and/or noisy, e.g., in the case of
gene expression profile data. In this paper, we present a reproducible method
for inferring PBNs directly from real gene expression data measurements taken
when the system was at a steady state. The steady-state dynamics of PBNs is of
special interest in the analysis of biological machinery. The proposed
approach does not rely on reconstructing the state evolution of the network,
which is computationally intractable for larger networks. We demonstrate the
method on samples of real gene expression profiling data from a well-known
study on metastatic melanoma. The pipeline is implemented using Python and we
make it publicly available.

_Keywords_ steady-state data samples, network structure, dynamics,
discretisation, predictor sets, perceptron, complex networks

##  1 Introduction

Rapid progress in the development of next-generation sequencing technologies
for genomics has provided valuable insights into complex biological systems
[1]. Modelling single-cell or gene networks is becoming increasingly
important. The question of modelling complex molecular regulatory networks is
an important one for bioinformatics. The goal of systems biology is to
intervene on the state of the cell, using the dynamics of the underlying
regulatory network. A model that could accurately represent such dynamics
could be used for analysis, including control [2, 3, 4, 5, 6], steady-state
distribution [7, 8, 9, 10], observability [11, 12, 13]. Such analyses aid the
development of genetic therapies [14].

Boolean Networks (BNs) were introduced for this purpose by Kauffman [15]. In
brief, a BN comprises a set of Boolean variables, each variable representing
the on/off state of a gene, while interactions between genes are expressed by
Boolean functions. It was found that even randomly generated BNs exhibit
behaviour reminiscent of gene regulatory networks, with naturally arising
attractor states which represent cell types or the phenotype [16, 17]. This
explains the popularity of BNs for modelling gene interactions [18, 19].

However, with few exceptions, gene expression data suggests a number of
possible successor states to any given state in a BN, thereby refuting the
determinism inherent in BNs. Thus, a probabilistic BN (PBN) was introduced by
Shmulevich et al. [20] in which the definition of a BN was adapted such that
for each gene, at each time point, a Boolean function (and predictor gene set)
is chosen with some conditional probability [21].

Inferring the PBN representation of a gene regulatory network (GRN) is quite
involved. First, the directed graph expressing interactions between genes
needs to be constructed; then, the Boolean functions need to be determined;
followed by determining the probabilities of selecting a Boolean function as
well as the number of candidate functions on each gene. Existing work (cf
Section 2) tends to focus on inference from time-series gene expression data
as the temporal aspect reveals the transition structure of the corresponding
PBN. However, as already pointed out in [22], there are concerns over the
number of (typically expensive to obtain) observations needed in such gene
microarray data. Approaches based on ODEs (e.g., [23]) require lots of
observations to tune the large number of parameters of the model, while in
practice only a handful are available. More such observations are available
when the underlying gene network is at a steady state [7], e.g, see gene
expression profiles of melanoma by Bittner et al [24].

In this paper, we propose a systematic method for inferring PBNs directly from
real gene expression data measurements, collected using microarray technology,
when the system is at a steady-state. The steady-state (long-run) behaviour of
a PBN is of interest to system biology as it allows to determine the long-term
influence of a gene on another gene or determine the long-term joint
probabilistic behaviour of a few selected genes [7].

The key contribution of our paper is a reproducible pipeline for going from
gene (steady-state) data samples to the PBN representation of the long-run
behaviour of the underlying genetic network. We use a predictor gene set
rather than temporal data to infer the ”transition structure”. Unlike other
proposals, our method does not require the construction of the probability
transition matrix, whose size grows exponentially on the number of nodes, and
hence becomes computationally intractable for larger networks [25].

The remainder of the paper is structured as follows. Section 2 outlines
related work. Preliminary background knowledge is presented in Section 3. The
main algorithm for our inference method is in Section 4. PBNs are produced in
Section 6 using the process described in Section 5. Concluding remarks are in
Section 7.

##  2 Related Work

There have been various methods for PBN inference, focusing on causality,
using different types of gene data [26]. Previous work on PBN inference from
time series gene data includes [27], SCODE [23] with ODEs, and most recently
the Stochastic Conjunctive Normal Form (SCNF) -based method by Apostolopoulou
et al [28] which can address larger networks.

Previous work on inference from steady-state data samples is relatively
limited and goes back to Shmulevich et al [7]. A tool for computing the ssd
probabilities has been proposed in [29]. Melkman et al [30] infer threshold
PBNs, a particular version of PBNs where every input threshold function of a
node must have the same number of parameters and also satisfy certain
stringent conditions. Kobayashi et al [10] construct PBNs from BNs by casting
inference as an integer linear programming problem and constructs a PBN that
fits the given steady-state distribution.

Kim et al. [31] use steady-state gene data samples from the study on
metastatic melanoma by Bittner et al. [24] (we use the same data here). They
choose the genes for their PBN using a combination of Coefficient of
Determination (COD) analysis and biological background knowledge (we do not
assume any prior knowledge). For the functions, they ternarise their data, and
construct Lookup Tables in place of the functions for each gene. They also
analyse the PBNs produced by analysing the steady-state distribution (ssd) of
the resulting network.

Shmulevich et al. [20], who introduced PBNs, describe a method for determining
functions for nodes in a PBN. This requires finding sets of input genes which
have high COD with the target gene, and using the predictive model used for
the calculation of the COD as the function for the particular set of input
genes. The probability for choosing the particular input gene set is
proportional to the COD of the input gene set.

Since discretisation of gene data is an important factor for inference. Chen
et al. [32] describe a method for quantising gene data using the expressions
of housekeeping genes within the dataset. Housekeeping genes are genes which
keep a constant expression, as they perform important functions within the
cell. Since they have a constant expression, they can be used to estimate the
probability distribution function (PDF) of the gene expressions within a
microarray. The constructed PDF can be used for using a hypothesis test to
determine whether or not a gene is over- or under-expressed. However, this
method hinges on knowledge of which of the genes are housekeeping genes and
this typically is not readily available.

As discussed in the introductory section, we focus on constructing PBNs from
real, microarray gene data samples, collected while the system is in a
sterady-state, instead of simulated, time-series data or starting from BNs. We
present a reproducible method to perform such a task.

##  3 Preliminaries

###  3.1 Boolean Networks

A BN [15] is a directed graph, $G=\\{V,E\\}$, comprised of vertices $V$ and
edges $E$. The vertices $v\in V$ represent the Boolean variables, which in
this case represent genes in a gene regulatory network. The directed edges
$\\{v_{i},v_{j}\\}=e_{i,j}\in E$ represent that one variable, $v_{i}$,
influences another, $v_{j}$. Each vertex is associated with a Boolean function
$f_{i}$ given by $f_{i}:\\{0,1\\}^{n_{in}}\mapsto\\{0,1\\}$. The input for
$f_{i}$ is a Boolean vector of length $n_{in}$, which represents the states of
all of the input vertices, and the output is a single Boolean value, which is
then used as the next state of the variable $v_{i}$. For a vertex $i$, the
input vertices are the vertices from which all incoming edges originate, given
by $\\{v_{j}|\exists\\{v_{j},v_{i}\\}\\}=e_{j,i}\in E$.

###  3.2 Probabilistic Boolean Networks

Probabilistic Boolean networks are an extension of Boolean networks. They are
directed graphs $G$, as in Boolean networks, except each function $f_{i}$ for
each node $i$ in the case of Boolean networks is replaced by a set of Boolean
functions $F_{i}=\\{f_{i}^{1},f_{i}^{2},\dots,f_{i}^{l_{i}}\\}$, and
probabilities $c_{i}=\\{c_{i}^{1},c_{i}^{2},\dots,c_{i}^{l_{i}}\\}$. Hence,
the logical function $f_{i}$ has $l_{i}$ possibilities, each with a
corresponding conditional probability of being selected at every time step.

More formally, during run time, a function $f_{i}^{j}$ for the node $v_{i}$ is
chosen with probability $c_{i}^{j}$, $j\in[1,l_{i}]$. PBNs are an extension to
BNs in the sense that if each node within a PBN has a single function, it
becomes identical to the BN.

###  3.3 State Transition Graphs

For each PBN there exists a state transition graph (STG). An STG is a directed
graph $G=\\{V,E\\}$, where the vertices $v_{i}\in V$ represent the possible
states of the PBN, and the edges $\\{v_{i},v_{j}\\}=e_{i,j}\in E$ represent
the possibility of a transition from state $v_{i}$ to $v_{j}$. Since the
probability of getting to another state $v_{j}$ only depends on the current
state $v_{i}$, we can say that the STG is a Markov chain.

By saying that the PBN has a steady state distribution (SSD), we mean that the
STG of the PBN has a steady state distribution. For an STG to have an SSD, it
needs to be ergodic \- that is, every state can reached from every other
state. To guarantee that the STG is ergodic, random perturbations with low
probability are introduced to the PBN.

###  3.4 Microarray Gene Data Samples

The data used to infer a PBN in our work was taken from the study of
metastatic melanoma found in Bittner et al [24], which has been extensively
studied in the literature [31, 33, 34, 6]. The study extracts and analyses the
gene expression profiles of 31 melanoma cells using microarray technology. To
make sure that the gene expression levels used in inferring the corresponding
PBN are those of genes when the network is in a steady state, the Kolmogorov-
Smirnov (KS) statistic is applied, as discussed in more detail in Section 5.

To utilise a particular gene in DNA, see [32], assuming the cell is at a
steady-state, the relevant segment of the molecule must first be transcribed,
producing messenger RNA (mRNA) which is accessible to the rest of the
proteins. The quantity of mRNA in a cell signifies the degree of protein
production associated with a particular gene.

DNA microarrays measure the presence of mRNA within a cell. The microarrays
consist of a surface with an array of robotically placed complementary DNA for
the genes to be analysed. mRNA tightly bonds with complementary DNA, hence the
microarray can be used to isolate different mRNA molecules. The process is
known as hybridisation.

The quantity of mRNA within a cell is measured by tagging the mRNA with
fluorescent molecules, hybridising them with a microarray, and exciting the
fluorescent molecules. The emitted brightness is proportional to the amount of
mRNA present.

Since the amount of mRNA differs depending on the gene, the data is normalised
by dividing the values recorded by the values recorded from a reference probe.
Since values recorded are non-negative, the ratio values are in the range of
$[0,\infty)$. Furthermore, since we would expect the values of within the
reference probe and the sample to not be different, the median for the ratio
values is expected to be 1. These are the values provided by Bittner et al.
[24] in the form of a matrix of size 8,150 (number of genes) by 31 (number of
samples). A small sample of the raw data is shown later in Fig. 1(a).

For demonstrating our method of inferring a PBN, we work with the subset of
melanoma genes analysed by Datta et al. [35], which are extensively studied in
the literature [31, 33, 10, 34, 6], namely WNT5A, pirin, S100P, RET1, MART1,
HADHB and STC2. This offers straightforward validation for our approach since
it produces the same PBN.

It is worth noting that larger PBNs may be constructed following the pipeline
described in this paper, and we have constructed the 28 node PBN given in [34]
as well as 70 node PBN which includes the 28 nodes already studied in [34]
padded with the 42 nodes with the highest weighting of importance, using
discriminative weights, which determine how a gene changes during the
experiment compared to the control cells [24]].

###  3.5 Coefficient of Determination

Coefficients of Determination (CODs) were described by Kim et al. [36] as a
method to determine which gene determines the state of which other gene. A COD
of a target variable, $Y$, with regards to an input variable, $X$, is a
measure on how well the target variable can be predicted using the input
variable. A predictive model $f$ is used to predict the value of the target
variable with and without the input variable, and compute the errors $\bar{e}$
and $e$ respectively. The relative change of error of the predictive model is
the COD $\theta$, given by eq. 1:

$\theta=\frac{\bar{e}-e}{\bar{e}}$ (1)

There are no constraints on what can be used as a predictive model. We opted
for a perceptron. This is because there exists a closed-form solution for
linear regression of the perceptron, described by Kim et al. [36], which can
be used instead of training. This aids in lowering the computation time.

The weights of a perceptron, $A$, can be computed using the closed form
solution:

$\begin{split}A&=R^{+}\cdot C\\\ R&=X\cdot X^{T}\\\ C&=X\cdot Y\\\
\end{split}$ (2)

###  3.6 Discretisation

Since PBNs use discrete values, the gene data which consists of real values
has to be discretised. Discretisation is a process where values get mapped
from the real value domain to the integer domain. For the problem at hand,
since genes can be in one of two states, the range of the function should be
either 0 or 1. Hence the function should take the form of:

$f:G\to G_{d},x\geq 0,\forall x\in G,y\in\\{0,1\\}\forall y\in G_{d}$ (3)

Such a method is described in detail in [37]. It consists of deciding upon a
threshold value $t$ with which all real values are compared. Each value then
gets mapped to 0 if it is below the threshold, and to 1 otherwise, as given by
eq. 4.

$G_{d}[x,y]=\left\\{\begin{array}[]{ll}0&G[x,y]<t\\\ 1&G[x,y]\geq t\\\
\end{array}\right.$ (4)

The threshold may be any metric. Common metrics are means or medians. The
threshold may also be the boundary between the top $x\%$ of entries and the
rest.

Shmulevich et al. [21] describe a process of using k-means clustering to
cluster the data, and assigning values to the data points depending on the
cluster they belong to. However, since half of the data points lie in the
range $(0,1)$, and the other half is in the range $(1,\infty)$, the lower
cluster ends up larger, resulting in a larger threshold that produces more
zeros. This can be remedied by performing k-means clustering on the logarithms
of the data points. This makes the ranges of both halves the same, producing
more representative clusters.

##  4 Inference of PBNs

In this section we describe the inference method and how it can be
implemented. Our approach to inferring a PBN starts with the real gene
expression data in the form of a matrix $G$ as input (see Fig. 1(a)), and
produces a PBN (see Fig. 1(b)). The input matrix is of size $m\times n$, where
$m$ is the number of genes and $n$ is the number of samples.

The method we apply for inferring PBNs draws upon work done by Shmulevich et
al. [20]. First, it requires the dataset to be discretised (recall Section
3.6). This process is performed in Algorithm 1.

Data: Input array G, discretisation axis, calculation method

Result: Discretised array G

for  _row g in G over axis_ do

t = method(G);

for  _element in g_ do

if  _element $<$ t_ then

element = 0;

else

element = 1;

end if

end for

return G ;

end for

Algorithm 1 Discretisation algorithm

Given a target gene, $n_{p}$ sets of genes with the highest CODs are found.
This is done following Algorithm 2.

Data: Input array G, number of input genes per tuple k, number of input tuples
n

Result: List of predictor-COD-input tuples

tuples = [];

for  _targetIndex in $1\dots$ G.geneAxisSize_ do

indexCombinations = generateAllCombinations($1\dots G.geneAxisSize\setminus
targetIndex$, k);

// Would return triplets of all possible k-combinations of the gene indexes
within the array.

buffer.init(n);

for  _inputCombo in indexCombinations_ do

COD, weights = calcCOD(G[targetIndex],G[inputCombo]);

if  _COD $>$ min(buffer.COD)_ then

buffer.add(COD, weights, inputCombo);

buffer.removeSmallestCOD() ;

end if

end for

tuples[targetIndex] $\leftarrow$ buffer;

end for

return buffer;

Algorithm 2 Tuple list generator

A buffer of size $n_{p}$ is initialised, and each possible combination of
input genes have their CODs calculated. If a combination of inputs has a COD
higher than at least one saved in the buffer, the buffer entry with the lowest
COD gets replaced by the new combination of inputs. This results in a buffer
full of input combinations with the highest CODs. One such buffer is
initialised per target gene, resulting in $n_{p}$ input combinations per
target gene.

During run-time, a set of input genes is chosen with probability proportional
to the COD of the set, and the next state is governed by the state of those
input genes in conjunction with the predictive model that was saved. For all
intents and purposes, the list with input gene, perceptron weights and
probabilities are enough to construct a PBN, as the input genes convey the
connectivity, and the perceptron weights convey the logic for that set of
input genes. The process is summarised in Algorithm 3.

Data: Input array G, discretisation axis, threshold calculation method, number
of input genes per tuple k, number of input tuples n

Result: List of predictor-probability-input tuples

$\hat{G}$ = discretise(G, axis, method);

predictorList = genPredictorList($\hat{G}$, k, n);

CODsum = 0 ;

for  _tuple in predictor list_ do

CODsum $\leftarrow$ CODsum + tuple[COD] ;

end for

for  _tuple in predictor list_ do

tuple[probability] $\leftarrow\frac{tuple[COD]}{CODsum}$ ;

end for

return predictor list;

Algorithm 3 General algorithm for generating predictor lists Figure 1: Input
and output for the inference method.

(a) The input, gene expression values for 31 samples taken at a steady-state,
given in the form of a matrix

(b) The output is a PBN

##  5 Analysis

The analysis of the generated PBNs in our approach are based on steady-state
distribution, which is fairly standard, e.g., see [31]. The PBN is run for $T$
steps in order to get it within a steady state. Then it is run for the next
$N$ steps, recording the state it is at. To confirm whether or not the PBN is
in a steady state after $T$ steps, the Kolmogorov-Smirnov (KS) statistic is
calculated for the two halves of $N$.

The entries recorded in $N$ are split in to two halves - one containing states
$[0,\frac{N}{2}]$, the other containing $[\frac{N}{2}+1,N]$. The entries are
subsampled with the interval $G$. The histograms are converted to cumulative
distribution functions, and the maximum vertical distance between them is
found, which is the KS statistic.

The significance test shows the probability of the two CDFs being drawn from
the same distribution. If the PBN had not reached a steady state after $T$
steps, the halves of $N$ would be drawn from different distributions, which
would be indicated by the KS test. The recorded states are a string of binary
values. Therefore, for ease of analysis, they are used as gray-coded integers,
and displayed on a histogram (cf. Fig. 2). This makes the horizontal distance
on the histogram proportional to the Hamming distance between two network
states.

##  6 Evaluation

We have implemented the pipeline using Python $3$ and made it publicly
available on https://github.com/UoS-PLCCN/pbn-inference.

We have constructed PBNs of size 7 from data produced by Bittner et al. [24]
using different thresholds for the quantisation methods. The thresholds were
(a) average of a gene expression; (b) median of a gene expression, and (c)
k-means clustering of a gene expression. The data was quantised on a per-gene
basis, with each gene having 10 triplets of input genes.

For the construction and validation of the histograms representing the steady-
state distribution, we have chosen the parameters to be $T=10^{6}$, $N=4\cdot
10^{6}$, $G=10$ and $R=100$. On a laptop with 32 GB of RAM and an Intel® Core™
i7-7700HQ Processor, each histogram took around 9 hours to produce. The
results are shown in Fig. 2.

Figure 2: SSDs of PBNs generated using different quantisation methods.
States on the x-axis; SSD probability on the y-axis

(a) Mean

(b) Median

(c) K-Means clustering

It can be seen that the average and the median quantisation methods produce
very similar histograms, with three peaks each, and the latter two peaks being
in similar positions. The histogram generated using the PBN constructed from
k-means clustering only has one prominent state, which can also be observed in
the other two PBNs. It may be constructive to note that the few very prominent
states in the histograms show in Fig. 2 agrees with the assumption claimed by
Kim et al. [31] that gene regulatory networks found in nature only occupy a
small fraction of the possible state space.

For the purposes of direct comparison, we have trialled the proposed method in
the DREAM (Dialogue on Reverse Engineering Assessment and Methods)
challenge111http://dreamchallenges.org/project/dream-3-in-silico-network-
challenge/ which offers a benchmark for network inference (DREAM 3) [38] and
scored 8th (out of 29).

##  7 Conclusion

In this work we described the inference a PBN directly from real gene data,
collected using microrarray technology, which were taken when the system was
at a steady-state. This kind of gene profiling is typically less costly to
obtain than time series data, and includes more data points. Using the
evaluation methods described in the literature, e.g., by Kim et al. [31], we
have concluded that the pipeline works well for the examples provided.
However, it is subject to fine-tuning the parameters. We have provided the
method in a systematic pipeline which can be reproduced and made it publicly
available on github https://github.com/UoS-PLCCN/pbn-inference.

We note that the method scored 8th (out of 29) in the DREAM challenge and has
been used to infer large PBNs (N=200).

It is worth noting is that the proposed method does not require a state
transition probability matrix to be produced. It can be extracted from the
PBN, however, the time required grows exponentially with the size of the PBN.
This means that conventional mathematical methods in the literature that make
use of the transition probability matrix may not always be applicable.

One concern is that the transitions get fitted to the quantised dataset. It is
widely accepted that the states observed in the dataset are steady states of
the cells. Since the transition rules get fitted to the steady states of the
cells, the resulting PBN will be driven towards the steady states observed
within the data. However, while it is certain that the method captures the
long-run behaviour (steady-state) of the underlying gene regulatory network,
there is little certainty that the PBN will behave with biological accuracy
between the observed steady states. This concern could possibly be addressed
by using time-series gene data to augment the method presented here, as this
type of data captures the change of gene expression levels with respect to
time. This promises to capture the behaviour at and between steady states,
without reconstruction of the state evolution of the PBN, and is certainly
worth exploring further in future work.

## References

  * [1] C. Gawad, W. Koh, and S. Quake.  Single-cell genome sequencing: current state of the science.  Nat Rev Genet, 17:175–188, 2016.
  * [2] Y-Y. Liu, J-J. Slotine, and Albert-László Barabási.  Controllability of complex networks.  Nature, 473(7346):167, 2011.
  * [3] M. R. Karlsen and S. Moschoyiannis.  Evolution of control with learning classifier systems.  Applied Network Science, 3(1):1–30, 2018.
  * [4] G. Papagiannis and S. Moschoyiannis.  Learning to control random Boolean networks: A deep reinforcement learning approach.  In Complex Networks 2019, volume 881, pages 721–734. Springer, Cham, 2019.
  * [5] Y. Wu and T. Shen.  Policy iteration algorithm for optimal control of stochastic logical dynamical systems.  IEEE Transactions on Neural Networks and Learning Systems, 29(5):2031–2036, 2019.
  * [6] Georgios Papagiannis and Sotiris Moschoyiannis.  Deep reinforcement learning for control of probabilistic boolean networks.  In Complex Networks 2020, volume 944, pages 361–371. Springer, 2020\.
  * [7] I. Shmulevich et al.  Steady-state analysis of genetic regulatory networks modelled by probabilistic Boolean networks.  Comp Funct Genomics, 4(6):601–608, 2003.
  * [8] S. Moschoyiannis and M.W. Shields.  A set-theoretic framework for component composition.  Fundamenta Informaticae, 59(4):373–396, 2004.
  * [9] W.-K. Ching, M. K. Ng Zhang, and T. Akutsu.  An approximation method for solving the steady-state probability distribution of probabilistic Boolean networks.  Bioinformatics, 23(12):1511–1518, 2007.
  * [10] K. Kobayashi and K. Hiraishi.  Design of probabilistic Boolean networks based on network structure and steady-state probabilities.  IEEE Transactions on Neural Networks and Learning Systems, 28(8):1966–1971, 2017.
  * [11] Qunxi Zhu, Yang Liu, Jianquan Lu, and Jinde Cao.  Controllability and observability of boolean control networks via sampled-data control.  IEEE Trans. Control. Netw. Syst., 6(4):1291–1301, 2019.
  * [12] Kuize Zhang and Karl Henrik Johansson.  Efficient verification of observability and reconstructibility for large boolean control networks with special structures.  IEEE Transactions on Automatic Control, 65(12):5144–5158, 2020\.
  * [13] S. Savvopoulos and S. Moschoyiannis.  Impact of removing nodes on the controllability of complex networks.  In Complex Networks 2017, 2017.
  * [14] H. F. Fumia and M. L. Martins.  Boolean network model for cancer pathways: Predicting carcinogenesis and targeted therapy outcomes.  PLoS ONE, 8(7):e69008, 2013.
  * [15] S. A. Kauffman.  Metabolic stability and epigenesis in randomly constructed genetic nets.  Journal of Theoretical Biology, 22(3):437 – 467, 1969.
  * [16] Dimitrios Voukantsis, Kenneth Kahn, Martin Hadley, Rowan Wilson, and Francesca M Buffa.  Modeling genotypes in their microenvironment to predict single- and multi-cellular behavior.  GigaScience, 8(3), 01 2019.
  * [17] Evangelos Chatzaroulas, Vytenis Sliogeris, Pedro Victori, Francesca M. Buffa, Sotiris Moschoyiannis, and Roman Bauer.  A structural characterisation of the mitogen-activated protein kinase network in cancer.  Symmetry, 14(5), 2022.
  * [18] R. Albert and H. G. Othmer.  The topology of the regulatory interactions predicts the expression pattern of the segment polarity genes in Drosophila melanogaster.  Journal of Theoretical Biology, 223(1):1–18, 2003.
  * [19] M. Davidich and S. Bornholdt.  The transition from differential equations to Boolean networks: a case study in simplifying a regulatory network model.  J Theor Biol, 255(3):269–77, 2008.
  * [20] I. Shmulevich, E. R. Dougherty, S. Kim, and W. Zhang.  Probabilistic Boolean networks: a rule-based uncertainty model for gene regulatory networks.  Bioinformatics, 18 2:261–74, 2002.
  * [21] I. Shmulevich and E. R. Dougherty.  Probabilistic Boolean Networks: The Modeling and Control of Gene Regulatory Networks.  SIAM, 2010.
  * [22] Z. Bar-Joseph.  Analyzing time series gene expression data.  Bioinformatics, 20(16):2493–2503, 2004.
  * [23] H. Matsumoto et al.  SCODE: an efficient regulatory network inference algorithm from single-cell RNA-Seq during differentiation.  Bioinformatics, 33(15):2314–2321, 2017.
  * [24] M. Bittner et al.  Molecular classification of cutaneous malignant melanoma by gene expression profiling.  Nature, 406:536–40, 09 2000.
  * [25] Tatsuya Akutsu et al.  Control of boolean networks: Hardness results and algorithms for tree structured networks.  J. of Theoretical Biology, 244(4):670–679, 2007.
  * [26] C. Glymour, K. Zhang, and P. Spirtes.  Review of causal discovery methods based on graphical models.  Frontiers in Genetics, 10(524), 2019.
  * [27] A. Silescu and V. Honavar.  Temporal boolean network models of genetic networks and their inference from gene expression time series.  Complex Systems, 13(2001):61–78, 2001.
  * [28] I. Apostolopoulou and D. Marculescu.  Tractable learning and inference for large-scale probabilistic Boolean networks.  IEEE Transactions on Neural Networks and Learning Systems, 30(9), 2019.
  * [29] Andrzej Mizera, Jun Pang, and Qixia Yuan.  Assa-pbn: An approximate steady-state analyser of probabilistic boolean networks.  In Automated Technology for Verification and Analysis, pages 214–220, Cham, 2015. Springer International Publishing.
  * [30] Avraham A. Melkman, Xiaoqing Cheng, Wai-Ki Ching, and Tatsuya Akutsu.  Identifying a probabilistic boolean threshold network from samples.  IEEE Transactions on Neural Networks and Learning Systems, 29(4):869–881, 2018.
  * [31] S. Kim, E. R. Dougherty, Y. Chen, M. Bittner, and E. Suh.  Can markov chain models mimic biological regulation?  Journal of Biological Systems, 10, 03 2003.
  * [32] Y. Chen.  Ratio-based decisions and the quantitative analysis of cDNA microarray images.  Journal of Biomedical Optics, 2(4):364, 1997.
  * [33] R. Pal, A. Datta, and E. R. Dougherty.  Optimal infinite-horizon control for probabilistic Boolean networks.  IEEE Transactions on Signal Processing, 54(6):2375–2387, 2006.
  * [34] U. Sirin, F. Polat, and R. Alhajj.  Employing batch reinforcement learning to control gene regulation without explicitly constructing gene regulatory networks.  In International Joint Conference on Artificial Intelligence (IJCAI), pages 2042–2048, 08 2013.
  * [35] A. Datta, A. Choudhary, M. L. Bittner, and E.R. Dougherty.  External control in markovian genetic regulatory networks.  Machine Learning, 52:3614 – 3619 vol.4, 07 2003.
  * [36] S. Kim, E. Dougherty, M. Bittner, Yidong Chen, Krishnamoorthy Sivakumar, Paul Meltzer, and Jeffrey Trent.  General nonlinear framework for the analysis of gene interaction via multivariate expression arrays.  Journal of Biomedical Optics, 5:411–24, 11 2000.
  * [37] C. Velarde, C. Rubio-Escudero, and R. Romero-Zaliz.  Boolean networks: a study on microarray data discretization.  ESTYLF08, Cuencas Mineras (Mieres-Langreo), pages 17–19, 2008.
  * [38] Daniel Marbach, Robert J. Prill, Thomas Schaffter, Claudio Mattiussi, Dario Floreano, and Gustavo Stolovitzky.  Revealing strengths and weaknesses of methods for gene network inference.  Proceedings of the National Academy of Sciences, 107(14):6286–6291, 2010.