# Learning to search for and detect objects in foveal images using deep
learning

Beatriz Paula  Instituto Superior Técnico
Universidade de Lisboa
Portugal
<EMAIL_ADDRESS>Plinio Moreno  Insitituto de Sistemas e
Robótica & Instituto Superior Técnico
Universidade de Lisboa
Lisboa, Portugal
<EMAIL_ADDRESS>
###### Abstract

The human visual system processes images with varied degrees of resolution,
with the fovea, a small portion of the retina, capturing the highest acuity
region, which gradually declines toward the field of view’s periphery.
However, the majority of existing object localization methods rely on images
acquired by image sensors with space-invariant resolution, ignoring biological
attention mechanisms.

As a region of interest pooling, this study employs a fixation prediction
model that emulates human objective-guided attention of searching for a given
class in an image. The foveated pictures at each fixation point are then
classified to determine whether the target is present or absent in the scene.
Throughout this two-stage pipeline method, we investigate the varying results
obtained by utilizing high-level or panoptic features and provide a ground-
truth label function for fixation sequences that is smoother, considering in a
better way the spatial structure of the problem.

Finally, we present a novel dual task model capable of performing fixation
prediction and detection simultaneously, allowing knowledge transfer between
the two tasks. We conclude that, due to the complementary nature of both
tasks, the training process benefited from the sharing of knowledge, resulting
in an improvement in performance when compared to the previous approach’s
baseline scores.

###### Index Terms:

Visual Search, Object Detection, Scanpath Prediction, Foveal Vision, Deep
Learning

##  I Introduction

Computer Vision (CV) is an interdisciplinary area that connects several lines
of study such as Computer Science, Biology, Psychology, Engineering,
Mathematics and Physics. From an engineering standpoint, it attempts to
understand and model the visual system in order to interpret the visual world
[1]. It has many applications such as object detection, face recognition and
action and activity recognition, in a multitude of fields from military to
medicine.

A fundamental difference between the human visual system and current
approaches to object search is the acuity of the image being processed [2].
The human eye captures an image with a very high resolution in the fovea, a
small region of the retina, and a decrease in sampling resolution towards the
periphery of the field of view. This biological mechanism is crucial for the
real-time image processing of the rich visual data that reaches the eyes
(0.1-1 Gbits), since visual attention prioritizes interesting and visually
distinctive areas of the scene, known as salient regions, and directs the gaze
of the eyes. In contrast, image sensors, by default, are designed to capture
the world with equiresolution in a homogeneous space invariant lattice [3],
and current solutions to vision system performance rely on the increase of the
number of pixels. This limits real-time applications due to the processing
bottleneck and the excessive amount of energy needed by state-of-the-art
technologies.

In the last decade, Deep Learning (DL) has shown tremendous success when
compared to traditional Machine Learning (ML) approaches. First proposed in
the 80s [4], the Convolutional Neural Network (CNN), a DL technique, is
inspired by the human visual processing system. In 2012, in the ImageNet Large
Scale Visual Recognition Challenge, the winner, Alex Krizhevsky, introduced a
CNN implementation solution [5] that showed its massive power as a training
architecture. The AlexNet shows a lot of similarities with LeNet [6] that was
published in 1998. However, by scaling up both the data and the computational
power, we are only now seeing the true potential of Deep Learning in Computer
Vision.

Nevertheless, it remains challenging to replicate and model the human visual
system. In an attempt to solve this, there have been advancements in combining
these DL techniques with image foveation and saliency detection models. It has
been shown that a ”foveated object detector can approximate the performance of
the object detector with homogeneous high spatial resolution processing while
bringing significant computational cost savings” [7].

The primary objective of this thesis is to utilize goal-guided scanpath data
for object detection in images with foveated context. The following is a
description of the problem we intend to solve: given an input image and an
object category, indicate the presence or absence of instances of that class
in the scene while adjusting the acuity resolution to mimic the human visual
system.

Our contributions include: (i) Benchmark of recent approaches based on Deep
Learning, which are able to predict fixations, on a recent large-scale
dataset; (ii) a ground-truth label function for fixation sequences that is
smoother, considering in a better way the spatial structure of the problem;
(iii) evaluation of two alternative visual representations (conventional high-
level features from VGG and a more elaborate multi-class presence
description); and (iv) the introduction of a novel dual task approach that
simultaneously performs fixation and target detection.

##  II Related Work

Human attention is driven by two major factors, bottom-up and top-down factors
[8]. While bottom-up is driven by stimulus properties such as distinct
features in the field of vision, which means saliency detection is executed
during the pre-attentive stage, top-down factors are influenced by
behaviourally relevant stimuli, such as prior knowledge, expectations and
goals [9]. This latter factor can be detected in the 1967 experiment where
individuals where asked to analyze a scene with different goals [10]. In this
experiment individuals observed a picture of a family in a room with a visit,
and were to answer different questions, such as to freely examine, to
determine the people’s ages, the material circumstances of the family, etc.
According to Yarbus “depending on the task in which a person is engaged, i.e.,
depending on the character of the information which he must obtain, the
distribution of the points of fixation on an object will vary correspondingly,
because different items of information are usually localized in different
parts of an object”.

The CV task of Gaze Prediction aims to predict fixation patterns made by
people in image viewing and can have a spatial representation, in fixation
density maps, and an added temporal representation when predicting scanpaths.
In this area of study, most work focuses on free-viewing, which, as mentioned,
is led by bottom-up attention.

In the method proposed by Ngo and Manjunath [11], CNNs are used for feature
extraction and feature maps compilation which are then used in a Long Short
Term Network (LSTM) responsible for modeling gaze sequences during free-
viewing. This latter network was introduced in [12] as a solution to the
vanishing gradient problem of RNNs that prevented the understanding of long
dependencies since the correlation structure died down during the backward
pass. These networks have a more complex structure which tweak the hidden
states with an additive interaction, instead of a linear transformation
previous recurrent networks performed, allowing the gradient to fully
backpropagate all the way to the first iteration.

However, as previously mentioned, human scanpaths during search tasks vary
depending on the target items they are trying to gain information from,
therefore guided search cannot be predicted based on free-viewing knowledge,
where there were no explicit goals. Goal-directed attention is additionally
relevant due to the human search efficiency in complex scenes that accounts
for scene context and target spatial relations [13].

In [14], a similar approach to the free-viewing scanpath predictor presented
in [11] was taken. However, they leveraged a Convolutional Long Short Term
Memory (ConvLSTM) architecture, and introduced a foveated context to the input
images on top of an additional input encoding the search task, which found
human fixation sequences to be a good foundation for object localization. The
ConvLSTM had been previously introduced in [15] as a variant of LSTMs better
suited for 3-dimensional inputs, such as images. This adaptation still
contains the same two states: a hidden state, $h$, and a hidden cell state,
$c$; and the same four intermediate gates: the input gate $i$, forget gate
$f$, output gate $o$ and candidate input $\tilde{c}$; as the LSTM
architecture. However, a convolution is performed during the computation of
the gates instead of the previous product operations, as seen in the following
equations:

$i_{t}=\sigma(W_{i}*x_{t}+U_{i}*h_{t-1})$ (1)
$f_{t}=\sigma(W_{f}*x_{t}+U_{f}*h_{t-1})$ (2)
$o_{t}=\sigma(W_{o}*x_{t}+U_{o}*h_{t-1})$ (3)
$\tilde{c}_{t}=\tanh{(W_{c}*x_{t}+U_{c}*h_{t-1})}$ (4) $c_{t}=f_{t}\odot
c_{t-1}+i_{t}\odot\tilde{c}_{t}$ (5) $h_{t}=o_{t}\odot tanh(c_{t})$ (6)

where $\odot$ denotes an element wise product, $*$ denotes a convolution, and
$W$ and $U$ are the weight matrices of each gate that operate over the hidden
states.

The limited amount of available data containing human scanpaths in visual
search was, however, identified as a significant obstacle in [14]. Since then
a new large-scale dataset has been introduced in [16], which has shown
promising results in [17], where an inverse reinforcement learning algorithm
was able to detect target objects by predicting both the action (fixation
point selection) and state representations at each time step, therefore
replicating the human attention transition state during scanpaths. This
approach additionally utilized features extracted from a Panoptic FPN model
[19], that performs panoptic segmentation which is the unification of ”the
typically distinct tasks of semantic segmentation (assign a class label to
each pixel) and instance segmentation (detect and segment each object
instance)” [18].

Lastly, recent works in CV have begun experimenting with a Transformer-styled
design, as it has been demonstrated that such architectures offer top
performance in the field of Natural Language Processing [21]. For example,
visual attention-driven transformers have been applied to the medical field to
aid in diagnostics, as seen in [23] and [24] where scanpath data provided
task-guided visual attention that aided in disease classification in
histopathology images and chest radiographs, respectively.

##  III Systems Overview

In this section, we present the architecture of the two strategies used in
this study: a two-stage pipeline system consisting of a gaze fixations
predictor and an image classifier, and a dual-task model that conducts
scanpath prediction and target detection simultaneously.

###  III-A Fixation Prediction Module

We undertake tests on high-level and panoptic image feature inputs within the
fixation prediction module. Despite configuring a unique architecture for each
model, they both share the same structure. At each time-step $T=t$, the Input
Transformation Section aggregates the features of the foveated pictures at
each fixation location since the beginning of the gaze sequence, $T\in
0,...,t$, as well as the task encoding of the target object. This combined
input is then sent to the Recurrent Section, which uses ConvLSTM layers to
emulate human-attention through its hidden states. The Recurrent Section then
outputs its final hidden state $h_{t+1}$ to the model’s Output Section, which
predicts the next scanpath fixation as a discrete location in an image grid
with dimensions $H\times W$. We will now present the architecture of each
model in detail.

####  III-A1 Fixation Prediction from High-Level Features

In this model, we utilized the high-level features retrieved from the
ImageNet-trained VGG16 model [25, 26] with dimensions $H\times W\times Ch$,
and it is composed by the following sections:

  * •

Input Transformation: To aggregate the feature maps and the task encoding, we
opted for a Multiply layer which executes the element-wise multiplication of
these inputs. In addition, depending on its format, the task encoding may be
transmitted through a Fully Connected (FC) Layer with $Ch$ units and a tanh
activation, followed by a Dropout Layer with a rate of $r_{Dropout}$ in order
to prevent overfitting.

  * •

Recurrent Section: This portion mainly consists on a ConvLSTM layer composed
of $F$ filters, which correspond to the dimensionality of its output, a kernel
size of $K$ x $K$, a stride of $S$, and a left and right padding of $P$. This
layer has a tanh activation, and the recurrent step utilizes a hard sigmoid
activation, a piece-wise linear approximation of the sigmoid function, for
faster computation. Subsequently, to prevent overfitting we perform batch
normalization, where we normalize the vector $h$ with the batch mean $\mu$ and
batch variance $\sigma$ of the current batch during training, as seen in
equation 7, where $\epsilon$ is a small constant, and $\alpha$ and $\beta$ are
learned scaling and offset factors, respectively. During inference, the arrays
are normalized with a moving mean $M$ and moving variance $V$, which are non-
trainable variables that were determined throughout training with equation 8,
where $\mu$ and $\sigma$ are the mean and variance of the test batch and
$\gamma$ is the momentum hyper-parameter.

$y=\alpha\cdot\frac{x-\mu}{\sqrt{\sigma+\epsilon}}+\beta$ (7)
$M=M\cdot\gamma+\mu\cdot(1-\gamma)\qquad V=V\cdot\gamma+\sigma\cdot(1-\gamma)$
(8)

  * •

Output Section: We first perform a flattening operation to each temporal slice
of the input with the help of a Time Distributed wrapper. The flattened array
is then fed to a FC layer and has $H\times W$ units and a softmax activation
function.

####  III-A2 Fixation Prediction from Panoptic Features

To compute these new features we resorted to the Panoptic FPN model,
introduced in [19]. Due to the higher inference time of this model, of 53.0
ms/image, when compared to the VGG16, of 4.2 ms/image, with a GPU, we pre-
computed the belief maps for each image at a high resolution, with the
original image at full acuity, and a low resolution, after applying a uniform
Gaussian blur filter with a blur radius of $\sigma$. This way we simulate the
human foveal system by computing the belief map $B$ as:

$B(t)=M_{t}\odot H+(1-M_{t})\odot L$ (9)

where $H$ and $L$ are the belief maps for the high and low resolution images,
respectively, and $M_{t}$ is the binary fixation mask at time step $t$, of
size $H\times W$, where every element is set to $0$ except the grid cells at
an euclidean distance shorter than $r$ from the current fixation point.

To duplicate the belief maps in [17], the task encoding is now a one-hot
encoding with dimensions $H\times W\times Cl$, where each row of the axis $Cl$
corresponds to an object class and the two-dimensional map $H\times W$ is all
set to one for the target class and zero for the others. This input is
subsequently transmitted to the Input Transformation stage, where it is
concatenated with the image feature maps.

The recurrent section of the model is composed of $d$ ConvLSTM layers, each
followed by a Batch Normalization layer. Every ConvLSTM is constructed with
the same hyper-parameters: each one has $F_{LSTM}$ filters with a kernel size
of $K_{LSTM}\times K_{LSTM}$, a stride of $S_{LSTM}$, a padding of $P_{LSTM}$,
a tanh activation function and a hard sigmoid activation during the recurrent
step.

Finally, in the output section, we conducted experiments over two different
setups. The first one is composed of a 3d-Convolutional layer with a sigmoid
activation followed by a time distributed flattening operation. The second
setup is comprised of the same 3d-Convolutional layer, but with a ReLu
activation, and a flattening layer followed by a Fully Connected layer with
softmax activation.

###  III-B Target Detection Module

In the last stage of the pipeline, the model finally gets to evaluate, at each
time-step of the sequence, if the fixation point coincides with the location
of the target object. To achieve this, we chose to once again leverage a state
of the art classifier, the VGG16, and develop 18 binary classifiers, one for
each task. In addition, as a baseline, we utilize a complete VGG16 trained on
the ImageNet dataset to perform classification on our data, replicating an
approach taken in [27].

In order to fine-tune the already pre-trained VGG16 model, we substituted its
classification layers, after the last convolutional layer, with three fully
connected layers, $FC_{i}$ with $i\in\\{1,2,3\\}$, each with $U_{i}$ units,
where $FC_{1}$ and $FC_{2}$ were followed by a ReLu activation function and
$FC_{3}$ was followed by a sigmoid action function. During training, only the
parameters of these last $FC$ layers were updated.

###  III-C Dual Task Model

Our multitasking approach to the problem of attention-guided object
localization considers three different architectural alternatives. The primary
distinction between them is the employment of a ConvLSTM layer in both task
branches in the first two models, but only in the fixation prediction branch
in the third. All models receive as input the high-level feature maps, with
dimensions $H\times W\times Ch$ and a one-hot task encoding array, of size
$Cl$, which are then aggregated. Similar to the subsection III-A1, the
grouping of both of these inputs is accomplished by passing the task encoding
through a fully connected layer with $Ch$ units and tanh activation, and
conducting an element-wise multiplication with the foveated image’s feature
maps. After this shared module, the models branch off to complete each
specific task using the following architectures:

  * •

Architecture A After performing the input transformation, where we aggregate
the feature maps and task encoding, the array $x_{t}$ is fed to two ConvLSTM
layers, each containing $F$ filters of size $K\times K$ with a stride of $S$,
and an output tanh activation and a hard sigmoid activation in the recurrent
step. However, following each iteration of the fixation prediction recurrent
module, its internal states $h_{t}^{fix}$ and $c_{t}^{fix}$ are passed to the
detection branch as the internal states, $h_{t-1}^{det}$ and $c_{t-1}^{det}$,
of the preceding time step, as illustrated in figure 1(a). Subsequently, the
compositions of the fixation prediction and target detection branches resemble
those described in sections III-A1 and III-B, respectively. In the first
branch, a temporal flattening operation is performed to $h^{fix}$, followed by
an output layer consisting of a FC layer with softmax activation. In the
second we classify each temporal slice of $h^{det}$ by employing the same
structure of three FC layers $FC_{i}\in\\{1,2,3\\}$, with $U_{i}$ units,
respectively, where the first two layers have a ReLu activation while the
output layer has a sigmoid activation.

  * •

Architecture B Similar to the previous architecture, each task branch employs
ConvLSTM layer. The sole difference is that we now conduct the iterations of
the detection module first, and send the internal states $h_{t}^{det}$ and
$c_{t}^{det}$ to the fixation prediction module for the preceding time steps
$t-1$

  * •

Architecture C This design contains a fixation prediction branch that is
identical to that of architecture A. Nevertheless, the target detection branch
no longer has a ConvLSTM layer. Instead, at each time step $t$, the combined
input $x_{t}$ computed by the shared module is concatenated with the output of
the ConvLSTM layer, $h_{t+1}^{fix}$ of the fixation prediction task, as
illustrated in figure 1(b). Finally, this concatenation is followed by the
same three FC layers utilized by the previous architectures.

(a) Architecture A.

(b) Architecture C.

Figure 1: Information flow across the fixation prediction and target detection
branch in architectures A, on the top, and C, on the bottom.

##  IV Implementation

In this section we present the dataset, in addition to the conditions in which
our models were trained and tested.

###  IV-A Dataset

The dataset utilized in this study was the COCO-Search18 dataset, introduced
in [16]. In terms of the amount of images, target classes, and fixations, it
remains the largest dataset on goal-directed behavior, to our knowledge. This
dataset consists of 6,202 images from the Microsoft COCO dataset [28], evenly
split between target-present and target-absent images, of 18 target categories
(bottle, bowl, car, chair, analogue clock, cup, fork, keyboard, knife, laptop,
microwave, mouse, oven, potted plant, sink, stop sign, toilet and tv), with
eye movement recordings from 10 individuals.

As humans were able to fixate the target object with their gaze within their
first six saccades 99% of the time, as illustrated in figure 2, fixation
sequences with length greater than that were discarded. Additionally, the
sequences were padded with a repeated value of the last fixation point to
achieve a fixed length of 7, including the initial center fixation. This was
done to replicate the procedure of a similar work [27], where participants
were instructed to fixate their gaze on the target object, once they found
them, during search tasks.

To train the fixation prediction module and the dual task model, we used a
random dataset split of 70% train, 10% validate and 20% test over each class
category and all images were resized to 320 x 512 which resulted in feature
maps with $10\times 16$ spacial dimensions. For the fine-tuning of the binary
classifiers in the target detection module, due to the cluttered nature of the
images in the dataset, we utilize a cropped region of the images surrounding
each fixation point, with its size proportional to the median dimensions
within each class.

Figure 2: Cumulative probability of target fixation throughout sequences of
the COCO-Search18 dataset train split

###  IV-B Training

During the training phase, all our models were optimized with the Adam
algorithm [31] and a learning rate of $lr=0.001$, for a maximum of 100 epoch
with an early stopping mechanism activated when the validation loss stops
improving after a duration of 5 epochs. Additionally, every dropout is
performed with $r_{Dropout}=0.5$, the batch normalization uses
$\epsilon=0.001$ and $\gamma=0.99$, and we use a batch size of 256 in every
module apart from the fixation prediction performed with high-level features.

  * •

Fixation Prediction from High-Level Features During training, we estimate the
best weight and bias parameters that minimize the loss between the predicted
output $\hat{y}$ and the ground truth label $y$, with the cross entropy
function computed for every fixation time step $t$ for every sequence $s$ of
each mini-batch $b$:

$L_{CE}=-\sum_{s=1}^{S}\sum_{t=0}^{T}\sum_{i=1}^{H\times
W}y_{i}*log(\hat{y}_{i})$ (10)

where $S$ corresponds to the batch size, $T$ corresponds to the sequence
length which is set to 6 (in addition to the initial fixation point at $t=0$)
and $H\times W$ is the output size which is set to 160.

We set with $F=5$ filters, a kernel size of $K=4$ and a stride of $S=2$, and,
during training, we varied the batch size between 32, 64, 128 and 256 and we
conducted an ablation study over theses additional hyper-parameters and
settings:

    * –

Fovea size: In this work we utilized the same real-time foveation system as in
[14], and assessed the model performance when varying the fovea size, which
defines the radius of the region with highest visual acuity with the values of
50, 75 and 100 pixels.

    * –

Task encoding: In order to encode the object class being searched, we
experimented with several representations. The first was a one-hot encoding
array of size 18 representing each class. The second was a normalized heat map
of fixations made during the observations of that same task, compiled
exclusively with training data. This encoding was represented both with a two
dimensional array of size 10 x 16 and a one dimensional array of the flattened
grid with size 160. When utilizing a one dimensional encoding, this array was
passed through a FC layer with 512 units during the Input Transformation,
while the two dimensional encoding was fed directly into the Multiply layer.

    * –

Ground truth: We consider both a one-hot encoding representation of the
ground-truth label and a two dimensional Gaussian function with the mean set
to the cell coordinates of the actual fixation location and the variance set
to 1.

  * •

Fixation Prediction from Panoptic Features During the training of this setup
of the fixation prediction module we once again aim to minimize the loss
function in equation 10. Additionally, the feature maps used have dimension
$10\times 16\times 134$ in order to replicate the scale of our grid shaped
output.

In regards to the model architecture, we fixed the ConvLSTM layers’
configuration to $F_{LSTM}=10$ filters with square kernels of size
$K_{LSTM}=3$, a stride of $S_{LSTM}=1$ and a padding of $P_{LSTM}=1$ to
maintain the features spatial resolution. In the output section, we set the
$F_{Conv}=1$ filter of the 3d-Convolutional layer to have a kernel size of
$K_{Conv}=2$, a stride of $S_{Conv}=1$ and padding of $P_{Conv}=1$. The second
setup of this section is configured to have a Fully Connected layer with 160
units.

For this model, we additionally varied the depth of the recurrent section with
$d\in\\{1,3,5\\}$, and altered the structure of the output section to utilize
both a sigmoid and softmax as the final activation function. Concerning the
data representation, we once again evaluated the impact of having a one-hot or
a Gaussian ground truth encoding, and explored several belief maps settings:
we varied the radius $r$ of the mask, $M_{t}$, with values $r\in\\{1,2,3\\}$
(each to emulate a corresponding fovea size of 50, 75 and 100); and
experimented with a cumulative mask configuration, $M^{\prime}_{t}$, where the
binary mask utilized in equation 9, in addition to the information of the
current time step, accumulates the high acuity knowledge of all previous time
steps. All panoptic feature maps were computed with a low resolution map $L$
extracted from a blurred input image with a Gaussian filter with radius
$\sigma=2$.

  * •

Target Detection The binary classifiers were implemented with each fully
connected layer having $U_{1}=512$, $U_{2}=256$ and $U_{3}=1$ units, a dropout
rate of $r_{Dropout}=0.5$, and we varied the fovea size between 50, 75 and 100
pixels. They were trained with a loss function defined as:

$L_{BCE}=-\frac{1}{N}\sum_{i=1}^{N}y_{i}\cdot log(\hat{y}_{i})+(1-y_{i})\cdot
log(1-\hat{y}_{i})$ (11)

  * •

Dual Task Concerning the architecture, the ConvLSTM layers are configured with
$F=5$ filters of size $K=4$ to execute the convolutional operations with
stride $S=2$, while the fully connected layers of the detection branch are
configured with $U_{1}=64$, $U_{2}=32$ and $U_{3}=1$ units.

During optimization we aim to minimize the combined loss $L_{Dual}$ defined
as:

$L_{Dual}=w_{fix}\cdot L_{fix}+(1-w_{fix})\cdot L_{det}$ (12)

where $L_{fix}$ and $L_{det}$ correspond to the loss of the fixation and
detection prediction, respectively. The first is calculated using the
categorical cross entropy function stated in equation 10, while the second is
calculated using the weighted binary cross entropy defined in equation 14.

$w=y\cdot w_{1}+(1-y)\cdot w_{0}$ (13)
$L_{det}=w\cdot[y*log(\hat{y})+(1-y)\cdot log(1-\hat{y})]$ (14)

This final decision was made to address the imbalance in the detection data,
as the target is absent in half of the images, and the detection is positive
in a limited section of the scanpath sequences in the remaining images. To
prevent the model from biasing its prediction on the predominant class, we set
the relevance of the loss computed on the positive and negative detections
with the weights $w_{1}$ and $w_{0}$ set to $1.6$ and $0.7$, respectively. We
attained these values by, for the case of $w_{1}$ computing the multiplicative
inverse of the ratio of positive detections on the total number of detections,
and dividing it by 2; and, similarly, in the case of $w_{O}$ we perform the
same computations, but with the ration of negative detections.

To determine the optimal configuration for each model’s architecture, an
ablation study was done over the fovea size and the degree of importance
$w_{fix}$ of the fixation loss on the total loss. The former variable had its
values set to 50, 75 and 100 pixels, while the latter had its values set to
0.10, 0.25, 0.50, 0.75, 0.90.

###  IV-C Prediction

During the testing phase, both the fixation prediction module of the two-stage
pipeline and the dual task model were used to predict a scanpath sequence of
fixed length $l=7$, based on the training data, and we set the fixation point
at $t=0$ as the center cell of the discretized grid. Additionally, they were
both deployed with a beam search algorithm that selects the best $m$ fixation
points at each time step, where $m$ is the beam width hyper-parameter which we
set to 20, and appends them to the sequences they were generated from, while
saving the target detection prediction in the case of the dual task model. In
the next time step, the model runs for each of these $m$ predicted sequences,
to select the next best $m$ predictions. This algorithm selects the optimal
solution when considering the context given by the previously saved sequences,
while still having a linear complexity of $\mathcal{O}(l)=m\times l$.

In regards to the target presence detector, all models were deployed on the
scanpaths produced by the highest performing scanpath predictor. The baseline
classifier was tested similarly to the binary models, except that the input
cropped images are now constrained to a consistent size of $224\times 224$.
Due to the differences between the datasets, we adapt the ImageNet classes to
our targets by grouping some sibling sub-classes (e.g. the task bowl
corresponds to the union of the sub-classes mixing bowl and soup bowl), and we
do not classify the tasks fork, sink, and oven because they do not exist in
the ImageNet dataset. We classify a target as present when the corresponding
class has the highest classification score and as absent otherwise.

##  V Results

###  V-A Two Stage-Pipeline

####  V-A1 Fixation Prediction

Regarding the fixation prediction task, we evaluate our models on several
metrics. The first is Search Accuracy which is computed as the ratio of
sequences in which a fixation point selects a grid cell that intersects the
target’s bounding box. The second is Target Fixation Cumulative Probability
(TFP), which is plotted in figure 5, and presents the search accuracy attained
by each time step. From this plot we can additionally attain the TFP - Area
Under Curve (TFP-AUC) and the Probability Mismatch. The former computes the
area under the TFP curve and the latter is the sum of absolute differences
between the model’s TFP and the human’s observable data. Finally, we
additionally compute the Scanpath Ratio as the ratio between the sum of
euclidean distances between each fixation point and the distance from the
initial fixation to the center of the target’s bounding box.

Through the ablation study we conducted, we found that the high-level features
scanpath predictor from the two stage pipeline was able to achieve higher
search accuracy scores when using a one-hot task encoding and a Gaussian
ground-truth, as seen in figure 3, where the search scores are depicted in
box-plots grouped by each training setting. This new label representation was
able to outperform, with a search accuracy of 0.690, the traditional one-hot
encoding that attained a score of 0.650, bench-marking the model presented in
[27] on the COCO-Search18 dataset, which, previously, had only been able to
achieve a score of 0.621.

(a) Batch Size

(b) Fovea Size

(c) Task Encoding

(d) Ground Truth

Figure 3: Search accuracy box-plots of the single task High-level features’
models, grouped by their training settings, with the mean values and outliers
depicted as triangles and circles, respectively.

In figure 5 we have plotted the TFP curve of several highest performing model
configurations, in addition to the human search behavior111The Human TFP
refers to the human behavior observed in the training data set due to the fact
that the search fixations during the test dataset were withheld in the COCO-
Search18. All the remaining performances depicted were computed for the test
data split., and a random scanpath baseline model, which selected a random
human sequence from the train data split for the same search target class. We
observed a decrease in performance throughout time for all models through the
slopes of each function. Apart from the randomized generator, the models were
able to detect the majority of the targets by the second fixation step. Our
models, however, reach a performance plateau much sooner than the human
behavior, and barely detected any new targets throughout the last four
fixation points.

With the intent to better understand these low performance ceilings, we first
theorized that our models stopped searching after the second fixation point,
and throughout the remainder of the scanpath predicted solely the neighboring
grid cells. As depicted in figure 6, the predicted sequences do converge
throughout the time steps, however, this particularity is simply a replica of
the training data behavior.

When evaluating our model against two models presented in [17], which also use
panoptic features, the IRL and the BC-LSTM, we found that our models
outperform the latter across all criteria, although the former approach seems
to imitate human task-guided attention more accurately. The IRL utilizes
inverse reinforcement learning, while the BC-LSTM utilizes the same ConvLSTM
structure as our work, but instead of using it to directly predict the
scanpaths they use it in an intermediate step to update the feature maps.

(a) Mask Radius.

(b) Mask configuration.

(c) Depth of recurrent section.

(d) Activation Function.

Figure 4: Search accuracy box-plots of the single task Panoptic features’
models grouped by their training settings, with the mean values and outliers
depicted as triangles and circles, respectively.

In turn, the panoptic features’ fixation prediction model was found, through
figure 4, to reach higher accuracies with a single ConvLSTM layer in its
recurrent section, indicating an overfitting of the model as the number of
network parameters increases. In addition, while the highest search accuracy
score of 0.686 was attained with a cumulative mask of radius $r=1$, we could
not detect a direct correlation between these settings and the model’s
performance.

Surprisingly, we can also observe that using a final sigmoid activation
function greatly improved the model’s results, despite the fact that this
activation function is typically applied to single-class or multi-labeled
classifiers. This setup is only conceivable because the ground truth is
Gaussian. We interpreted the results as the optimizer perceiving the loss
function as the weighted sum of the prediction’s logarithm, as defined by the
equation 10, which caused it to boost the scores nearer the human fixation
positions. In contrast, when a final softmax activation is utilized, the model
turns the scores of the last hidden layers into class probabilities. Due to
the ground truth encoding, there is a saturation of the loss when every grid
cell is considered, as opposed to only examining the probability of the true
class in a one-hot encoding configuration, resulting in the model’s poor
performance.

Additionally, when comparing the single-task model leveraging panoptic
features to the models presented previously, we found, in figure 5, the models
with high-level feature maps are able to fixate targets much sooner even
though the panoptic features can lead to a similar search accuracy, leading to
a smaller TFP-AUC score of 3.259 and a higher probability mismatch of 1.514.
In addition, this approach is much less efficient as scanpaths generated
through panoptic features travel a much greater distance, as depicted in
figure 6, leading to a a scanpath ratio score of only 0.463.

TABLE I:  Performance evaluation of several best performing models (rows) based on several Fixation Prediction metrics (columns). | Search | TFP-AUC $\uparrow$ | Probability | Scanpath
---|---|---|---|---
| Accuracy $\uparrow$ | Mismatch $\downarrow$ | Ratio $\uparrow$
Human | 0.990 | 5.200 | - | 0.862
High-Level Features - One-hot GT | 0.650 | 3.068 | 1.727 | 0.753
High-Level Features - Gaussian GT | 0.690 | 3.413 | 1.360 | 0.727
Panoptic Features | 0.686 | 3.259 | 1.514 | 0.463
Dual - Architecture A | 0.719 | 3.496 | 1.263 | 0.808
Dual - Architecture C | 0.701 | 3.446 | 1.320 | 0.791
IRL | N/A | 4.509 | 0.987 | 0.826
BC-LSTM | N/A | 1.702 | 3.497 | 0.406
Random Scanpath | 0.235 | 1.150 | 3.858 | -
Figure 5: Search accuracy achieved by each model throughout scanpaths, with
means and standard errors first computed over target classes. Figure 6:
Euclidean distances between fixation points.

####  V-A2 Target Detection

For the target detection task we used accuracy, precision and recall as
metrics. The fine-tuned classifiers with foveation radius of 50 pixels had the
maximum performance for all measures, with a mean accuracy, precision and
recall of 82.1%, 86.9% and 75.4%, respectively, whereas the configuration of
75 pixels achieved the lowest accuracy and recall. Regarding each individual
task, the models with the best performance were those detecting the bottle and
stop sign classes, with the first dominating in accuracy and recall with
scores of 92.7% and and 93.2%, respectively, and the second in precision
scoring 100%, all averaged across fovea settings. Additionally, we may remark
that the knife, microwave, and potted plant jobs had the most variation in
metrics. The first model, with a fovea of 50 pixels, has a recall score of
80.0% but an accuracy score of 48.9%, meaning that the model inaccurately
anticipated the existence of the target more than half of the time while
rarely reporting positive values as target present. Contrarily, the potted
plant model with a fovea of 50 pixels scored 97.9% for precision and 52.3% for
recall.

In turn, the baseline pre-trained model attained incredibly high ranges across
metrics for the vast majority of the classes, as seen in figure 7(b). Although
its average precision is similar to the previous models, its recall is
incredibly low, averaging the worst in the foveation size with 50 pixels
context at 28.1%. We attribute this to the classes groupings and the
differences between the training and testing datasets. We additionally noted a
slight increase of the mean scores with the increase of the foveation radius,
across all metrics, with its highest accuracy, precision and recall scores
being 64.0%, 85.7% and 34.3%, respectively.

(a) Fine-Tuned Classifier

(b) Baseline Classifier

Figure 7: Performances of the fine-tuned and baseline target detectors in
terms of accuracy, precision and recall, for the fovea size setting of 50
pixels.

Additionally, to better comprehend the performance of our model, we compare it
to various state-of-the-art benchmarks. The VGG-16 model, trained on the
ImageNet dataset, is able to reach a 71.3% top-1 accuracy when classifying
1,000 classes on the validation split of the same dataset [32], which our
fine-tuned model outperforms in this problem. This increase in performance is
likely attributable to the nature of our task, as we are performing binary
classification whereas the benchmark presented for the VGG-16 was performed on
multi-class classification of one thousand classes with some classes having
highly nuanced visual distinctions, such as mud turtle and box turtle being
two distinct classes.

In [27], three separate approaches were taken, the first being a K nearest
neighbors approach, where the feature maps of each fixation point, during
testing, are evaluated against the labeled feature maps obtained from the
training dataset of the same task. The prediction is then computed according
to the majority label of the closest neighbors. The second approach is a
simple method in which a comparison is made between the test sample and the
feature maps averaged across the training data. This comparison is made
through a loss, which, in the cases it falls bellow a pre-defined threshold,
detects the target as present. Lastly, the same approach as our baseline was
taken. These three achieved accuracies of 77.8%, 30.5% and 70.24%,
respectively. The first two results may be attributable, once more, to the
lower size of the dataset used in their study, while the last model
outperforms our baseline probably due to a better match between class
equivalencies, resulting in a higher degree of similarity between the training
and testing datasets.

###  V-B Dual Task

Finally, we evaluate the dual task models in their performances in both tasks.
In figure 8(a), we can see that architecture B fared the worst, and,
architecture A outperforms more than half of the models with architecture C,
in addition to achieving the highest search accuracy of 73.4% when configured
with a fovea size of 100 pixels and a fixation loss’ weight of 0.75.

In terms of its detection performance, the dual task model with architecture
C, a fovea size of 50, and $w_{fix}$ set to 0.90 achieved the highest target
presence detection rate of 68.7%. In figure 8(c), considering the quartiles
and upper limit of its performance, the architecture of design C is deemed to
be the most effective. Regarding the weight of the fixation loss, we can also
see that models trained with bigger values obtained a larger interquartile
range than models trained with smaller values.

(a) Architecture.

(b) Fixation loss’ weight.

(c) Architecture.

(d) Fixation loss’ weight.

Figure 8:  Box-plots of the Search accuracy, on the top, and Detection
accuracy, on the bottom, grouped by their configuration, with the mean values
and outliers depicted as triangles and circles, respectively.

In addition, we noted that the best scanpath prediction model only achieved a
detection accuracy of 49.7% while the best target presence predictor only
achieved a search accuracy of 63.9%. In light of this, we also considered the
average of both of these metrics to determine which models were more capable
of performing both tasks simultaneously. The majority of the time, design A
earned a higher score than design C, while design B ranked the lowest.
Regarding the remaining parameters, we can only observe that a bigger fovea
radius led to higher average scores, and a higher fixation loss’ weight
resulted in a better top score, with the exception of setting $w_{f}ix=0.25$.
The model configured with architecture C, a fovea size of 75 pixels, and
$wfix$ set to 0.90 achieved the top score of 67.7% with search and detection
accuracies of 70.1% and 65.3%, respectively.

During testing, we discovered a trend in the target detection predictions from
models with architecture A, which is exemplified in figure 9, which plots the
confusion matrices of the temporal predictions of the models with the highest
overall scores for architectures A and C. During the second half of the gaze
sequence, the former model hardly predicts that the target is absent,
indicating that the recurrent nature of the detection branch biases the model
to output the pattern learned from the training data that, for target present
images, the object is detected a vast majority of the time, as stated
previously. In contrast, we can see that the model with architecture C seeks
to differentiate between target present and target missing fixations at later
time steps.

(a) Architecture A.

(b) Architecture C.

Figure 9: Target detection confusion matrices, in percentage, separated by
time step for the best overall models with fixation first, on top, and
concatenation, bellow, architectures.

Finally, when comparing the fixation prediction performance of this approach
to that of previously presented models, the dual approach led to a higher
search accuracy, as seen in figure I, resulting in a higher TFP-AUC score of
3.496 and 3.446 and a lower probability mismatch of 1.263 and 1.320 for the
overall best models with architectures A and C, respectively. In addition, the
two models exhibit a higher search efficiency with scanpath ratios of 0.808
and 0.791, respectively, because, as shown in 6, the model with architecture C
travels a shorter distance during the first time step and the model with
architecture A is able to maintain an efficient search throughout the
remainder of the scanpath.

##  VI Conclusions

In this thesis, we described two methods for predicting the presence or
absence of a target in an image with foveated context. In the first, a two-
stage pipeline system, the fixation prediction module produced the best
results while using high-level feature maps as image representations, both in
terms of search accuracy and search efficiency. In addition, we discovered
that the usage of a Gaussian ground truth label encoding, which is our first
contribution, enhanced search accuracy. This novel representation captures the
spatial structure of the problem by not only encouraging predictions to the
exact discretized human fixation positions, but also favoring attempts to
cells near these locations.

In the fixation prediction models with panoptic features, a cumulative mask
representation and a final sigmoid activation function were shown to be
advantageous. The latter is only conceivable when using a Gaussian ground
truth label. This setup imposes that the categorical cross entropy loss
performs a weighted sum of the loss contribution of each cell, where positions
closer to the correct location are assigned a higher weight. However, while
these image representations achieved similar search accuracy to the high-level
feature maps, they presented a lower search efficiency. This suggests that the
ConvLSTM is better suited for high-level feature inputs, but the IRL learning
strategy was able to harness the full potential of the panoptic features by
predicting both the action and state representations at each time step.

Two classifiers performed the target presence prediction portion of the two-
stage pipeline model. A pre-trained VGG-16 was used as a baseline for the dual
task model, and a fine-tuning approach set the performance ceilings for the
detection task.

The final contribution of our thesis is a dual-task model that executes both
tasks concurrently while enabling information sharing between them by
executing a common input transformation and establishing linking channels
throughout each task branch. This multi-task approach improved search
precision when the task prediction branch initiated the predictions, i.e. in
designs A and C. However, the former suffered a reduction in detection
accuracy, whilst the latter achieved the maximum score when compared to our
baseline method. We additionally concluded that the proposed approach of
employing a ConvLSTM layer during the detection branch resulted in the model
basing its predictions on the temporal pattern of target detection, as opposed
to utilizing the input visual features as the basis for its prediction.

##  VII Future Work

The gap in search accuracy performance between the human performance ceiling
and our best results reflects the challenge of duplicating task-guided
attention at the start of the scanpath. As the model’s target fixation reaches
a plateau earlier than its training data, future research should concentrate
on increasing the fixation prediction of the first time-steps and not so much
the last. In the testing phase, it is important to test a new fixation
prediction algorithm that gives greater weight to sequences with higher scores
at the beginning of the scanpath.

Regarding the target’s presence prediction during the dual-task approach, we
have demonstrated that the use of a recurrent layer biases the model towards
temporal patterns of target detection. An alternative solution would be for
the task branch to generate the input image simultaneously, simulating an
encoder-decoder, so as to require the model to maintain its knowledge of high-
level features in its hidden states.

Future research should also investigate a visual transformer-based design, as
it has shown promising results in similar image classification and goal-
directed search tasks.

## References

  * [1] Szeliski, R., 2011. Computer vision algorithms and applications. Springer, London; New York.
  * [2] Borji, A. and Itti, L., 2012. State-of-the-art in visual attention modeling. IEEE transactions on pattern analysis and machine intelligence, 35(1), pp.185-207.
  * [3] Bandera, C. and Scott, P.D., 1989, November. Foveal machine vision systems. In Conference Proceedings., IEEE International Conference on Systems, Man and Cybernetics (pp. 596-599). IEEE.
  * [4] Fukushima, K., 1988. Neocognitron: A hierarchical neural network capable of visual pattern recognition. Neural networks, 1(2), pp.119-130.
  * [5] Krizhevsky, A., Sutskever, I. and Hinton, G.E., 2017. Imagenet classification with deep convolutional neural networks. Communications of the ACM, 60(6), pp.84-90.
  * [6] LeCun, Y., Bottou, L., Bengio, Y. and Haffner, P., 1998. Gradient-based learning applied to document recognition. Proceedings of the IEEE, 86(11), pp.2278-2324.
  * [7] Akbas, E. and Eckstein, M.P., 2017. Object detection through search with a foveated visual system. PLoS computational biology, 13(10), p.e1005743.
  * [8] James, W., 1890. The principles of psychology, Vol. 1. Henry Holt and Co.
  * [9] Corbetta, M. and Shulman, G.L., 2002. Control of goal-directed and stimulus-driven attention in the brain. Nature reviews neuroscience, 3(3), pp.201-215.
  * [10] Yarbus, A.L., 2013. Eye movements and vision. Springer.
  * [11] Ngo, T. and Manjunath, B.S., 2017, September. Saccade gaze prediction using a recurrent neural network. In 2017 IEEE International Conference on Image Processing (ICIP) (pp. 3435-3439). IEEE.
  * [12] Graves, A., 2012. Long short-term memory. Supervised sequence labelling with recurrent neural networks, pp.37-45.
  * [13] Kreiman, G. and Zhang, M., 2018. Finding any Waldo: zero-shot invariant and efficient visual search.
  * [14] Nunes, A., Figueiredo, R. and Moreno, P., 2020, June. Learning to search for objects in images from human gaze sequences. In International Conference on Image Analysis and Recognition (pp. 280-292). Springer, Cham.
  * [15] Shi, X., Chen, Z., Wang, H., Yeung, D.Y., Wong, W.K. and Woo, W.C., 2015. Convolutional LSTM network: A machine learning approach for precipitation nowcasting. Advances in neural information processing systems, 28.
  * [16] Chen, Y., Yang, Z., Ahn, S., Samaras, D., Hoai, M., and Zelinsky, G. (2021). COCO-Search18 Fixation Dataset for Predicting Goal-directed Attention Control. Scientific Reports, 11 (1), 1-11, 2021.
  * [17] Yang, Z., Huang, L., Chen, Y., Wei, Z., Ahn, S., Zelinsky, G., Samaras, D., and Hoai, M. (2020). Predicting Goal-directed Human Attention Using Inverse Reinforcement Learning. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (pp. 193-202).
  * [18] Kirillov, A., He, K., Girshick, R., Rother, C. and Dollár, P., 2019. Panoptic segmentation. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (pp. 9404-9413).
  * [19] Kirillov, A., Girshick, R., He, K. and Dollár, P., 2019. Panoptic feature pyramid networks. In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition (pp. 6399-6408).
  * [20] He, K., Gkioxari, G., Dollár, P. and Girshick, R., 2017. Mask r-cnn. In Proceedings of the IEEE international conference on computer vision (pp. 2961-2969).
  * [21] Chen, C., Wu, Y., Dai, Q., Zhou, H.Y., Xu, M., Yang, S., Han, X. and Yu, Y., 2022. A Survey on Graph Neural Networks and Graph Transformers in Computer Vision: A Task-Oriented Perspective. arXiv preprint arXiv:2209.13232.
  * [22] Jonnalagadda, A., Wang, W.Y., Manjunath, B.S. and Eckstein, M.P., 2021. Foveater: Foveated transformer for image classification. arXiv preprint arXiv:2105.14173.
  * [23] Chakraborty, S., Gupta, R., Ma, K., Govind, D., Sarder, P., Choi, W.T., Mahmud, W., Yee, E., Allard, F., Knudsen, B. and Zelinsky, G., 2022. Predicting the Visual Attention of Pathologists Evaluating Whole Slide Images of Cancer. In International Workshop on Medical Optical Imaging and Virtual Microscopy Image Analysis (pp. 11-21). Springer, Cham.
  * [24] Bhattacharya, M., Jain, S. and Prasanna, P., 2022. RadioTransformer: A Cascaded Global-Focal Transformer for Visual Attention-guided Disease Classification. arXiv preprint arXiv:2202.11781.
  * [25] Simonyan, K. and Zisserman, A., 2014. Very deep convolutional networks for large-scale image recognition. arXiv preprint arXiv:1409.1556.
  * [26] Deng, J., Dong, W., Socher, R., Li, L.J., Li, K. and Fei-Fei, L., 2009, June. Imagenet: A large-scale hierarchical image database. In 2009 IEEE conference on computer vision and pattern recognition (pp. 248-255). Ieee.
  * [27] B. Cabarrão, Learning to search for objects in foveal images using deep learning, Master’s thesis, Universidade de Lisboa \- Instituto Superior Técnico, 2022
  * [28] Lin, T.Y., Maire, M., Belongie, S., Hays, J., Perona, P., Ramanan, D., Dollár, P. and Zitnick, C.L., 2014, September. Microsoft coco: Common objects in context. In European conference on computer vision (pp. 740-755). Springer, Cham.
  * [29] Cerf, M., Harel, J., Einhäuser, W. and Koch, C., 2007. Predicting human gaze using low-level saliency combined with face detection. Advances in neural information processing systems, 20.
  * [30] Judd, T., Ehinger, K., Durand, F. and Torralba, A., 2009, September. Learning to predict where humans look. In 2009 IEEE 12th international conference on computer vision (pp. 2106-2113). IEEE.
  * [31] Kingma, D.P. and Ba, J., 2014. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980.
  * [32] Keras applications - performances, https://keras.io/api/applications/, accessed: 2022-09-13.


# Strong regulatory graphs

Patric Gustafsson
Department of Information Technologies This work was initiated as part of
Patric Gustafsson’s master thesis project at Åbo Akademi University (ÅA) in
2019. Patric is no longer affiliated with ÅA. Åbo Akademi University  Finland
<EMAIL_ADDRESS>Ion Petre
Department of Mathematics and Statistics  University of Turku  Finland
National Institute for Research and Development for Biological Sciences
Romania
<EMAIL_ADDRESS>
###### Abstract

Logical modeling is a powerful tool in biology, offering a system-level
understanding of the complex interactions that govern biological processes. A
gap that hinders the scalability of logical models is the need to specify the
update function of every vertex in the network depending on the status of its
predecessors. To address this, we introduce in this paper the concept of
strong regulation, where a vertex is only updated to active/inactive if all
its predecessors agree in their influences; otherwise, it is set to ambiguous.
We explore the interplay between active, inactive, and ambiguous influences in
a network. We discuss the existence of phenotype attractors in such networks,
where the status of some of the variables is fixed to active/inactive, while
the others can have an arbitrary status, including ambiguous.

###### keywords:

Biomodeling, interaction networks, regulatory graphs, Boolean networks,
phenotype attractors.

Strong regulatory graphs

##  1 Introduction

Logical modeling is a well established mathematical approach to biology,
offering a simple, intuitive way to understand complex systems. Biological
networks can be represented as graphs and experimental data integrated in them
to yield a system-level view of key regulators, dynamical patterns, and
response to various changes. There are many successful applications of logical
modeling such as the identification of drug targets and treatments [1],
modeling of regulatory networks [2], of signaling pathways [3], of cell growth
and apoptosis [4], applications in immunology [5], and many others. A logical
model is discrete (often Boolean), with its dynamics defined as the result of
the interplay of the influences among its variables. Such influences are
modeled as directed edges between the variables, sometimes with the annotation
of the activation/inhibition nature of the interaction, see [6] for a review
on Boolean modeling in systems biology.

Constructing large-scale models has been made possible in recent years by the
availability of data respositories such as KEGG ([7]), OmniPath ([8]),
InnateDB ([9]), SIGNOR ([10]) and DrugBank ([11]). Using such resources has
made it possible to build models consisting of thousands of variables and
interactions [12, 13]. Integrating such large scale data into logical or
Boolean models is problematic because of the need to identify the update
functions for each variable of the model. These functions specify the precise
conditions under which a variable changes its status, e.g. between active and
inactive, depending on the status of its predecessors/regulators in the
network. The difficulty is that the data is almost never enough to identify
the precise nature of these functions. Instead, the modeler often postulates
the type of update functions in the network (propositional formulas ([14]),
threshold conditions ([15]), multi-valued functions ([16]), etc.) The
specifics of how to choose each update function are based on previous
literature and to some extent, to the choice of the modeler, subject to
experimentation with various options.

The approach we take in this paper is motivated by the concept of controlling
a network: choose a user intervention in the network (e.g., through fixing
some variables to some constants), so that the model converges to a certain
desired attractor. Network controllability has found interesting applications
in biology, e.g. in drug repurposing [12, 13]. Tracking the effect of the user
interventions throughout the network is however difficult when the variables
of the model are under conflicting influences, with some regulators pushing
for their activation, while others for their inhibition. We propose in this
paper a notion of _strong regulatory graphs_ where the update of a node’s
status is determined to be active/inactive _only if its predecessors concur in
their influences_. Otherwise it is set to be _ambiguous_ , meaning that it
could be both active and inactive, depending on the precise (numerical) setup
of its predecessors and of their influences, which is in practice very
difficult to determine. This leads to an intricate interplay between
ambiguous, negative and positive influences in our framework. Interesting
questions about the spread of ambiguity in the network can be asked, e.g. in
terms of the existence of phenotype attractors. We introduce in this paper a
simple mathematical formalization of the concept of _strong regulation_. We
discuss the phenotype attractor problem for strong regulatory graphs: whether
attractors exist in which some of the variables are constant (they follow a
given active/inactive phenotype), while the others can be arbitrary, even
ambiguous.

##  2 Strong regulatory graphs (SRG)

The model proposed here builds on the regulatory graphs structure [17, 18, 19,
16] and adds to them the concept of strong regulation. Choosing how to update
the status of a node is obvious if its active regulators concur in their
influences: if they all exert an activation regulation on the node, then it
gets _activated_ (denoted as $1$), and if all exert an inhibition regulation
on it, then it gets _inhibited_ (denoted as $-1$). A node under conflicting
regulation, with some active regulators trying to activate it and others to
inhibit it, is less clear how to update. The typical approach is to postulate
an update function, Boolean or multi-valued, that sets its status to active or
inactive, based on a specific schema of the status of its predecessors. The
main difficulty of these approaches is that the experimental data and the
level of detail in the model are rarely sufficient to specify such detailed
update functions. Our proposal instead is that under such conflicting the
status of the node should simply be recognized as being _ambiguous_ (denoted
as $0$). This means that the node may potentially be active or inactive, but
the level of detail in the model is not enough to specify it with confidence
either way. The ambiguity of a node may well cascade down through the system,
in that the regulation it exerts on its successor in the regulatory graph may
be active or not. Still, an ambiguous regulation, may well be resolved by
another _active_ regulation on the same node. This opens the possibility to
reason about the dynamic interplay between activation, inhibition, and
ambiguity, where the semantics of activation and inhibition is in its
_strongest_ possible form to mean that the node has been _unanimously_
activated/inhibited.

###### Definition 2.1

A strong regulatory graph (SRG) is an edge-labeled graph $G=(V,E)$ where
$V=\\{v_{1},v_{2},\ldots,v_{n}\\}$ is the finite set of _vertices_ and
$E\subseteq V\times V$ is the set of directed edges. The labeling of the edges
is done through the partition $\\{E_{+},E_{-}\\}$ of $E$. We say that edge
$(u,v)\in E_{+}$ is an _activating edge_ and that $u$ is an _activator_ of
$v$. We say that edge $(u,v)\in E_{-}$ is an _inhibiting edge_ and that $u$ is
an _inhibitor_ of $v$.

A state of the SRG is a function $\alpha:V\rightarrow\\{-1,0,1\\}$ whose
intended meaning is that $\alpha(v)=-1$ if $v$ is inactive in state $\alpha$,
$\alpha(v)=1$ if $v$ is active, and $\alpha(v)=0$ if $v$ is ambiguous. If
$\alpha(u)=1$ ($\alpha(u)=0$), we say that $u$ is an _activator_ (_potential
activator_ , resp.) of $v$ in $\alpha$ if $(u,v)\in E_{+}$ and that it is an
_inhibitor_ (_potential inhibitor_ , resp) of $v$ in $\alpha$ if $(u,v)\in
E_{-}$.

The dynamics of a strong regulatory graph $G=(V,E)$ is defined through a
state-transition system. The state is given by the -1/0/1 activation status of
all nodes in the graph. The set of states is thus $\\{-1,0,1\\}^{|V|}$, also
represented as a function from $V$ to $\\{-1,0,1\\}$. Given a state, a
transition will indicate the change in the activation status of the nodes in
the graph. Depending on whether we consider a synchronous update of the
activation status (all nodes are simultaneously updated) or an asynchronous
update (only one node is updated at a time), we obtain either a deterministic
or a non-deterministic finite state transition system. For simplicity, we only
discuss in this paper the synchronous activation update and the dynamics of
strong regulatory graph in terms of deterministic state-transition systems.
The asynchronous update strategy leads to a similar discussion.

In defining the update rule for the status of a node, we follow proposal for
the concept of _strong regulation_ , where a node is activated/inhibited if
all its regulators concur in their influences on the node. In case of
conflicting influences, the status of the node is set to ‘ambiguous’. The
influence of a regulator whose status is ambiguous is quite subtle, depending
on the status of the other regulators and of the node itself. We discuss the
intuition of the update rule below and then give it a formal definition.

If there are no potentially active inhibitors and no potentially active
activators, then the vertex preserves its activation status. A node $v$ is set
to be active after in two situations:

  * •

$v$ was active and none of its inhibitors were potentially active, i.e.,
active or ambiguous (in this case, the node remains active, absent any
potential or active inhibition regulation), or

  * •

at least one activator of $v$ was active and none of its inhibitors were
potentially active (in this case, the regulation on the node is un-ambiguous
towards activation).

Similarly, a node $v$ becomes inactive after an application of the update rule
in two situations:

  * •

$v$ was inactive and none of its activators were potentially active, i.e.,
active or ambiguous, or

  * •

at least one inhibitor of $v$ was active and none of its activators were
potentially active.

If the node is under conflicting influences, its status becomes ambiguous.
Here are four cases when this happens:

  * •

at least one activator and at least one inhibitor of $v$ were potentially
active; this reflects a node under potentially conflicting regulation;

  * •

$v$ was active and one of its inhibitors was ambiguous (regardless of the
status of its activators); this reflects the ambiguity of whether there is an
inhibitory influence on $v$ or not, which may either leave the node active, or
switch to inactive;

  * •

$v$ was inactive and one of its activators was ambiguous (regardless of the
status of its inhibitors); this reflects the ambiguity of whether there is an
activating influence on $v$ or not, which may either leave the node inactive,
or switch it to active;

  * •

$v$ was ambiguous and none of its regulators is active (i.e., they are either
inactive or ambiguous); this reflects the situation where the status of $v$
cannot be clarified because the status of its regulators is either inactive or
ambiguous.

For a vertex $v\in V$ we define its set
$\operatorname{\mathsf{Reg}}_{+}(\alpha,v)$ of activators and its set
$\operatorname{\mathsf{Reg}}_{-}(\alpha,v)$ of inhibitors in state
$\alpha:V\rightarrow\\{-1,0,1\\}$ as follows:

$\displaystyle\operatorname{\mathsf{Reg}}_{+}(\alpha,v)$
$\displaystyle=\\{\alpha(u)\mid(u,v)\in E_{+},\alpha(u)\in\\{0,1\\}\\},$
$\displaystyle\operatorname{\mathsf{Reg}}_{-}(\alpha,v)$
$\displaystyle=\\{\alpha(u)\mid(u,v)\in E_{-},\alpha(u)\in\\{0,1\\}\\}.$

Set $\operatorname{\mathsf{Reg}}_{+}(\alpha,v)$
($\operatorname{\mathsf{Reg}}_{-}(\alpha,v)$) contains $1$ if $v$ has at least
one active activator (inhibitor, resp.) in state $\alpha$ and it contains $0$
if $v$ has at least one ambiguous activator (inhibitor, resp.) in state
$\alpha$.

We also define the reflexive extension of $\operatorname{\mathsf{Reg}}$ to
include the status of the vertex itself. The reflexive extension captures a
sort of “inertia” of the dynamics (defined below) where a node may preserve
its status in the absence of active regulators.

$\displaystyle\overline{\operatorname{\mathsf{Reg}}}_{+}(\alpha,v)$
$\displaystyle=\operatorname{\mathsf{Reg}}_{+}(\alpha,v)\cup\\{\alpha(v)\mid\alpha(v)\in\\{0,1\\}\\},$
$\displaystyle\overline{\operatorname{\mathsf{Reg}}}_{-}(\alpha,v)$
$\displaystyle=\operatorname{\mathsf{Reg}}_{-}(\alpha,v)\cup\\{-\alpha(v)\mid\alpha(v)\in\\{-1,0\\}\\},$

Note that an inactive node $v$, i.e. $\alpha(v)=-1$, contributes $1$ to
$\overline{\operatorname{\mathsf{Reg}}}_{-}(\alpha,v)$, consistent with our
suggestion of “inertia”: absent any active regulation on node $v$, it will
remain remain inactive.

###### Definition 2.2

Let $G=(V,E)$ be a strong regulatory graph. Its dynamics is given through the
state-transition system $(S_{G},T_{G})$, where
$S_{G}=\left\\{\gamma:V\rightarrow\\{-1,0,1\\}\right\\}$, and $T_{G}\subseteq
S_{G}\times S_{G}$ defined as follows: $(\alpha,\beta)\in T_{G}$ (denoted also
$\alpha\rightarrow_{G}\beta$) if

$\beta(v)=\begin{cases}1,\text{ if
}1\in\overline{\operatorname{\mathsf{Reg}}}_{+}(\alpha,v)\text{ and
}\operatorname{\mathsf{Reg}}_{-}(\alpha,v)=\emptyset,\\\ -1,\text{ if
}1\in\overline{\operatorname{\mathsf{Reg}}}_{-}(\alpha,v)\text{ and
}\operatorname{\mathsf{Reg}}_{+}(\alpha,v)=\emptyset,\\\ 0,\text{
otherwise.}\end{cases}$

|  |
---|---|---
(a) |  | (b)
Figure 1: Two simple regulatory graphs [16]. The vertices are shown with
rectangles, the activation edges with green pointed arrows and the inhibition
edges with red blunt arrows.

###### Example 2.3

We discuss two regulatory graphs introduced in [16] and shown in Figure 1. The
set of vertices is $V=\\{A,B,C\\}$, the activation edges are shown with green
pointed arrows and the inhibition edges with red blunt arrows.

For the graph in Figure 1(a), consider the state $\alpha_{1}=(-1,1,1)$. Then

  * •

$\operatorname{\mathsf{Reg}}_{+}(\alpha_{1},A)=\\{1\\}$ (since $(C,A)$ is an
activation edge and $\alpha_{1}(C)=1$) and
$\operatorname{\mathsf{Reg}}_{-}(\alpha_{1},A)=\\{1\\}$ (since $(B,A)$ is an
inhibition edge and $\alpha_{1}(B)=1$). Also,
$\overline{\operatorname{\mathsf{Reg}}}_{+}(\alpha_{1},A)=\\{1\\}$ and
$\overline{\operatorname{\mathsf{Reg}}}_{-}(\alpha_{1},A)=\\{1\\}$ (since
$\alpha_{1}(A)=-1$).

  * •

$\operatorname{\mathsf{Reg}}_{+}(\alpha_{1},B)=\emptyset$ (since the only
incoming activation edge into $B$ is $(A,B)$ but $\alpha_{1}(A)=-1$) and
$\operatorname{\mathsf{Reg}}_{-}(\alpha_{1},B)=\emptyset$ (since there are no
incoming inhibition edges into $B$). Also,
$\overline{\operatorname{\mathsf{Reg}}}_{+}(\alpha_{1},B)=\\{1\\}$ and
$\overline{\operatorname{\mathsf{Reg}}}_{-}(\alpha_{1},B)=\emptyset$ (since
$\alpha_{1}(B)=1$).

  * •

$\operatorname{\mathsf{Reg}}_{+}(\alpha_{1},C)=\operatorname{\mathsf{Reg}}_{-}(\alpha_{1},C)=\emptyset$
(since there are no incoming edges into $C$). Also,
$\overline{\operatorname{\mathsf{Reg}}}_{+}(\alpha_{1},$ $C)=\\{1\\}$ and
$\overline{\operatorname{\mathsf{Reg}}}_{-}(\alpha_{1},C)=\emptyset$ (since
$\alpha_{1}(C)=1$).

This means that $\alpha_{1}\rightarrow(0,1,1)$. This corresponds well to our
intuition for the meaning of the strong regulatory graphs: $A$ is under
conflicting regulation and it becomes set to _ambiguous_ , $B$ and $C$ are
under no active regulation and they preserve their status. Denote
$\alpha_{2}=(0,1,1)$. Then

  * •

$\operatorname{\mathsf{Reg}}_{+}(\alpha_{2},A)=\\{1\\}$ (since $(C,A)$ is an
activation edge and $\alpha_{2}(C)=1$) and
$\operatorname{\mathsf{Reg}}_{-}(\alpha_{1},A)=\\{1\\}$ (since $(B,A)$ is an
inhibition edge and $\alpha_{2}(B)=1$). Also,
$\overline{\operatorname{\mathsf{Reg}}}_{+}(\alpha_{2},A)=\\{0,1\\}$ and
$\overline{\operatorname{\mathsf{Reg}}}_{-}(\alpha_{1},A)=\\{0,1\\}$ (since
$\alpha_{2}(A)=0$).

  * •

$\operatorname{\mathsf{Reg}}_{+}(\alpha_{2},B)=\\{0\\}$ (since $(A,B)$ is an
activation edge and $\alpha_{2}(A)=0$) and
$\operatorname{\mathsf{Reg}}_{-}(\alpha_{2},B)=\emptyset$ (since there are no
incoming inhibition edges into $B$). Also,
$\overline{\operatorname{\mathsf{Reg}}}_{+}(\alpha_{2},B)=\\{0,1\\}$ and
$\overline{\operatorname{\mathsf{Reg}}}_{-}(\alpha_{2},B)=\emptyset$ (since
$\alpha_{2}(B)=1$).

  * •

$\operatorname{\mathsf{Reg}}_{+}(\alpha_{2},C)=\operatorname{\mathsf{Reg}}_{-}(\alpha_{2},C)=\emptyset$
(since there are no incoming edges into $C$). Also,
$\overline{\operatorname{\mathsf{Reg}}}_{+}(\alpha_{2},$ $C)=\\{1\\}$ and
$\overline{\operatorname{\mathsf{Reg}}}_{-}(\alpha_{2},C)=\emptyset$ (since
$\alpha_{2}(C)=1$).

Consequently, $\alpha_{2}\rightarrow\alpha_{2}$. The intuition for the update
on $A$ and $C$ remains the same. For $B$, there is a potential activation
regulation coming from the ambiguous status of its activator $A$, but since
$B$ is active in $\alpha_{2}$, it will remain so in the next state.

For the graph in Figure 1(b), the following are valid state transitions:

  * •

$(-1,1,-1)\rightarrow(1,1,-1)\rightarrow(1,1,-1)$,

  * •

$(1,-1,-1)\rightarrow(1,1,-1)$,

  * •

$(-1,-1,1)\rightarrow(-1,-1,1)$ and

  * •

$(1,-1,1)\rightarrow(-1,1,1)\rightarrow(0,1,1)\rightarrow(0,1,1)$.

The state transition graphs for these two examples are in Figure 2.

---
(a)
(b)
Figure 2: The state transition graphs of the regulatory graphs in Example 2.3
and Figure 1.

The following observation shows that a vertex preserves its activation status,
in the absence of potentially active predecessors.

###### Lemma 2.4

Let $G=(V,E)$ be a strong regulatory graph, $v\in V$ in state $\alpha$, with
$\operatorname{\mathsf{Reg}}_{+}(\alpha,v)=\operatorname{\mathsf{Reg}}_{-}(\alpha,v)=\emptyset$.
If $\alpha\rightarrow_{G}\beta$, then $\beta(v)=\alpha(v)$.

###### Proof 2.5

If $\alpha(v)=-1$, then
$\overline{\operatorname{\mathsf{Reg}}}_{-}(\alpha,v)=\\{1\\}$ and so,
$\beta(v)=-1$. If $\alpha(v)=1$, then
$\overline{\operatorname{\mathsf{Reg}}}_{+}(\alpha,v)=\\{1\\}$ and so,
$\beta(v)=1$. If $\alpha(v)=0$, then
$\overline{\operatorname{\mathsf{Reg}}}_{-}(\alpha,v)=\overline{\operatorname{\mathsf{Reg}}}_{+}(\alpha,v)=\\{0\\}$
and so, $\beta(v)=0$.

We clarify now the conditions under which a vertex $v$ is updated to an
‘ambiguous’ state. Following Lemma 2.4, we only focus on the case when the
vertex is under some (potential) regulators, i.e., when either
$\operatorname{\mathsf{Reg}}_{+}(\alpha,v)\neq\emptyset$, or
$\operatorname{\mathsf{Reg}}_{-}(\alpha,v)\neq\emptyset$, as otherwise $v$
maintains its status unchanged. The theorem shows that our definition covers
exactly the intuition offered just before Definition 2.1.

###### Theorem 2.6

Let $G=(V,E)$ be an SRG and $v\in V$ in state $\alpha$, with
$\operatorname{\mathsf{Reg}}_{+}(\alpha,v)\cup\operatorname{\mathsf{Reg}}_{-}(\alpha,v)\neq\emptyset$.
If $\alpha\rightarrow_{G}\beta$, we have $\beta(v)=0$ if and only if

  1. (i)

$\operatorname{\mathsf{Reg}}_{+}(\alpha,v)\neq\emptyset$ and
$\operatorname{\mathsf{Reg}}_{-}(\alpha,v)\neq\emptyset$, or

  2. (ii)

$\overline{\operatorname{\mathsf{Reg}}}_{+}(\alpha,v)=\\{0\\}$ and
$\operatorname{\mathsf{Reg}}_{-}(\alpha,v)=\emptyset$, or

  3. (iii)

$\overline{\operatorname{\mathsf{Reg}}}_{-}(\alpha,v)=\\{0\\}$ and
$\operatorname{\mathsf{Reg}}_{+}(\alpha,v)=\emptyset$.

###### Proof 2.7

$\displaystyle\beta(v)=0\Leftrightarrow\beta(v)\neq 1\land\beta(v)\neq-1$
$\displaystyle\Leftrightarrow\ $
$\displaystyle(1\not\in\overline{\operatorname{\mathsf{Reg}}}_{+}(\alpha,v)\lor\operatorname{\mathsf{Reg}}_{-}(\alpha,v)\neq\emptyset)\land(1\not\in\overline{\operatorname{\mathsf{Reg}}}_{-}(\alpha,v)\lor\operatorname{\mathsf{Reg}}_{+}(\alpha,v)\neq\emptyset)$
$\displaystyle\Leftrightarrow\ $
$\displaystyle(1\not\in\overline{\operatorname{\mathsf{Reg}}}_{+}(\alpha,v)\land
1\not\in\overline{\operatorname{\mathsf{Reg}}}_{-}(\alpha,v))$ (1)
$\displaystyle\lor(1\not\in\overline{\operatorname{\mathsf{Reg}}}_{+}(\alpha,v)\land\operatorname{\mathsf{Reg}}_{+}(\alpha,v)\neq\emptyset)$
(2)
$\displaystyle\lor(\operatorname{\mathsf{Reg}}_{-}(\alpha,v)\neq\emptyset\land
1\not\in\overline{\operatorname{\mathsf{Reg}}}_{-}(\alpha,v))$ (3)
$\displaystyle\lor(\operatorname{\mathsf{Reg}}_{-}(\alpha,v)\neq\emptyset\land\operatorname{\mathsf{Reg}}_{+}(\alpha,v)\neq\emptyset)$
(4)

Because of (4) and the hypothesis that
$\operatorname{\mathsf{Reg}}_{+}(\alpha,v)\cup\operatorname{\mathsf{Reg}}_{-}(\alpha,v)\neq\emptyset$,
we can consider in (1), (2) and in (3) that one of the sets
$\operatorname{\mathsf{Reg}}_{+}(\alpha,v)$ and
$\operatorname{\mathsf{Reg}}_{-}(\alpha,v)$ is non-empty, while the other is
empty.

In (2), we observe that
$1\not\in\overline{\operatorname{\mathsf{Reg}}}_{+}(\alpha,v)$ is equivalent
to $\alpha(v)\neq 1\land 1\not\in{\operatorname{\mathsf{Reg}}}_{+}(\alpha,v)$,
which together with $\operatorname{\mathsf{Reg}}_{+}(\alpha,v)\neq\emptyset$
is equivalent with $\alpha(v)\neq
1\land\operatorname{\mathsf{Reg}}_{+}(\alpha,v)=\\{0\\}$. As we noted above,
we can add to the conjunction also the term
$\operatorname{\mathsf{Reg}}_{-}(\alpha,v)=\emptyset$ since
$\operatorname{\mathsf{Reg}}_{+}(\alpha,v)\neq\emptyset$. This gives us the
clause $\alpha(v)\neq
1\land\operatorname{\mathsf{Reg}}_{+}(\alpha,v)=\\{0\\}\land\operatorname{\mathsf{Reg}}_{-}(\alpha,v)=\emptyset$.
Under the hypothesis of the theorem that
$\operatorname{\mathsf{Reg}}_{+}(\alpha,v)\cup\operatorname{\mathsf{Reg}}_{-}(\alpha,v)\neq\emptyset$,
this clause is equivalent with

$\displaystyle\overline{\operatorname{\mathsf{Reg}}}_{+}(\alpha,v)=\\{0\\}\land\operatorname{\mathsf{Reg}}_{-}(\alpha,v)=\emptyset.$
(2’)

Using a symmetric argument, we can conclude that (3) can be replaced in the
disjunction with

$\displaystyle\overline{\operatorname{\mathsf{Reg}}}_{-}(\alpha,v)=\\{0\\}\land\operatorname{\mathsf{Reg}}_{+}(\alpha,v)=\emptyset.$
(3’)

Clause (1) is equivalent with $\alpha(v)=0\land
1\not\in\operatorname{\mathsf{Reg}}_{+}(\alpha,v)\land
1\not\in\operatorname{\mathsf{Reg}}_{-}(\alpha,v)$. Since
$\operatorname{\mathsf{Reg}}_{+}(\alpha,v)$,
$\operatorname{\mathsf{Reg}}_{-}(\alpha,v)\subseteq\\{0,1\\}$ and, as noted
above, we can assume that one is empty, while the other is not, it follows
that clause (1) can be replaced with

$\displaystyle\left(\alpha(v)=0\land\operatorname{\mathsf{Reg}}_{+}(\alpha,v)=\\{0\\}\land\operatorname{\mathsf{Reg}}_{-}(\alpha,v)=\emptyset\right)\lor$
$\displaystyle\left(\alpha(v)=0\land\operatorname{\mathsf{Reg}}_{-}(\alpha,v)=\\{0\\}\land\operatorname{\mathsf{Reg}}_{+}(\alpha,v)=\emptyset\right).$

This is equivalent with

$\left(\alpha(v)=0\land\overline{\operatorname{\mathsf{Reg}}}_{+}(\alpha,v)=\\{0\\}\land\operatorname{\mathsf{Reg}}_{-}(\alpha,v)=\emptyset\right)\lor$
---
$\left(\alpha(v)=0\land\overline{\operatorname{\mathsf{Reg}}}_{-}(\alpha,v)=\\{0\\}\land\operatorname{\mathsf{Reg}}_{+}(\alpha,v)=\emptyset\right).$
(1’)

These two conditions are absorbed within the disjunction under the clauses
(2’) and (3’). This proves the claim of the lemma.

The following result is a simple consequence of Definition 2.2 and Theorem 2.6
and it offers more insight into the dynamics of strong regulatory graphs.

###### Lemma 2.8

Let $G=(V,E)$ be a strong regulatory graph, $v\in V$ and $\alpha$, $\beta$ two
states of $G$. If $\alpha\rightarrow_{G}\beta$, then

  * •

$\beta(v)=\alpha(v)$, if $\alpha(v)=1$ and
$\operatorname{\mathsf{Reg}}_{-}(\alpha,v)=\emptyset$, or $\alpha(v)=-1$ and
$\operatorname{\mathsf{Reg}}_{+}(\alpha,v)=\emptyset$, or $\alpha=0$ and
$1\not\in\operatorname{\mathsf{Reg}}_{+}(\alpha,v)\cup\operatorname{\mathsf{Reg}}_{-}(\alpha,v)$;

  * •

$\beta(v)=1$, if $1\in\operatorname{\mathsf{Reg}}+-(\alpha,v)$ and
$\operatorname{\mathsf{Reg}}_{-}(\alpha,v)=\emptyset$;

  * •

$\beta(v)=-1$, if $1\in\operatorname{\mathsf{Reg}}_{-}(\alpha,v)$ and
$\operatorname{\mathsf{Reg}}_{+}(\alpha,v)=\emptyset$;

  * •

$\beta(v)=0$, if $\operatorname{\mathsf{Reg}}_{+}(\alpha,v)\neq\emptyset$ and
$\operatorname{\mathsf{Reg}}_{-}(\alpha,v)\neq\emptyset$.

##  3 Phenotype attractors

We are interested in attractors that are defined through a fixed configuration
(also called phenotype in [20]) on some (possibly not all) of their variables.
In other words, we are interested in minimal cycles of the state transition
systems, where some of the variables of the models are constant. This is
similar to the concept of a target set in the partial controlability of
complex networks [21]. The notion of phenotype was introduced for Boolean
networks in [20] and discussed in connection with the control of Boolean
networks in [22].

###### Definition 3.1

Let $G=(V,E)$ be a strong regulatory graph. We define a _trap set_ of its
state transition graph in the usual way as a set $S$ of vertices such that for
any $\alpha\rightarrow\beta$, if $\alpha\in S$, then $\beta\in S$. An
_attractor_ is defined as (non-empty) minimal trap set under set inclusion.

A target $T$ is a set of vertices $T\subseteq V$ and a _$T$ -phenotype_ is a
function $\alpha_{T}:T\rightarrow\\{-1,1\\}$ that gives an active/inactive
status assignment of the vertices in $T$. An _$\alpha_{T}$ -phenotype
attractor_ ${\cal A}_{\alpha_{T}}$ is an attractor of the state transition
graph that has the phenotype $\alpha_{T}$ on $T$, i.e., for any $\beta\in{\cal
A}_{\alpha_{T}}$, $\beta(t)=\alpha_{T}(t)$, for all $t\in T$.

The _phenotype problem_ is to decide whether for a given phenotype $\alpha$,
there exists an $\alpha$-phenotype attractor. We prove that the problem can be
solved in polynomial time. Even more, we give a simple characterization of
phenotype attractors. We say that $u$ is an _inhibition-predecessor_
(_activation-predecessor_) of $v$ if $(u,v)\in E_{-}$ ($(u,v)\in E_{+}$,
resp.) We say that $u$ is an _inhibition-ancestor_ (_activation-ancestor_) of
$v$ if there are vertices $x_{1},x_{2},\dots,x_{n}$, $n\geq 1$ with
$(x_{i},x_{i+1})\in E_{-}$ ($(x_{i},x_{i+1})\in E_{+}$, resp.) for all $1\leq
i\leq n-1$, $x_{1}=u$ and $x_{n}=v$.

###### Theorem 3.2

Let $T\subseteq V$ be a target set and a phenotype
$\alpha_{T}:T\rightarrow\\{-1,1\\}$ on $T$. There exists an
$\alpha_{T}$-attractor if and only if the following two conditions hold:

  1. (a)

for all $v\in T$ with $\alpha_{T}(v)=1$ and for all $u\in T$ inhibition-
predecessors of $v$, we have $\alpha_{T}(u)=-1$;

  2. (b)

for all $v\in T$ with $\alpha_{T}(v)=-1$ and for all $u\in T$ activation-
ancestors of $v$, we have $\alpha_{T}(u)=-1$.

###### Proof 3.3

Assume first that there is an $\alpha_{T}$-attractor ${\cal A_{\alpha_{T}}}$.
This means that for all $\alpha\in{\cal A}_{\alpha_{T}}$ and for all $t\in T$,
$\alpha(t)=\alpha_{T}(t)$. Consider $v\in T$. If $\alpha_{T}(v)=1$, then for
all $\alpha\in{\cal A}_{\alpha_{T}}$,
$\operatorname{\mathsf{Reg}}_{-}(\alpha,v)=\emptyset$, i.e., its inhibition
predecessors $v\in T$ have $\alpha_{T}(v)=-1$. If $\alpha_{T}(v)=-1$, then for
all $\alpha\in{\cal A}_{\alpha_{T}}$,
$\operatorname{\mathsf{Reg}}_{+}(\alpha,v)=\emptyset$, i.e., its activation
predecessors $v$ have $\alpha_{T}(v)=-1$. This argument can be iterated
throughout all activation predecessors, yielding the second claim.

Consider now a phenotype $\alpha_{T}$ satisfying properties (a) and (b). We
can construct an $\alpha_{T}$-attractor ${\cal A}_{\alpha_{T}}$ in the
following way. We visit the graph going against the edges, starting from the
target set $T$, and marking them with -1 or 1 as we visit them. Let $S$ be the
set of nodes still to be explored. We start the exploration with $S=T$, whose
marking is already set through $\alpha_{T}$. Take all vertices in $S$ marked
with $1$, consider all their inhibition predecessors, mark them with -1 and
add them to $S$. Because of property (a), none of them was marked with $1$ and
so, this marking will not contradict any previously set marking. Take all
vertices in $S$ marked with $-1$, consider all their activation predecessors
and mark them with -1. Because of property (b), none of them was marked with
$1$ and so, this marking will not contradict any previously set marking.
Iterate through this step until set $S$ stops growing. To construct an
$\alpha_{T}$-attractor, consider a state $\alpha_{1}:V\rightarrow\\{-1,0,1\\}$
defined on $S$ through the markings set above, and taking an arbitrary choice
from $\\{-1,0,1\\}$ for the vertices in $V\setminus S$. We then consider the
transitions
$\alpha_{1}\rightarrow\alpha_{2}\rightarrow\alpha_{3}\rightarrow\cdots$.
Obviously, there will eventually be a repetition of states
$\alpha_{i}=\alpha_{j}$, $i<j$ in this sequence. The $\alpha_{T}$-attractor is
$\\{\alpha_{i},\alpha_{i+1},\dots,\alpha_{j-1}\\}$.

The following is a reformulation of Theorem 3.2. A vertex $t$ is an active
(inactive) target if $t\in T$ and $\alpha_{T}(t)=1$ ($\alpha_{T}(t)=-1$,
resp.) An activation path is a path in a strong regulatory graph made of
activation edges only.

###### Theorem 3.4

Let $T\subseteq V$ be a target set. For phenotype
$\alpha_{T}:T\rightarrow\\{-1,1\\}$ there is an $\alpha_{T}$-attractor if and
only if

  1. (a)

there is no activation path from an active target to an inactive target and

  2. (b)

there is no activation path plus a final inhibition edge between two active
targets.

##  4 Application to a regulatory cancer network

We discuss in this section a strong regulatory graph model of some of the key
elements of the RTK (receptor tyrosine kinase) signaling through the MAPK-
((mitogen - activated protein kinase) and the PI3K/AKT- (lipid kinase
phoshoinositide-3-kinase) pathways. Our goal is to demonstrate that the strong
regulatory graphs are expressive enough to capture some interesting properties
of a biological model. The MAPK signaling pathway communicates signals from
outside the cell to the nucleus of the cell, it is involved in cell growth and
proliferation, it is often mutated in cancer and a popular target of cancer
treatments [23, 24, 25]. It interplays with the PI3K/AKT signaling pathway, a
master regulator of the cell, whose activation contributes to the development
of tumors and resistance to anticancer therapies [26, 27]. We consider here a
portion of the Boolean network model of [14]. The model is illustrated in
Figure 3, with the key regulators shown as vertices, with green pointed arrows
for the activation edges and red blunt arrows for the inhibition. The model
includes FOXO3, a protein known for its role in inducing cell death, often
inhibited in tumors, and AKT, with a role in cell survival and often activated
in tumors. Following [14], the activation of RTK depends on the presence of
growth factors, not included in the model, and so its status is set to be
constant $-1$. We consider the $(-1,1)$ status of (FOXO3, AKT) as an indicator
of uncontrolled proliferation, that of $(1,-1)$ an indicator of non-
proliferation, and that of $(-1,-1)$ an indicator of moderate proliferation.
(Because of the inactivation edge from AKT to FOXO3, an $(1,1)$-attractor is
not possible.) These are of course over-simplifications of a much more complex
interplay of interactions. We show that the strong regulatory graph model has
the same variety of outcomes as the Boolean network of [15], while adopting
the strong update rule proposed in this paper and the ambiguity status in case
of conflicting regulation.

Figure 3: The MAPK-PI3K/AKT signaling pathways.

In our discussion of the dynamics of the model, we write the variables in the
order (RTK, RAS, PI3K, MAPK, PIP3, FOXO3, AKT). The model has several
phenotype attractors associated with all three possibilities of the
(FOXO3,AKT) phenotype, there are $(-1,-1)$-, $(1,-1)$-, and
$(-1,1)$-attractors:

  * •

for $s_{1}=(-1,-1,-1,1,-1,{\bf-1},{\bf-1})$, $s_{1}\rightarrow s_{1}$;

  * •

for $s_{2}=(-1,-1,-1,-1,-1,{\bf 1},{\bf-1})$, $s_{2}\rightarrow s_{2}$;

  * •

for $s_{3}=(-1,1,1,1,1,{\bf-1},{\bf 1})$, $s_{3}\rightarrow s_{3}$.

A mutation on $PI3K$, making it constant active, gives proliferation as the
only stable set, consistent with the biological observations in [26] and with
the behavior of the Boolean network model of [15]: for $s_{4}=(-1,-1,{\bf
1},1,1,{\bf-1},{\bf 1})$, $s_{4}\rightarrow s_{4}$. Also consistent with the
Boolean network model of [15] is the observation that setting PI3K to inactive
(corresponding to the idea of targeting it with drug inhibitors), allows all
proliferation outcomes to be possible:

  * •

$(-1,-1,{\bf-1},-1,1,{\bf 1},{\bf-1})\rightarrow(-1,-1,{\bf-1},1,1,{\bf
1},{\bf 1})\rightarrow(-1,-1,{\bf-1},1,1,{\bf-1},{\bf 1})\rightarrow(-1,$
$-1,{\bf-1},1,1,{\bf-1},{\bf 1}$, with the final vertex begin a singleton
attractor;

  * •

$(-1,-1,{\bf-1},-1,-1,{\bf 1},{\bf-1})$ is a singleton attractor;

  * •

$(-1,-1,{\bf-1},1,-1,{\bf
1},{\bf-1})\rightarrow(-1,-1,{\bf-1},1,-1,{\bf-1},{\bf-1})$, with the final
vertex begin a singleton attractor.

##  5 Discussion

Our main motivation in this article was to offer a solution to the gap in the
current modeling of large network models. Current network models span between
structural modeling approaches (well-supported by interaction data, scalable,
but not dynamical) and detailed logical models (dynamical, but problematic to
scale up because the data is typically not detailed enough to specify the
update functions of each variable). We introduced the concept of strong
regulatory graph, a form of logical network model where the active/inactive
status of a vertex is set only if its regulators agree in their influences;
otherwise, the status of the vertex is set to ambiguous. Ambiguity here means
that the node may potentially be active or inactive, depending on the exact
numerical setup of its predecessors, a level of detail that is left out of
logical models. The ambiguity may have a cascading influence over the
successors of the node, but it may also cancel out through the clarifications
brought by other active vertices in the network. We discussed an update rule
that defines the interplay between active, inactive, and ambiguous regulators.
We also discussed the phenotype attractor problem and showed that it is easy
to decide whether an attractor of a given phenotype exists in a strong
regulatory graph.

It is straightforward to see that for any strong regulatory graph there is a
Boolean network that can simulate it. The state transition rule is indeed
based on simple logical tests, and so easily implementable through Boolean
functions. Also, the range of values $\\{-1,0,1\\}$ can easily be implemented
using two Boolean variables for each vertex of the strong regulatory graph.
This observation connects the strong regulatory graphs to the rich literature
on the control of Boolean networks.

We discussed only the synchronous version of the strong regulatory graphs,
where all vertices are updated simultaneously. Conceptually, the asynchronous
version is similar, but it leads to a non-deterministic state-transition
systems, whose properties should be interesting to characterize in a further
study. Also left for another study is the study of controlability of strong
regulatory graphs. Its targeted version can be easily defined in terms of
reaching a phenotype attractor from any state with some minimal interventions
in the state transition system.

## References

  * [1] Bloomingdale P, Nguyen VA, Niu J, Mager DE.  Boolean network modeling in systems pharmacology.  _Journal of Pharmacokinetics and Pharmacodynamics_ , 2018.  45(1):159–180.  10.1007/s10928-017-9567-4.  URL https://doi.org/10.1007/s10928-017-9567-4.
  * [2] Cacace E, Collombet S, Thieffry D.  Chapter Seven - Logical modeling of cell fate specification—Application to T cell commitment.  In: Peter IS (ed.), Gene Regulatory Networks, volume 139 of _Current Topics in Developmental Biology_ , pp. 205–238. Academic Press, 2020\.  https://doi.org/10.1016/bs.ctdb.2020.02.008.  URL https://www.sciencedirect.com/science/article/pii/S007021532030034X.
  * [3] Montagud A, Béal J, Tobalina L, Traynard P, Subramanian V, Szalai B, Alföldi R, Puskás L, Valencia A, Barillot E, Saez-Rodriguez J, Calzone L.  Patient-specific Boolean models of signalling networks guide personalised treatments.  _eLife_ , 2022.  11:e72626.  10.7554/eLife.72626.  URL https://doi.org/10.7554/eLife.72626.
  * [4] Sizek H, Hamel A, Deritei D, Campbell S, Ravasz Regan E.  Boolean model of growth signaling, cell cycle and apoptosis predicts the molecular mechanism of aberrant cell cycle progression driven by hyperactive PI3K.  _PLOS Computational Biology_ , 2019.  15(3):1–43.  10.1371/journal.pcbi.1006402.  URL https://doi.org/10.1371/journal.pcbi.1006402.
  * [5] Checcoli A, Pol JG, Naldi A, Noel V, Barillot E, Kroemer G, Thieffry D, Calzone L, Stoll G.  Dynamical Boolean Modeling of Immunogenic Cell Death.  _Frontiers in Physiology_ , 2020.  11.  10.3389/fphys.2020.590479.  URL https://www.frontiersin.org/articles/10.3389/fphys.2020.590479.
  * [6] Wang RS, Saadatpour A, Albert R.  Boolean modeling in systems biology: an overview of methodology and applications.  _Physical Biology_ , 2012.  9(5):055001.  10.1088/1478-3975/9/5/055001.  URL https://dx.doi.org/10.1088/1478-3975/9/5/055001.
  * [7] Kanehisa M, Sato Y, Kawashima M, Furumichi M, Tanabe M.  KEGG as a reference resource for gene and protein annotation.  _Nucleic Acids Research_ , 2015.  44(D1):D457–D462.  10.1093/nar/gkv1070.  https://academic.oup.com/nar/article-pdf/44/D1/D457/9482226/gkv1070.pdf, URL https://doi.org/10.1093/nar/gkv1070.
  * [8] Türei D, Korcsmáros T, Saez-Rodriguez J.  OmniPath: guidelines and gateway for literature-curated signaling pathway resources.  _Nature Methods_ , 2016.  13:966–967.  https://doi.org/10.1038/nmeth.4077.
  * [9] Breuer K, Foroushani AK, Laird MR, Chen C, Sribnaia A, Lo R, Winsor GL, Hancock REW, Brinkman FSL, Lynn DJ.  InnateDB: systems biology of innate immunity and beyond–recent updates and continuing curation.  _Nucleic Acids Research_ , 2013.  41(Database issue):D1228–D1233.  https://doi.org/10.1093/nar/gks1147.
  * [10] Licata L, Lo Surdo P, Iannuccelli M, Palma A, Micarelli E, Perfetto L, Peluso D, Calderone A, Castagnoli L, Cesareni G.  SIGNOR 2.0, the SIGnaling Network Open Resource 2.0: 2019 update.  _Nucleic Acids Research_ , 2020.  48(D1):D504–D510.  https://doi.org/10.1093/nar/gkz949.
  * [11] Wishart DS, Feunang YD, Guo AC, Lo EJ, Marcu A, Grant JR, Sajed T, Johnson D, Li C, Sayeeda Z, Assempour N, Iynkkaran I, Liu Y, Maciejewski A, Gale N, Wilson A, Chin L, Cummings R, Le D, Pon A, Knox C, Wilson M.  DrugBank 5.0: a major update to the DrugBank database for 2018.  _Nucleic Acids Research_ , 2017.  46(D1):D1074–D1082.  https://doi.org/10.1093/nar/gkx1037.
  * [12] Siminea N, Popescu V, Sanchez Martin JA, Florea D, Gavril G, Gheorghe AM, Iţcuş C, Kanhaiya K, Pacioglu O, Popa LI, Trandafir R, Tusa MI, Sidoroff M, Păun M, Czeizler E, Păun A, Petre I.  Network analytics for drug repurposing in COVID-19.  _Briefings in Bioinformatics_ , 2021.  23(1).  10.1093/bib/bbab490.  Bbab490, https://academic.oup.com/bib/article-pdf/23/1/bbab490/42508594/bbab490.pdf, URL https://doi.org/10.1093/bib/bbab490.
  * [13] Popescu VB, Kanhaiya K, Năstac DI, Czeizler E, Petre I.  Network controllability solutions for computational drug repurposing using genetic algorithms.  _Scientific Reports_ , 2022.  12(1):1437.  10.1038/s41598-022-05335-3.  URL https://doi.org/10.1038/s41598-022-05335-3.
  * [14] GT Zañudo J, Steinway SN, Albert R.  Discrete dynamic network modeling of oncogenic signaling: Mechanistic insights for personalized treatment of cancer.  _Current Opinion in Systems Biology_ , 2018.  9:1–10.  https://doi.org/10.1016/j.coisb.2018.02.002.  Mathematic modelling, URL https://www.sciencedirect.com/science/article/pii/S2452310017302032.
  * [15] Zanudo JG, Aldana M, Martínez-Mekler G.  Boolean threshold networks: Virtues and limitations for biological modeling.  In: Information Processing and Biological Systems, pp. 113–151. Springer, 2011.
  * [16] Naldi A, Thieffry D, Chaouiya C.  Decision Diagrams for the Representation and Analysis of Logical Models of Genetic Networks.  In: Calder M, Gilmore S (eds.), Computational Methods in Systems Biology. Springer Berlin Heidelberg, Berlin, Heidelberg.  ISBN 978-3-540-75140-3, 2007 pp. 233–247.
  * [17] Thomas R.  Regulatory networks seen as asynchronous automata: a logical description.  _Journal of Theoretical Biology_ , 1991.  153(1):1–23.
  * [18] Thomas R, Thieffry D, Kaufman M.  Dynamical behaviour of biological regulatory networks—I. Biological role of feedback loops and practical use of the concept of the loop-characteristic state.  _Bulletin of Mathematical Biology_ , 1995.  57(2):247–276.
  * [19] Chaouiya C, Remy E, Mossé B, Thieffry D.  Qualitative Analysis of Regulatory Graphs: A Computational Tool Based on a Discrete Formal Framework.  In: Benvenuti L, De Santis A, Farina L (eds.), Positive Systems. Springer Berlin Heidelberg, Berlin, Heidelberg.  ISBN 978-3-540-44928-7, 2003 pp. 119–126.
  * [20] Klarner H, Heinitz F, Nee S, Siebert H.  Basins of Attraction, Commitment Sets, and Phenotypes of Boolean Networks.  _IEEE/ACM Trans. Comput. Biol. Bioinformatics_ , 2020.  17(4):1115–1124.  10.1109/TCBB.2018.2879097.  URL https://doi.org/10.1109/TCBB.2018.2879097.
  * [21] Kanhaiya K, Czeizler E, Gratie C, Petre I.  Controlling Directed Protein Interaction Networks in Cancer.  _Scientific Reports_ , 2017.  7(1):10327.  10.1038/s41598-017-10491-y.  URL https://10.1038/s41598-017-10491-y.
  * [22] Cifuentes-Fontanals L, Tonello E, Siebert H.  Control in Boolean Networks With Model Checking.  _Frontiers in Applied Mathematics and Statistics_ , 2022.  8.  10.3389/fams.2022.838546.  URL https://www.frontiersin.org/articles/10.3389/fams.2022.838546.
  * [23] Yaeger R, Corcoran RB.  Targeting Alterations in the RAF–MEK Pathway.  _Cancer Discovery_ , 2019.  9(3):329–341.  10.1158/2159-8290.CD-18-1321.  https://aacrjournals.org/cancerdiscovery/article-pdf/9/3/329/2942719/329.pdf, URL https://doi.org/10.1158/2159-8290.CD-18-1321.
  * [24] Burotto M, Chiou VL, Lee JM, Kohn EC.  The MAPK pathway across different malignancies: A new perspective.  _Cancer_ , 2014.  120(22):3446–3456.  https://doi.org/10.1002/cncr.28864.  https://acsjournals.onlinelibrary.wiley.com/doi/pdf/10.1002/cncr.28864, URL https://acsjournals.onlinelibrary.wiley.com/doi/abs/10.1002/cncr.28864.
  * [25] Dhillon A, Hagan S, Rath O, Kolch W.  MAP kinase signalling pathways in cancer.  _Oncogene_ , 2007.  26:3279–3290.  https://doi.org/10.1038/sj.onc.1210421.
  * [26] Yang J, Nie J, Ma X, Wei Y, Peng Y, Wei X.  Targeting PI3K in cancer: mechanisms and advances in clinical trials.  _Mol Cancer_ , 2019.  18(26).  https://doi.org/10.1186/s12943-019-0954-x.
  * [27] He Y, Sun MM, Zhang GG, Yang J, Chen KS, Xu WW, Li B.  TTargeting PI3K/Akt signal transduction for cancer therapy.  _Sig Transduct Target Ther_ , 2021.  6(425).  https://doi.org/10.1038/s41392-021-00828-5.