11institutetext: Dipartimento di Sistemi, Informatica, Applicazioni,
Università degli Studi di Firenze, Firenze, Italy 11email:
<EMAIL_ADDRESS>Departamento de Matemática,
Universidade Federal do Ceará, Fortaleza, CE, Brazil 22email:
<EMAIL_ADDRESS>
# Königsberg Sightseeing: Eulerian Walks in Temporal Graphs

Andrea Marino 11 Ana Silva 22

###### Abstract

An Eulerian walk (or Eulerian trail) is a walk (resp. trail) that visits every
edge of a graph $G$ at least (resp. exactly) once. This notion was first
discussed by Leonhard Euler while solving the famous Seven Bridges of
Königsberg problem in 1736. What if Euler had to take a bus? In a temporal
graph $(G,\lambda)$, with $\lambda:E(G)\to 2^{[\tau]}$, an edge $e\in E(G)$ is
available only at the times specified by $\lambda(e)\subseteq[\tau]$, in the
same way the connections of the public transportation network of a city or of
sightseeing tours are available only at scheduled times. In this scenario,
even though several translations of Eulerian trails and walks are possible in
temporal terms, only a very particular variation has been exploited in the
literature, specifically for infinite dynamic networks (Orlin, 1984). In this
paper, we deal with temporal walks, local trails, and trails, respectively
referring to edge traversal with no constraints, constrained to not repeating
the same edge in a single timestamp, and constrained to never repeating the
same edge throughout the entire traversal. We show that, if the edges are
always available, then deciding whether $(G,\lambda)$ has a temporal walk or
trail is polynomial, while deciding whether it has a local trail is
$\NP$-complete even if it has lifetime 2. In contrast, in the general case,
solving any of these problems is $\NP$-complete, even under very strict
hypothesis.

##  1 Introduction

An Eulerian walk (or Eulerian trail) is a walk (resp. trail) that visits every
edge of a graph $G$ at least (resp. exactly) once. The Eulerian trail notion
was first discussed by Leonhard Euler while solving the famous Seven Bridges
of Königsberg problem in 1736, where one wanted to pass by all the bridges
over the river Preger without going twice over the same bridge. Imagine now a
similar problem, where you have a set of sights linked by possible routes. If
the routes themselves are also of interest, a sightseeing tourism company
might want to plan visits on different days that cover all the routes. One
could do that with no constraints at all (thus performing a walk), or with the
very strict constraint of never repeating a route (thus getting a trail), or
constraining oneself to at least not repeting the same route on the same day
(thus getting what we called a local trail). If we further assume that some
routes might not be always accessible, we then get distinct problems defined
on temporal graphs.

In a temporal graph $(G,\lambda)$, with $\lambda:E(G)\to 2^{[\tau]}$, an edge
$e\in E(G)$ is available only at the times specified by
$\lambda(e)\subseteq[\tau]$, in the same way the connections of the public
transportation network of a city or of sightseeing tours are available only at
scheduled times. In this scenario, paths and walks are valid only if they
traverse a sequence of adjacent edges $e_{1},\ldots,e_{k}$ such that for each
$i\in[k-1]$, $\lambda(e_{i})\leq\lambda(e_{i+1})$, i.e. whose time sequence is
non-decreasing (similarly, strictly increasing sequences, i.e. with
$\lambda(e_{i})<\lambda(e_{i+1})$, can be considered).

Several translations of Eulerian trails and walks are possible in temporal
terms, depending on the constraints we consider. In particular, we study the
following variations. Below, all the walks and trails are implicitly
considered to be temporal, as defined in the previous paragraph.

######  Problem 1.

Given a temporal graph $(G,\lambda)$, we consider the following problems:

  * •

Eulerian Walk: deciding whether $(G,\lambda)$ has an Eulerian walk, i.e. a
walk traversing each edge at least once.

  * •

Eulerian Local Trail: deciding whether $(G,\lambda)$ has an Eulerian local
trail, i.e. a walk traversing each edge at least once and at most once in each
time stamp.

  * •

Eulerian Trail: deciding whether $(G,\lambda)$ has an Eulerian trail, i.e. a
walk traversing each edge exactly once.

We also consider the related problems where the walks/trails are closed (first
vertex equal to the last one), respectively referring to them as Eulerian
Closed Walk, Eulerian Local Tour, and Eulerian Tour.

The research on temporal graphs have attracted a lot of attention in the past
decade (we refer the reader to the surveys [18, 16] and the seminal paper
[15]), and they appear also under different names, e.g. as time-varying graphs
[3], as evolving networks [2], and as link streams [16]. Even after all the
received attention, surprisingly enough none of the above problems have been
previously considered.

Concerning Eulerian Walk, one of the closest concepts is the one defined by
Orlin in [20], where he gives a polynomial algorithm to check the existence of
an Eulerian closed walk (i.e. a _tour_) in dynamic graphs. However, the
dynamic graph model is quite different from the temporal graph model used in
this paper, as pointed out in [18]. Indeed, looking at the corresponding time-
expanded graph, temporal edges can go back in time and the graph is infinite.
Nevertheless, the results presented there seemed to point towards the
polinomiality of the problems investigated here, as observed in [18]: “the
results proved for it [the dynamic graph model] are resounding and possibly
give some first indications of what to expect when adding to combinatorial
optimization problems a time dimension”. We found however that this is not the
case, as we will show that even Eulerian Walk turns out to be much harder on
temporal graphs (for more details about dynamic graphs and the result proved
by Orlin in [20] see Appendix 0.A). Taking inspiration in [20], we also define
a dynamic-based temporal graph as a temporal graph whose edges are always
available, and we analyze the complexity of the above problems on these
particular instances.

Still concerning Eulerian Walk, a closely related problem is the Temporal
Exploration problem (TEXP) [19], which consists of, given a temporal graph
$(G,\lambda)$, finding a temporal walk that visits all vertices in $G$
(possibly, more than once) whose arrival time is minimum. In [19], they prove
that this problem is -complete and even not approximable unless $\P=\NP$; this
is in stark contrast with the static version of the problem, which can be
trivially solved in linear time. However this is quite different wrt our
Eulerian Walk, as we are considering walks passing through all the edges
instead of all the vertices, and this difference is crucial as a
transformation from TEXP to Eulerian Walk does not seem to be easy. Indeed, by
simply transforming each vertex into an edge, we would get two types of edges,
connection-edges and vertex-edges, and a negative answer to Eulerian Walk
could be just motivated by the impossibility of visiting all the connection-
edges, which are actually not required to be visited to solve TEXP. This is
not surprising as we are trying to transform a Hamiltonian walk (i.e. a walk
passing by each vertex) into an Eulerian walk. We mention that, unlike the
problems investigated here, a lot of research has been devoted to temporal
node exploration, e.g. bounding the arrival time of such walks in special
instances [9, 8] and extending previous results in the case of non-strictly
increasing paths [10].

Concerning Eulerian Trail and Eulerian Local Trail, observe that, when
$\tau=1$, then both of them degenerate into the original formulation of the
Seven Bridges of Königsberg problem. This is why we think they appear to be
more natural adaptations of the static version of the problem. Nevertheless
they have never been investigated before, up to our knowledge.

#### Our results

Our results are summarized in Table 1 and detailed in Theorem 1.1.

######  Theorem 1.1.

Given a temporal graph $(G,\lambda)$

  1. 1.

Eulerian Walk is -complete, even if each snapshot of $(G,\lambda)$ is a forest
of constant size, while it is is polynomial for bounded $\tau$. It is also
polynomial if $(G,\lambda)$ is dynamic-based.

  2. 2.

Eulerian Local Trail is -complete for each $\tau\geq 2$ and in the case
$(G,\lambda)$ is dynamic-based.

  3. 3.

Eulerian Trail is -complete for each $\tau\geq 2$. It is polynomial if
$(G,\lambda)$ is dynamic-based.

Same applies to tours, i.e. Eulerian Closed Walk, Eulerian Local Tour, and
Eulerian Tour.

Theorem 1.1 gives a complete taxonomy of our problems, also focusing on the
possibility of getting polynomial algorithms when we have a small lifetime
$\tau$. In particular, for Eulerian Trail and Eulerian Local Trail, since they
become polynomial when $\tau=1$, the bound for $\tau$ is optimal, giving us a
complete dichotomy with respect to the lifetime of $(G,\lambda)$. In contrast
Eulerian Walk is easily solvable for every fixed $\tau$, showing that walks
are easier than trails even on the temporal context.

It is important to remark that none of the above variations immediately
implies any of the others. We will show indeed that the property of being
Eulerian for the static base graph $G$ is in general a necessary but not
sufficient condition for the existence of an Eulerian trail, becoming
sufficient only if we restrict to dynamic-based temporal graphs. In the case
of Eulerian local trail, we will see that this property is not even necessary.

Finally, as a by product of our reductions we get the following result about
static graphs, which can be of independent interest.

######  Corollary 1.

Given a graph $G$, deciding whether all the edges of $G$ can be covered with
two trails is -complete.

| Eulerian Walk | Eulerian Local Trail | Eulerian Trail
---|---|---|---
General |  -c for $\tau$ unbounded |  -c [from below] |  -complete
Poly for $\tau$ fixed |  | for $\tau=2$
Dynamic-Based | Poly† |  -c for $\tau=2$ | Poly⋆
Table 1: Our results concerning Problem 1. For general temporal graphs (first
row) and for dynamic-based temporal graphs (second row). † corresponds to
deciding whether $G$ is connected. ⋆ corresponds to deciding whether $G$ has
an Eulerian trail.

#### Further Related Work.

Other than the papers on TEXP previously mentioned, there is a vast literature
about finding special paths or walks in temporal graphs, and some interesting
papers include [4, 24, 21, 19]. Also, it is largely known that a static graph
$G$ has an Eulerian tour (trail) if and only if $G$ has at most one non-
trivial component and all the vertices have even degree (at most two vertices
have odd degree). A graph is called Eulerian if it has an Eulerian tour.

When considering dynamic-based temporal graphs, as edges are assumed to be
always available during the lifetime $\tau$, we could relate our problems to
several other problems on static graphs. A closely related one would be the
Chinese Postman problem, where the edges of the graph have positive weights
and one wants to find an Eulerian closed walk on $G$ with minimum weight; in
other words, one wants to add copies of existing edges in order to obtain an
Eulerian graph of minimum sum weight. Even if we regard the Chinese Postman
problem where the weights are all equal to 1, this is very different from our
approach since for us, repetition of a long common trail in different
snapshots does not make the solution worse, while it would considering the
Chinese Postman problem. It is easy to see though that the solution for the
Chinese Postman would give us an upper bound for the amount of time spent on
an Eulerian local tour of a dynamic-based graph, as we could start a new trail
on a new snapshot whenever an edge repetition was detected. The Chinese
Postman problem is largely known to be polynomial [13], and some variations
that take time into consideration have been investigated, mostly from the
practical point of view (see e.g. [5, 22, 23]), but none of which is
equivalent to our problem.

The problem of trying to obtain an Eulerian subgraph (as opposed to a
supergraph) has also been studied. In [6], the authors study a family of
problems where the goal is to make a static graph Eulerian by a minimum number
of deletions. They completely classify the parameterized complexity of various
versions of the problem: vertex or edge deletions, undirected or directed
graphs, with or without the requirement of connectivity. Also in [11], the
parameterized complexity of the following Euler subgraph problems is studied:
(i) Largest Euler Subgraph: for a given graph $G$ and integer parameter $k$,
does $G$ contain an induced Eulerian subgraph with at least $k$ vertices?; and
(ii) Longest Circuit: for a given graph $G$ and integer parameter $k$, does
$G$ contain an Eulerian subgraph with at least $k$ edges?

Eulerian Local Trail on dynamic-based graphs is actually more closely related
to the problem of covering the edges of a graph with the minimum number of
(not necessarily disjoint) trails, whereas the aforementioned problems are
more concerned with either minimizing edge repetitions or maximizing the
subgraph covered by a single trail. Even if the trail cover problem can be so
naturally defined and involve such a basic structure as trail, up to our
knowledge it has not yet been previously investigated. Note that Eulerian
Local Trail is slightly different from trail cover, since we also require that
together the trails form a walk. In any case, a small modification of our
proof of Theorem 1.1.2 implies that deciding whether the edges of a graph can
be covered with at most two trails is $\NP$-complete (Corollary 1).
Interestingly enough, the vertex version of this problem, namely the path
cover problem, has been largely investigated (see e.g. [1, 14, 17]).

#### Preliminaries.

We use and extend the notation in [18]. A temporal graph is a graph together
with a function on the edges saying when each edge is active; more formally, a
_temporal graph_ is a pair $(G,\lambda)$, where $\lambda:E(G)\rightarrow
2^{\mathbb{N}-\\{0\\}}$. Here, we consider only finite temporal graphs, i.e.,
graphs such that $\max\bigcup_{e\in E(G)}\lambda(e)$ is defined. This value is
called the _lifetime of $(G,\lambda)$_ and denoted by $\tau$. Given
$i\in[\tau]$, we define the _snapshot $G_{i}$_ as being the subgraph of $G$
containing exactly the edges active in time $i$; more formally,
$V(G_{i})=V(G)$ and $E(G_{i})=\\{e\in E(G)\mid i\in\lambda(e)\\}$.

Given vertices $v_{0},v_{k}$ in a graph $G$, a $v_{0},v_{k}$-walk in $G$ is an
alternating sequence $(v_{0},e_{1},v_{1},\ldots,e_{k},v_{k})$ of vertices and
edges such that $e_{i}$ goes from $v_{i-1}$ to $v_{i}$ for
$i\in\\{1,\ldots,k\\}$. We define a walk in a temporal graph similarly, except
that a walk cannot go back in time. More formally, given a temporal graph
$(G,\lambda)$ and a $v_{0},v_{k}$-walk
$W=(v_{0},e_{1},v_{1},\ldots,e_{k},v_{k})$, we say that $W$ is a temporal
$v_{0},v_{k}$-walk if
$\lambda(e_{1})\leq\lambda(e_{2})\leq\ldots\leq\lambda(e_{k})$. It is _closed_
if it starts and finishes on the same vertex of $G$, i.e., if $v_{0}=v_{k}$.

We say that a temporal walk $W$ is a _local trail_ if there are no two
occurrences of the same edge of $G$ in the same snapshot, i.e., if $W$
restricted to $G_{i}$ is a trail in $G$ for every $i\in[\tau]$. We say that
$W$ is a _trail_ if there are no two occurrences of the same edge of $G$ in
$W$. A closed (local) trail is also called a _(local) tour_. Finally, a
temporal walk $W$ is called _Eulerian_ if at least one copy of each edge of
$G$ appears at least once in $W$. Observe that, by definition, an Eulerian
trail visits every edge exactly once.

A _dynamic-based graph_ is a temporal graph $(G,\lambda)$ where the edges are
always available.111This is the reason why we use the term dynamic-based, as
they are similar to the dynamic networks used in [20] when studying Eulerian
trails, except that edges cannot go back in time and the lifetime is finite.
We denote a dynamic-based graph simply by $(G,[\tau])$ where $\tau$ is the
lifetime of the temporal graph.

##  2 Eulerian Walk

In this section we focus on Eulerian Walk, i.e. deciding if there is a
temporal walk passing by each edge at least once, proving the results in Item
1 in Theorem 1.1, summarized in the first column of Table 1.

In particular, a preliminary result, whose formal proof is reported in
Appendix 0.B.1, concerns the case where the lifetime $\tau$ is bounded. It
consists basically of checking whether there is a choice of connected
components $H_{1},\ldots,H_{\tau}$, one for each timestamp $i$, that together
cover all the edges of $G$ and is such that $H_{i}$ intersects $H_{i+1}$, for
each $i\in[\tau-1]$.

######  Lemma 1.

Given a temporal graph $(G,\lambda)$ with fixed lifetime $\tau$, solving
Eulerian Walk on $(G,\lambda)$ can be done in time $O((n+m)\cdot n^{\tau-1})$,
where $n=|V(G)|$ and $m=|E(G)|$.

In the following, we show that when $\tau$ is unbounded, deciding whether
$(G,\lambda)$ admits an Eulerian walk is $\NP$-complete by reducing from 3-.
This is best possible because of the above lemma.

######  Theorem 2.1.

Given a temporal graph $(G,\lambda)$, deciding whether $(G,\lambda)$ admits an
Eulerian walk is $\NP$-complete, even if each snapshot of $(G,\lambda)$ is a
forest of constant size.

###### Proof.

We make a reduction from 3-. Let $\phi$ be a 3-CNF formula on variables
$\\{x_{1},\cdots,x_{n}\\}$ and clauses $\\{c_{1},\cdots,c_{m}\\}$, and
construct $G$ as follows. For each clause $c_{i}$, add vertices
$\\{a_{i},b_{i}\\}$ to $G$ and edge $a_{i}b_{i}$. Now consider a variable
$x_{i}$, and let $c_{i_{1}},\cdots,c_{i_{p}}$ be the clauses containing
$x_{i}$ positively, and $c_{j_{1}},\cdots,c_{j_{q}}$ be the clauses containing
$x_{i}$ negatively. Add two new vertices $x_{i},\overline{x}_{i}$ to $G$, and
edges $\\{x_{i}a_{i_{k}}\mid k\in[p]\\}\cup\\{\overline{x}_{i}a_{j_{k}}\mid
k\in[q]\\}$; denote the spanning subgraph of $G$ formed by these edges by
$H_{i}$, and let $H^{\prime}_{i}$ be equal to $H_{i}$ together with edges
$\\{a_{i}b_{i}\mid i\in\\{i_{1},\cdots,i_{p},j_{1},\cdots,j_{q}\\}$. We can
suppose that $\\{i_{1},\cdots,i_{p}\\}\cap\\{j_{1},\cdots,j_{q}\\}=\emptyset$
as otherwise the clauses in the intersection would always be trivially valid;
thus we get that $H_{i},H^{\prime}_{i}$ are forests. Finally, add a new vertex
$T$ and make it adjacent to every vertex in $\\{x_{i},\overline{x}_{i}\mid
i\in[n]\\}$.

We now describe the snapshots of $(G,\lambda)$. See Figure 1 to follow the
construction. We first build 2 consecutive snapshots in $(G,\lambda)$ related
to $x_{i}$, for each $i\in[n]$. The first one is equal to $H^{\prime}_{i}$,
and the second one contains exactly the edges
$\\{Tx_{i},T\overline{x}_{i},Tx_{i+1},T\overline{x}_{i+1}\\}$ if $i<n$, and if
$i=n$, then the second snapshot is equal to $G-\\{a_{j}b_{j}\mid j\in[m]\\}$;
this can be done because this subgraph is connected. Denote by
$S_{i}^{1},S_{i}^{2}$ the first and second snapshot of $x_{i}$, for each
$i\in[n]$. Put these snapshots consecutively in timestamps $1$ through $2n$,
in the order of the indexing of the variables. For now, observe that only the
last snapshot might not be a forest; this will be fixed later. We now prove
that $\phi$ is a satisfiable formula if and only if $(G,\lambda)$ admits an
Eulerian walk.

$G_{1}$$G_{2}$$G_{3}$$x_{1}$$a_{1}$$a_{2}$$b_{1}$$b_{2}$$\overline{x}_{1}$$a_{3}$$a_{4}$$b_{3}$$b_{4}$$x_{1}$$\overline{x}_{1}$$x_{2}$$\overline{x}_{2}$$T$$x_{2}$$a_{1}$$a_{2}$$b_{1}$$b_{2}$$\overline{x}_{2}$$a_{3}$$a_{4}$$b_{3}$$b_{4}$
Figure 1: First three snapshots of the construction. In this example, we have
$c_{1}$ containing $(x_{1}\vee x_{2})$, $c_{2}$ containing
$(x_{1}\vee\overline{x}_{2})$, $c_{3}$ containing
$(\overline{x}_{1}\vee{x}_{2})$, and $c_{4}$ containing
$(\overline{x}_{1}\vee\overline{x}_{2})$.

Observe that, since we are dealing with a walk, we are allowed to repeat edges
as many times as we want; hence, if we visit any vertex inside a component of
a snapshot, we can also visit the entire component. At the same time, it is
not possible to visit more than one component within a snapshot. This is
implicitly used below.

First consider a satisfying assignment of $\phi$. Now, construct an Eulerian
walk as follows. Start by visiting all the edges in the component of
$H^{\prime}_{1}$ containing $x_{1}$, if $x_{1}$ is true, or the one containing
$\overline{x}_{1}$, otherwise. Then, in $S^{2}_{1}$, jump to $x_{2}$ if
$x_{2}$ is true, or to $\overline{x}_{2}$, otherwise. Repeat the process until
reaching $S^{1}_{n}$, and at the last snapshot, visit all the edges in
$G-\\{a_{j}b_{j}\mid j\in[m]\\}$. Because each clause $c_{i}$ contains at
least one true variable, we know that edge $a_{i}b_{i}$ is visited.

Now, consider an Eulerian walk $W$ of $(G,\lambda)$, and denote by $W_{i}^{j}$
the walk $W$ restricted to $S^{j}_{i}$, for each $i\in[n]$ and $j\in[2]$. We
set $x_{i}$ to true if $W_{i}^{1}$ contains $x_{i}$, and to false otherwise.
Now, consider a clause $c_{i}$ containing variables
$x_{k_{1}},x_{k_{2}},x_{k_{3}}$. Because $a_{i}b_{i}$ appears only in
snapshots $S^{1}_{k_{1}},S^{1}_{k_{2}},S^{1}_{k_{3}}$, and only in the
component containing the literal that appears in $c_{i}$, we get that at least
one of the three literals must be set to true.

Finally, observe that we could repeat the same pattern as the one in the first
$2n-1$ snapshots in order to visit the remaining edges in $G-\\{a_{j}b_{j}\mid
j\in[m]\\}$. As long as the edges in $\\{a_{i}b_{i}\mid i\in[m]\\}$ appear
only in the first $2n-1$ snapshots, the same argument as before still applies,
and we get the further constraint that each snapshot is a forest.
Additionally, they can also be considered to have constant size since 3-is
$\NP$-complete even if each variable appears at most three times [7]. ∎

Now, if we consider a dynamic-based graph $(G,\lambda)$, since all the edges
are active throughout its lifetime, we clearly have that there exists an
Eulerian walk if and only if $G$ is connected, as highlighted by the following
Lemma.

######  Lemma 2.

Eulerian Walk is polynomial for dynamic-based temporal graphs.

By Lemma 1, Theorem 2.1, and Lemma 2, we obtain Item 1 of Theorem 1.1.
Finally, note that if one is interested in closed walks instead, not only our
$\NP$-completeness reduction can be adapted in order to ensure that we can
always go back to the initial vertex, but also the complexity results still
hold.

##  3 Eulerian Local Tours and Trails

In this section we focus on Item 2 of Theorem 1.1. In the whole section, we
will focus on dynamic-based temporal graphs as the hardness results for
general temporal graphs are implied by the ones we prove for this restricted
class. After the preliminary result in Lemma 3, we focus on proving the
hardness result for the problem of deciding whether $(G,[2])$ has an Eulerian
local tour, explaining the construction behind our reduction from NAE 3-,
whose correctness is proved in Theorem 3.1. We also argue that, if $G$ is a
cubic graph, then being Hamiltonian is a necessary but not sufficient
condition for $(G,[2])$ to admit an Eulerian local tour, arguing the need of
an _ad hoc_ reduction for our problem. As the reduction in Theorem 3.1 focuses
on solving Eulerian Local Tour for $\tau=2$, in Corollary 2 we extend this
result to each fixed $\tau$ and to trails, thus completing the proof of Item 2
of Theorem 1.1.

The following Lemma, whose proof is given in Appendix 0.B.2, will help our
proof.

######  Lemma 3.

Let $G$ be a graph. If $(G,[2])$ has an Eulerian local tour $T$, then $T$
restricted to timestamp $i$ must pass by all vertices of odd degree in $G$,
for each $i\in[2]$.

A simple consequence of the above lemma is that, as previously said, if $G$ is
cubic, then $G$ must Hamiltonian in order for $(G,[2])$ to have an Eulerian
local tour. Since deciding whether a cubic graph is Hamiltonian is
$\NP$-complete [12], this hints towards the $\NP$-completeness of the problem.
However, since the other way around is not necessarily true (see e.g. the
graph in Figure 2), we need an explicit reduction. Indeed, the construction in
Figure 2 shows us that we might need an arbitrarily large lifetime in order to
be able to visit all the edges of $(G,[\tau])$ even if $G$ is a 2-connected
outerplanar cubic graph (which is trivially hamiltonian).

Figure 2: Example of outerplanar graph $G$ such that $(G,[2])$ does not have
an Eulerian local tour.

In the following we explain the construction behind our reduction from NAE 3-.
Let $\phi$ be a CNF formula on variables $\\{x_{1},\cdots,x_{n}\\}$ and
clauses $\\{c_{1},\cdots,c_{m}\\}$. We start by presenting a meta-
construction, in the sense that part of the constructed graph will be
presented for now as black boxes and the actual construction is done later, as
depicted in Figure 3. The meta part concerns the clauses; so for now, denote
by $C_{i}$ the black box related to clause $c_{i}$. Without going into
details, $C_{i}$ will contain exactly one entry vertex for each of its
literal. So, given a literal $\ell$ contained in $c_{i}$, denote by
$I_{i}(\ell)$ the entry vertex for $\ell$ in $C_{i}$. All defined three
vertices are distinct.

$C_{i_{1}}$$C_{i_{2}}$$\cdots$$C_{i_{p}}$$C_{j_{1}}$$C_{j_{2}}$$\cdots$$C_{j_{q}}$$I_{i}$$O_{i}$
Figure 3: Edge gadget with clause black boxes.

Now, for each variable $x_{i}$, let $c_{i_{1}},\cdots,c_{i_{p}}$ be the
clauses containing $x_{i}$ positively and $c_{j_{1}},\cdots,c_{j_{q}}$
containing $x_{i}$ negatively. Add two new vertices, $I_{i}$ and $O_{i}$
(these will be the entry and exit vertices for the variable gadget), and add
the following edges (these compose the paths shown in Figure 3):

$\begin{array}[]{ll}E_{i}=&\\{I_{i}I_{i_{1}}(x_{i}),I_{i}I_{j_{1}}(\overline{x}_{i}),I_{i_{p}}(x_{i})O_{i},I_{j_{q}}(\overline{x}_{i})O_{i}\\}\\\
&\cup\\{I_{i_{h}}(x_{i})I_{i_{h+1}}(x_{i})\mid h\in[p-1]\\}\\\
&\cup\\{I_{j_{h}}(\overline{x}_{i})I_{j_{h+1}}(\overline{x}_{i})\mid
h\in[q-1]\\}\end{array}$

The paths will function as a switch, telling us whether the variable is true
or false within the considered snapshot; we then denote by $P_{i}$ the set of
edges in the path $(I_{i},I_{i_{1}}(x_{i}),\cdots,I_{i_{p}}(x_{i}),O_{i})$,
and by $\overline{P}_{i}$ the set of edges in the path
$(I_{i},I_{j_{1}}(\overline{x}_{i}),\cdots,I_{j_{q}}(\overline{x}_{i}),O_{i})$.
Now, to link the variable gadgets and to construct the clause gadgets, we will
need a gadget that will function as an edge that must appear in the trail
performed in $G_{1}$ and the one performed in $G_{2}$. For this, we use Lemma
3 applied to the gadget in Figure 4(a); when adding such a gadget between a
pair $u,v$, we simply say that we are adding the _forced edge $uv$_.

Now, to link the variable gadgets, we add three new vertices $s_{1},s_{2},t$
and the following forced edges.

$E^{\prime}=\\{s_{i}t\mid
i\in[2]\\}\cup\\{tI_{1},O_{n}t\\}\cup\\{O_{i}I_{i+1}\mid i\in[n]\\}.$

The new vertices simply help us assume where the trail starts and finishes.
Now, let $T$ be an Eulerian local tour of $(G,[2])$ and denote by $T_{i}$ the
trail in $G$ defined by $T$ restricted to $G_{i}$, for $i\in[2]$. It is fairly
easy to see (and we will prove it shortly) that if we can ensure that $T_{1}$
uses $P_{i}$ if and only if $T_{2}$ uses $\overline{P}_{i}$, then we can prove
equivalence with NAE 3-. In other words, the clause gadget must be so that,
for every clause $c_{j}$ containing $x_{i}$ (or equivalently
$\overline{x}_{i}$), we get that either both edges incident to $I_{j}(x_{i})$
in $P_{i}$ (or equivalently $I_{j}(\overline{x}_{i})$ in $\overline{P}_{i}$)
are used, or none of them is used. Such a gadget is presented in Figure 4(b),
where the red edges are forced.

$u$$v$ (a) Gadget related to a forced edge $uv$.

$c^{i}_{2}$$a^{i}_{1}$$a^{i}_{2}$$b^{i}_{1}$$b^{i}_{2}$$c^{i}_{1}$$I_{i}(\ell_{1})$$I_{i}(\ell_{2})$$I_{i}(\ell_{3})$
(b) Gadget related to clause $c_{i}$. Red edges represent forced edges.

Figure 4: Gadgets for the reduction in Theorem 3.1.

######  Theorem 3.1.

Let $G$ be a graph, with degree is at most 4. Then Eulerian Local Tour is
$\NP$-complete on $(G,[2])$.

###### Proof.

Let $\phi$ and $G$ be as previously stated. First, consider a truth NAE
assignment $f$ to $\phi$. We construct $T_{1},T_{2}\subseteq E(G)$ and prove
that they form an Eulerian local tour of $G$. Start by putting $P_{i}$ in
$T_{1}$ and $\overline{P}_{i}$ in $T_{2}$ if $x_{i}$ is true, and the other
way around if $x_{i}$ is false. From now on, whenever we add a forced edge to
$T_{1}$ and $T_{2}$, we are actually adding the trails depicted in Figure 5.

$u$$v$ (a) Trail added to $T_{1}$ when forced edge $uv$ is added to $T_{1}$.

$u$$v$ (b) Trail added to $T_{2}$ when forced edge $uv$ is added to $T_{2}$.

Figure 5: Trails related to forced edges.

Now, add $E^{\prime}$ to both $T_{1}$ and $T_{2}$, and consider $c_{i}$ with
literals $\ell_{1},\ell_{2},\ell_{3}$. Suppose, without loss of generality,
that $\ell_{1}$ is true and $\ell_{2}$ is false. We then add to $T_{1}$ the
trail depicted in Figure 6(a), and to $T_{2}$ the one depicted in Figure 6(b).
Observe that all internal edges of $C_{i}$ are covered; also note that the
edges incident to $I_{i}(\ell_{3})$ that come from the variable must be either
in $T_{1}$ or in $T_{2}$ by construction. We know that the remaining edges are
also covered by $T_{1}\cup T_{2}$ by construction. Finally, notice that both
$T_{1}$ and $T_{2}$ touch all odd-degree vertices in a way that every vertex
(including the even-degree ones) has even degree in $T_{1}$ and in $T_{2}$,
except $s_{1},s_{2}$ which have degree exactly 1. Also note that they form a
connected graph; indeed they are formed by the cycle passing through the
variable gadgets and $t$, together with some pending trails passing by the
clause gadgets. Therefore, we can find an $s_{1},s_{2}$-trail passing by all
edges of $T_{1}$, and an $s_{2},s_{1}$-trail passing by all edges of $T_{2}$,
thus getting our Eulerian local tour.

$c^{i}_{2}$$a^{i}_{1}$$a^{i}_{2}$$b^{i}_{1}$$b^{i}_{2}$$c^{i}_{1}$$I_{i}(\ell_{1})$$I_{i}(\ell_{2})$$I_{i}(\ell_{3})$
(a) Added to $T_{1}$ when $\ell_{1}$ is true.

$c^{i}_{2}$$a^{i}_{1}$$a^{i}_{2}$$b^{i}_{1}$$b^{i}_{2}$$c^{i}_{1}$$I_{i}(\ell_{1})$$I_{i}(\ell_{2})$$I_{i}(\ell_{3})$
(b) Added to $T_{2}$ when $\ell_{2}$ is false.

Figure 6: Trails in $C_{i}$ related to a given NAE assignment.

Now, let $T$ be an Eulerian local tour of $(G,[2])$, and for each $i\in[2]$,
denote by $T_{i}$ the trail in $G$ defined by $T$ restricted to $G_{i}$. First
observe that Lemma 3 indeed ensures that $T_{1}$ and $T_{2}$ restricted to the
gadget related to a forced edge $uv$ must be exactly as the trails depicted in
Figure 5; so in what follows we treat them exactly like edges that must appear
in $T_{1}$ and $T_{2}$. Since there are 2 vertices of degree 1, namely $s_{1}$
and $s_{2}$, by Lemma 3 we can suppose that $T_{1}$ starts in $s_{1}$ and
finishes in $s_{2}$, while $T_{2}$ starts in $s_{2}$ and finishes in $s_{1}$.
Therefore, each of $T_{1}$ and $T_{2}$ contains a tour in
$G-\\{s_{1},s_{2}\\}$, and hence:

  * (I)

Every vertex $u\in V(G)\setminus\\{s_{1},s_{2}\\}$ has even degree in both
$T_{1}$ and $T_{2}$.

Now, we prove that for each $x_{i}$, if $T_{1}$ intersects $P_{i}$, then
$T_{1}$ does not intersect $\overline{P}_{i}$, the same holding for $T_{2}$.
For this, consider $c_{i}$ with literals $\ell_{1}$, $\ell_{2}$, $\ell_{3}$,
and for each $j\in[3]$, denote by $e^{1}_{i}(\ell_{j}),e^{2}_{i}(\ell_{j})$
the edges incident to $I_{i}(\ell_{j})$ not contained in $C_{i}$. We first
prove that, for each $j\in[3]$ and $k\in[2]$:

  * (II)

Edge $e^{1}_{i}(\ell_{j})$ is used in $T_{k}$ if and only if edge
$e^{2}_{i}(\ell_{j})$ is used in $T_{k}$.

Without loss of generality, assume $j=1$ and let $k\in[2]$; suppose by
contradiction that $e^{1}_{i}(\ell_{1})$ is in $T_{k}$, while
$e^{2}_{i}(\ell_{1})$ is not in $T_{k}$. By (I) we get that either
$a_{1}^{i}I_{i}(\ell_{1})$ or $a_{2}^{i}I_{i}(\ell_{1})$ is in $T_{k}$, say
$a_{1}^{i}I_{i}(\ell_{1})$. But then we get that $a_{1}^{i}a_{2}^{i}$ and
$a^{i}_{2}I_{i}(\ell_{1})$ are not in $T_{k}$, as otherwise either $a_{1}^{i}$
or $I_{i}(\ell_{1})$ would have odd degree in $T_{k}$. This is a contradiction
since we then get $a^{i}_{2}$ with degree 1 in $T_{k}$. The same argument can
be analogously applied when $a_{2}^{i}I_{i}(\ell_{1})$ is in $T_{k}$ or when
$j\in\\{2,3\\}$.

Consider now a variable $x_{i}$, and let $c_{i_{1}},\cdots,c_{i_{p}}$ be the
clauses containing $x_{i}$ positively and $c_{j_{1}},\cdots,c_{j_{q}}$
containing $x_{i}$ negatively. Because $I_{i}$ has degree 3 in $G$ and
$E^{\prime}\setminus\\{s_{i}t\mid i\in[2]\\}$ is contained in $T_{1}$, we get
that exactly one between $I_{1}I_{i_{1}}(x_{i})$ and
$I_{1}I_{j_{1}}(\overline{x}_{i})$ is contained in $T_{1}$. From (II), we then
get that either $P_{i}$ is contained in $T_{1}$ or $\overline{P}_{i}$ is
contained in $T_{1}$. Observe that this implies that $P_{i}$ is contained in
$T_{k}$, while $\overline{P}_{i}$ is contained in $T_{3-k}$, for some
$k\in[2]$. We then set $x_{i}$ to be true if and only if $T_{1}$ contains
$P_{i}$. Because the edges in
$E_{i}=\\{e^{1}_{i}(\ell_{j}),e^{2}_{i}(\ell_{j})\mid j\in[3]\\}$ separate
$C_{i}$ from the rest of the graph and by (I), we get that both $T_{1}$ and
$T_{2}$ must intersect $E_{i}$. Finally by (II) we get that this assignment is
a NAE truth assignment for $\phi$. ∎

Observe that if we add two new vertices of degree one adjacent to vertex $t$,
then we get a reduction to the problem of deciding whether the edges of $G$
can be covered by two trails, proving Corollary 1.

######  Corollary 2.

Eulerian Local Tour and Eulerian Local Trail are $\NP$-complete on temporal
graphs with lifetime $\tau$ for every fixed $\tau\geq 2$. This also holds on
dynamic-based graphs.

###### Proof.

We first make a reduction from Eulerian Local Tour on $(G,[2])$ to Eulerian
Local Trail on $(G^{\prime},[\tau])$. Given $(G,[2])$, let $G^{\prime}$ be
obtained from $G$ by adding a star on $\tau+1$ vertices and identifying one of
its leaves with a vertex $s\in V(G)$. We argue that $(G,[2])$ has an Eulerian
local tour starting and finishing in $s$ if and only if $(G^{\prime},[\tau])$
has an Eulerian local trail. The lemma follows because we can then obtain a
Turing reduction by building a distinct instance for each $s\in V(G)$. Denote
the vertices of the initial star by $u,v_{1},\cdots,v_{\tau+1}$, where $u$ is
the central vertex, and $v_{2}$ is the vertex where $G$ is pending. Let $T$ be
an Eulerian local tour of $(G,[2])$ starting and finishing in $s$, and
$T_{1},T_{2}$ be the trails in $G$ defined by $T$. Build an Eulerian local
trail of $(G^{\prime},[\tau])$ by visiting $v_{1}u$, $v_{2}u$ and $T_{1}$ in
$G^{\prime}_{1}$, then performing $T_{2}$ and visiting $v_{2}u$ and $uv_{3}$
in $G^{\prime}_{2}$, and finish to visit the remaining edges of the star in
the obvious way.

Now, let $T$ be an Eulerian local trail of $(G^{\prime},[\tau])$, and denote
by $T_{i}$ the trail in $G^{\prime}$ defined by $T$ restricted to
$G^{\prime}_{i}$, for each $i\in[\tau]$. Observe that because we have $\tau+1$
cut edges, we get that each $T_{i}$ contains at most 2 of them, and in case it
contains exactly 2, say $v_{1}u,uv_{2}$, then $T_{i+1}$ either does not
contain any cut edge, or must intersect $T_{i}$ in $v_{1}u,uv_{2}$. This means
that the best we can do in order to finish by time $\tau$ is to visit exactly
two of them in the first snapshot, and exactly one more in each of the
subsequent snapshots. We can therefore suppose, without loss of generality
that $T_{i}$ contains $v_{i}u,v_{i+1}u$ for each $i\in[\tau]$. Note that this
implies that every edge of $(G,[2])$ must be visited in $T_{1}$ and $T_{2}$,
with $T_{1}$ starting in $v_{2}$ and $T_{2}$ finishing in $v_{2}$, as we
wanted to prove.

Finally, note that $(G^{\prime},[\tau])$ constructed above has an Eulerian
local trail if and only if $(G^{\prime},\tau+1)$ has an Eulerian local tour.
This completes our proof. ∎

##  4 Eulerian Tours and Trails

We finally focus on Eulerian Trail and Eulerian Tour, proving that in the
general case they are both -complete, hence, proving Item 3 in Theorem 1.1. To
this aim, we make an adaptation of the construction in Theorem 3.1. Observe
that here the base graph needs to be Eulerian as otherwise the answer to
Eulerian Trail is trivially NO. This also implies that the problem restricted
to dynamic-based graphs is trivial: if the base graph is Eulerian, then the
answer to Eulerian Tour is YES; otherwise, then the answer is NO. The trick
now is to take advantage of the function $\lambda$ in order to enforce the
edges.

######  Theorem 4.1.

Eulerian Tour and Eulerian Trail are $\NP$-complete, even on temporal graphs
with fixed lifetime $\tau\geq 2$.

###### Proof.

We first prove the case $\tau=2$. For this, we simply replace the gadget to
enforce an edge $uv$ in the construction of Section 3 by two paths of length
2, $P^{1}_{uv}$ and $P^{2}_{uv}$, where the edges in $P^{i}_{uv}$ are active
only in snapshot $G_{i}$, for each $i\in[2]$. Because the arguments used in
Section 3 depended only on the fact of $uv$ be indeed an enforced edge, we can
apply the same arguments here. The only difference is that the trails in
$G_{1}$ and $G_{2}$ now cannot intersect, which indeed is the case since the
intersection between $T_{1}$ and $T_{2}$ in Section 3 is exactly the set of
forced edges, and since here each appearance of a forced edge $uv$ is actually
related either to $P^{1}_{uv}$ or to $P^{2}_{uv}$.

Now, in order to prove the $\NP$-completeness for higher values of $\tau$, we
can simply add new vertices $v_{3},\cdots,v_{\tau}$ and edges
$\\{s_{1}v_{3}\\}\cup\\{v_{i}v_{i+1}\mid i\in\\{3,\cdots,\tau\\}\\}$. This
gives us that Eulerian Trail is $\NP$-complete on $(G,\lambda)$ with lifetime
$\tau$ for every fixed $\tau\geq 2$. And if we want a closed trail, it
suffices to identify $v_{\tau}$ with $s_{1}$, if $\tau\geq 4$, and if
$\tau=3$, we add a new vertices $v_{4}$ and edges $v_{3}v_{4},v_{4}s_{1}$
active in snapshot $G_{3}$. This concludes our proof. ∎

## References

  * [1] S Arumugam, I Hamid, and VM Abraham.  Decomposition of graphs into paths and cycles.  Journal of Discrete Mathematics, 2013, 2013.
  * [2] Pierre Borgnat, Eric Fleury, Jean-Loup Guillaume, Clémence Magnien, Céline Robardet, and Antoine Scherrer.  Evolving networks.  In Mining Massive Data Sets for Security, pages 198–203, 2007.
  * [3] Arnaud Casteigts, Paola Flocchini, Walter Quattrociocchi, and Nicola Santoro.  Time-varying graphs and dynamic networks.  International Journal of Parallel, Emergent and Distributed Systems, 27(5):387–408, 2012.
  * [4] Arnaud Casteigts, Anne-Sophie Himmel, Hendrik Molter, and Philipp Zschoche.  Finding temporal paths under waiting time constraints.  In Yixin Cao, Siu-Wing Cheng, and Minming Li, editors, 31st International Symposium on Algorithms and Computation, ISAAC 2020, December 14-18, 2020, Hong Kong, China (Virtual Conference), volume 181 of LIPIcs, pages 30:1–30:18. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2020.
  * [5] Merve Kayacı Çodur and Mustafa Yılmaz.  A time-dependent hierarchical chinese postman problem.  Central European Journal of Operations Research, 28(1):337–366, 2020.
  * [6] Marek Cygan, Dániel Marx, Marcin Pilipczuk, Michał Pilipczuk, and Ildikó Schlotter.  Parameterized complexity of eulerian deletion problems.  Algorithmica, 68(1):41–61, 2014.
  * [7] Elias Dahlhaus, David S. Johnson, Christos H. Papadimitriou, Paul D. Seymour, and Mihalis Yannakakis.  The complexity of multiterminal cuts.  SIAM Journal on Computing, 23(4):864–894, 1994.
  * [8] Thomas Erlebach, Michael Hoffmann, and Frank Kammer.  On temporal graph exploration.  In Magnús M. Halldórsson, Kazuo Iwama, Naoki Kobayashi, and Bettina Speckmann, editors, Automata, Languages, and Programming - 42nd International Colloquium, ICALP 2015, Kyoto, Japan, July 6-10, 2015, Proceedings, Part I, volume 9134 of Lecture Notes in Computer Science, pages 444–455. Springer, 2015.
  * [9] Thomas Erlebach and Jakob T Spooner.  Faster exploration of degree-bounded temporal graphs.  In 43rd International Symposium on Mathematical Foundations of Computer Science (MFCS 2018). Schloss Dagstuhl-Leibniz-Zentrum fuer Informatik, 2018.
  * [10] Thomas Erlebach and Jakob T. Spooner.  Non-strict temporal exploration.  In Andrea Werneck Richa and Christian Scheideler, editors, Structural Information and Communication Complexity - 27th International Colloquium, SIROCCO 2020, Paderborn, Germany, June 29 - July 1, 2020, Proceedings, volume 12156 of Lecture Notes in Computer Science, pages 129–145. Springer, 2020.
  * [11] Fedor V Fomin and Petr A Golovach.  Long circuits and large euler subgraphs.  SIAM Journal on Discrete Mathematics, 28(2):878–892, 2014.
  * [12] Michael R Garey, David S. Johnson, and R Endre Tarjan.  The planar hamiltonian circuit problem is np-complete.  SIAM Journal on Computing, 5(4):704–714, 1976.
  * [13] MEI Gu GLJAN.  Graphic programming using odd or even points, chinese j.  Mufh, 1:273–277, 1962.
  * [14] Renzo Gómez and Yoshiko Wakabayashi.  Covering a graph with nontrivial vertex-disjoint paths: existence and optimization.  In International Workshop on Graph-Theoretic Concepts in Computer Science, pages 228–238. Springer, 2018.
  * [15] David Kempe, Jon Kleinberg, and Amit Kumar.  Connectivity and inference problems for temporal networks.  In STOC ’00: Proceedings of the thirty-second annual ACM symposium on Theory of computing, 2000.
  * [16] Matthieu Latapy, Tiphaine Viard, and Clémence Magnien.  Stream graphs and link streams for the modeling of interactions over time.  Social Network Analysis and Mining, 8:61, 2018.
  * [17] Paul Manuel.  Revisiting path-type covering and partitioning problems.  arXiv preprint arXiv:1807.10613, 2018.
  * [18] Othon Michail.  An introduction to temporal graphs: An algorithmic perspective.  Internet Mathematics, 12(4):239–280, 2016.
  * [19] Othon Michail and Paul G Spirakis.  Traveling salesman problems in temporal graphs.  Theoretical Computer Science, 634:1–23, 2016.
  * [20] James B Orlin.  Some problems on dynamic/periodic graphs.  In Progress in Combinatorial Optimization, pages 273–293. Elsevier, 1984.
  * [21] Nicola Santoro, Walter Quattrociocchi, Paola Flocchini, Arnaud Casteigts, and Frédéric Amblard.  Time-varying graphs and social network analysis: Temporal indicators and metrics.  CoRR, abs/1102.0629, 2011.
  * [22] Jinghao Sun, Guozhen Tan, and Honglei Qu.  Dynamic programming algorithm for the time dependent chinese postman problem.  Journal of Information and Computational Science, 8:833–841, 2011\.
  * [23] Hsiao-Fan Wang and Yu-Pin Wen.  Time-constrained chinese postman problems.  Computers & Mathematics with applications, 44(3-4):375–387, 2002\.
  * [24] Huanhuan Wu, James Cheng, Silu Huang, Yiping Ke, Yi Lu, and Yanyan Xu.  Path problems in temporal graphs.  Proceedings of the VLDB Endowment, 7(9):721–732, 2014.

##  Appendix 0.A A brief summary about Eulerian tours in dynamic graphs

A _dynamic graph_ is a pair $(G,T)$ where $G$ is a finite digraph and $T$ is a
function $T:E(G)\rightarrow{\mathbb{Z}}$, called _transit time function_. A
dynamic graph can also be seen as a special type of infinite digraph ${\cal
G}$, where $V({\cal G})=V(G)\times{\mathbb{Z}}$, and $(u,i)(v,j)\in
E(\cal{G})$ if and only if $uv\in E(G)$ and $T(uv)=j-i$. Observe that the
transit time of an arc can also be negative, and therefore there might exist
arcs going from a vertex $(u,i)$ to a vertex $(v,j)$ with $j<i$, which in the
temporal graph context would be considered as going back in time. An _Eulerian
trail_ in $(G,T)$ is a trail that passes through all the edges of ${\cal G}$.
More formally, it is a function $f:{\mathbb{Z}}\rightarrow
V(G)\times{\mathbb{Z}}$ such that $f(i)f(i+1)\in E({\cal G})$ for every
$i\in{\mathbb{Z}}$, and for every $(u,i)(v,j)\in E({\cal G})$, there exists a
unique $\ell$ such that $(u,i)(v,j)$ is equal to $f(\ell)f({\ell+1})$.

Recall that a digraph $G$ is Eulerian if and only: (i) $G$ has at most one
non-trivial component; and (ii) the indegree of $u$, denoted with $d^{-}(u)$,
is equal to its outdegree, denoted as $d^{+}(u)$, for every $u\in V(G)$.
Observe that these conditions are also trivially necessary for the infinite
case. Also, note that, given $u\in V(G)$, each in-arc $uv$ of $G$ incident to
$u$ gives rise to exactly one in-arc $(u,i)(v,i+T(u,v))$ incident to $(u,i)$,
for every $i\in{\mathbb{Z}}$. The same clearly holds for every out-arc.
Therefore, one can see that ${\cal G}$ satisfies (i) and (ii) if and only if
$G$ satisfies (i) and (ii). However, as proved in [20], these are not the only
necessary conditions. Nevertheless, a characterization is still possible, with
an additional, also easy to test, condition.

######  Theorem 0.A.1 ([20]).

Let $(G,T)$ be a dynamic graph. Then $(G,T)$ has an Eulerian trail if and only
the following conditions hold:

  1. 1.

$d^{-}(u)=d^{+}(u)$ for every $u\in V(G)$;

  2. 2.

$G$ is connected; and

  3. 3.

$\sum_{e\in E(G)}T(e)\in\\{-1,1\\}$

Even if the necessary part of the proof is more technical, the sufficiency
part is quite natural, because condition (3) tells us that, given an eulerian
tour $T=(v_{1},\cdots,v_{m},v_{1})$ of $G$ and fixing a time $i$, we can use
$T$ to traverse all the edges incident to $(v_{1},i)$ in a way that we arrive
in $(v_{1},i+1)$ (or $(v_{1},i-1)$ if the sum is $-1$) just in time to apply
the same process to $(v_{1},i+1)$. Because $(v_{1},i)$ is chosen arbitrarily,
we are ensured to visit all edges of ${\cal G}$.

##  Appendix 0.B Proofs

###  0.B.1 Proof of Lemma 1

###### Proof.

Let $G_{1},\cdots,G_{\tau}$ be the snapshots of $G$; note first that if
$E(G_{i})$ is empty, then this snapshot can be suppressed. Our problem reduces
to deciding whether there is a choice of connected components
$H_{1},\ldots,H_{\tau}$, one for each timestamp $i$, that together cover all
the edges of $G$ and is such that $H_{i}$ intersects $H_{i+1}$, for each
$i\in[\tau-1]$. As for each $i\in[\tau]$, there are at most $n$ nodes in the
intersections, there are at most $O(n^{\tau-1})$ choices. For each choice the
test can be done in $O(\tau(n+m))$, obtaining $O(\tau(n+m)n^{\tau-1})$ running
time, which is $O((n+m)n^{\tau-1})$. ∎

###  0.B.2 Proof of Lemma 3

###### Proof.

For each $i\in[2]$, denote by $T_{i}$ the trail in $G$ equal to $T$ restricted
to timestamp $i$, and suppose, by contradiction, that $u\in V(G)$ is a vertex
with odd degree not contained in $T_{1}$. Because $T$ is a temporal tour,
observe that $T_{1}$ is a trail in $G$ starting at some $s$ and finishing at
some $t$, and $T_{2}$ is a trail in $G$ starting at $t$ and finishing at $s$,
with possibly $s=t$. This means that the subgraph of $G$ formed by the edges
of $T_{2}$ is such that every $x\in V(G)\setminus\\{s,t\\}$ has even degree.
This is a contradiction because, since no edge incident to $u$ is visited in
$T_{1}$, we get that all the edges incident to $u$ must be visited in $T_{2}$,
i.e., $u$ would have odd degree in $T_{2}$. The same argument holds in case
$u$ is not in $T_{2}$, and the lemma follows. ∎


Interpreting artificial neural networks to detect genome-wide association
signals for complex traits

Burak Yelmen1,*, Maris Alver1, Estonian Biobank Research Team1, Flora Jay2,
Lili Milani1

1 Estonian Genome Centre, Institute of Genomics, University of Tartu, Tartu,
Estonia

2 CNRS, INRIA, LISN, Paris-Saclay University, Orsay, France

<EMAIL_ADDRESS>
## Abstract

Investigating the genetic architecture of complex diseases is challenging due
to the highly polygenic and interactive landscape of genetic and environmental
factors. Although genome-wide association studies (GWAS) have identified
thousands of variants for multiple complex phenotypes, conventional
statistical approaches can be limited by simplified assumptions such as
linearity and lack of epistasis models. In this work, we trained artificial
neural networks for predicting complex traits using both simulated and real
genotype/phenotype datasets. We extracted feature importance scores via
different post hoc interpretability methods to identify potentially associated
loci (PAL) for the target phenotype. Simulations we performed with various
parameters demonstrated that associated loci can be detected with good
precision using strict selection criteria, but downstream analyses are
required for fine-mapping the exact variants due to linkage disequilibrium,
similarly to conventional GWAS. By applying our approach to the schizophrenia
cohort in the Estonian Biobank, we were able to detect multiple PAL related to
this highly polygenic and heritable disorder. We also performed enrichment
analyses with PAL in genic regions, which predominantly identified terms
associated with brain morphology. With further improvements in model
optimization and confidence measures, artificial neural networks can enhance
the identification of genomic loci associated with complex diseases, providing
a more comprehensive approach for GWAS and serving as initial screening tools
for subsequent functional studies.

#### Keywords:

Deep learning, interpretability, genome-wide association studies, complex
diseases

## Introduction

Understanding the genetic component of complex diseases is challenging due to
high polygenicity and confounding factors defining these traits. Genome-wide
association studies (GWAS) have been crucial for the discovery of thousands of
disease-related genomic loci in the past two decades [1], and they continue to
be relevant thanks to ever-expanding biobanks with deep phenotype/genotype
datasets, allowing the detection of the most minuscule but statistically
significant signals. Despite the immense success and broad applications over
all types of traits, conventional GWAS methodology still relies on linear
models (i.e., logistic and linear regression) for detecting such signals. More
recently, linear mixed models have been utilised to incorporate random and
fixed effects, especially to account for relatedness and population structure
without data pruning [2, 3]. Although such improvements seem to improve GWAS
outcomes generally, they do not diverge from the fundamental assumptions of
linear and additive genetic models for complex traits. A natural alternative
for the nonlinear modelling of such traits would be machine learning, and
especially deep learning, since sufficiently large neural networks can
approximate any bounded continuous function as defined in the universal
approximation theorem [4]. Hence, not surprisingly, various models have been
developed for the prediction of disease status from genetic data following the
success of deep learning applications in many domains [5, 6]. Despite the
state-of-the-art performance of these models in prediction, a main drawback of
utilising such complex architectures is the inherent difficulty of
interpretability due to the high number of parameters at various abstraction
levels. In this context, the development and integration of interpretability
methodology are not only important for understanding model biases and
improving architectures but also for detecting novel features (e.g., genomic
loci, haplotypes, or variants) associated with the target phenotype. Various
approaches have been proposed to tackle interpretability for the potential
discovery of causal genomic loci, ranging from intrinsically interpretable
architectures based on domain knowledge, such as functional gene annotations
[7, 8, 9], to post hoc analysis of trained models [10, 11, 12, 13]. Despite
these important advancements, multiple challenges remain for the widespread
use of interpretable deep learning in genomics. For instance, many of the
proposed methods rely on either specific architectures, specific
interpretability methods, or a priori biological knowledge. In this work, (i)
we propose a general framework for obtaining potentially associated loci (PAL)
from neural network models trained for phenotype prediction, (ii) compare the
performance of different post hoc and model-agnostic interpretability methods
and logistic regression using simulations and (iii) apply our approach to the
Estonian Biobank (EstBB) [14] schizophrenia (SCZ) cohort demonstrating its
utility for GWAS and reporting potentially novel loci not reported in
literature.

## Materials and Methods

### Data

We assembled the SCZ cohort from the EstBB dataset so that all patients (ICD
codes F20-F29) would have at least one antipsychotic prescription and have
their first diagnosis after 2006 with the diagnosis age between 14 and 60. For
SCZ cases, we first assembled a set of individuals without any mental,
behavioural and neurodevelopmental disorders (no ICD F* diagnosis) and without
any antipsychotic prescription. From this filtered set, controls were selected
to match the age, sex and BMI distribution of the case cohort with four
controls per case using MatchIt [15]. We used exact matching (method =
“exact”) for sex and nearest neighbour matching (method = “nearest”) for BMI
and age, resulting in 1814 cases (642 male, 1172 female) and 7325 controls
(2586 male, 4739 female). In addition, the whole dataset was pruned to
eliminate relatedness based on identity by descent (pi_hat $<$ 0.2). We
curated the genotype data of the EstBB SCZ case/control cohort derived from
the GSA array (Illumina’s GSA-MD-24v1, GSA-MD-24v2, ESTchip1_GSA-MD-24v2,
ESTchip2_GSA-MD-24v3, 2022-09-14 snapshot) with a final set of 290,522 SNPs
over 22 autosomal chromosomes after removing SNPs unique to at least one
version of the array, indels and filtering for bi-allelic regions.

### Simulations

We obtained the simulated phenotypes from the same combined SCZ case/control
genotype dataset. Instead of simulating genotypes, we simulated phenotypes to
keep genome distribution as realistic as possible. More in detail, let $X$
represent the genotype matrix where $X_{ij}$ is the genotype at position $j$
for individual $i$ with possible values -1 (AA), 0 (AB), 1 (BB); $\beta$ the
vector of genetic weights sampled from standard normal distribution $N$(0,1);
$\epsilon_{i}$ the noise term sampled from $N$(0, $k\sigma^{2}$) where $k$ is
a scaling factor and $\sigma^{2}$ is the variance of the total genetic
effects. Then the continuous phenotype $Y$ for individual $i$ is defined as:

$Y_{i}=\sum_{j\in D}f^{D}_{j}(X_{ij})+\sum_{j\in
R}f^{R}_{j}(X_{ij})+\sum_{(j,k)\in I}f_{jk}(X_{ij},X_{ik})+\sum_{(j,k,l)\in
I}f_{jkl}(X_{ij},X_{ik},X_{il})+\epsilon_{i}$

where the initial four terms define the genetic effects based on dominant,
recessive, two-way (between 2 SNPs) and three-way (between 3 SNPs)
interactions; $D$, $R$ and $I$ are disjoint sets of genomic positions for
dominant, recessive and interaction effects. For dominant and recessive effect
positions, contribution to the phenotype is defined by functions $f_{j}^{D}$
and $f_{j}^{R}$:

$f^{D}_{j}(X_{ij})=\begin{cases}\beta_{j}\cdot X_{ij}&\text{if
}X_{ij}\neq\text{-1}\\\ 0&\text{otherwise}\end{cases}$

$f^{R}_{j}(X_{ij})=\begin{cases}\beta_{j}\cdot X_{ij}&\text{if
}X_{ij}=\text{1}\\\ 0&\text{otherwise}\end{cases}$

Two-way and three-way interaction effects are defined by:

$f_{jk}(X_{ij},X_{ik})=\beta_{jk}\cdot X_{ij}\cdot X_{ik}$
$f_{jkl}(X_{ij},X_{ik},X_{il})=\beta_{jkl}\cdot X_{ij}\cdot X_{ik}\cdot
X_{il}$

Given a defined threshold $\tau$, the continuous phenotype is then binarized
by:

$B_{i}=\begin{cases}1&\text{if }Y_{i}>\tau\\\ 0&\text{otherwise}\end{cases}$

In practice, we chose $\tau$ to have approximately the same number of cases
($\sim$1814) as the real SCZ cohort. We simulated a total of 6 scenarios with
different noise scaling factors ($k$ = 1, $k$ = 2, $k$ = 3) and two different
causal SNP sizes (n = 100 and n = 1000) where all four genetic effect types
(dominant, recessive, two-way and three-way interactive) had positions with
ratios 5:5:8:2, respectively (i.e., for n = 100: 40, 10, 25, 25 and for n =
1000: 400, 100, 250, 250). Genotypes were represented by -1, 0, 1 instead of
0, 1, 2 to provide equal importance to either allele in a given position
(considering the neural network weights are initialised with Kaiming uniform
distribution, including negative and positive values uniformly). Nevertheless,
our preliminary trials with 0, 1, 2 representation provided similar results in
terms of model training.

### Neural network model

Models trained for binary trait classification were all feed forward fully
connected neural networks with an input layer, two hidden layers and an output
layer. The specific architecture consists of an initial dropout layer to
randomly mask a proportion of features (p = 0.99) at each batch training, a
linear layer with input size = 290,522 and output size = 290, a dropout layer
(p = 0.6) followed by ReLU activation, a linear layer with input size = 290
and output size = 29, a dropout layer (p = 0.6) followed by GELU activation, a
linear layer with input size = 29 and output size = 1 followed by sigmoid
activation function. Training was performed using cross-entropy loss function
and Adam optimizer with learning rate = 1e-5 and weight decay = 1e-3 for
regularisation. Initially, we tried a weighted loss function to compensate for
the class imbalance in the real dataset but this did not work (i.e., no
training run where the validation loss decreases), so we instead duplicated
the case samples and reduced the sample size of the controls by half with
random subsampling to match the case/control numbers [16]. We additionally
added random noise $N$(0,0.1) both to class labels (case label = 1, control
label = 0) and input genotypes to increase generalizability. All models were
trained for 1000 epochs with batch size = 256. All models were coded with
python-3.9 and pythorch-1.11 [17].

### Interpretability methods

There is a plethora of different interpretability methods for deep learning
models, but we focused on three post hoc and model agnostic approaches. Two of
them, gradients w.r.t input (saliency maps, SM) [18] and integrated gradients
(IG) [19], are based on gradients and applicable to any differentiable model.
SM creation is a commonly used and straightforward approach where gradients
w.r.t input are obtained to assess which features in a given input are
important for the output probability. For a binary classification neural
network model $f$ and input vector $x$, the saliency (with absolute values)
for the $j^{th}$ feature is defined as:

$SM_{j}(x)=\left|\frac{\partial f(x)}{\partial x_{j}}\right|$

IG is an improvement over the simpler gradient approach where gradients of the
model’s output with respect to the input are integrated along a path from a
baseline input to the actual input. For a neural network model $f$, input
vector $x$, baseline $x’$ and scalar $\alpha\in[0,1]$, integrated gradient for
the $j^{th}$ feature is defined as:

$IG_{j}(x)=(x_{j}-x^{\prime}_{j})\int_{0}^{1}\frac{\partial
f(x^{\prime}+\alpha(x-x^{\prime}))}{\partial x_{j}}\,d\alpha$

As the authors of the method emphasise, the choice of baseline vector is
important so that a meaningful interpolation can be achieved in the given
domain. In our application, we tested various baselines and decided on using a
vector of zeros (where the input is a case genotype) which in preliminary
analyses, provided better outcomes compared to (i) [input = case and control
genotypes] - [baseline = vector of average allele dosages for the whole
dataset] and (ii) [input = case genotypes] - [baseline = control genotypes]
settings. In this regard, our choice of genotype representations with -1, 0, 1
might also have been helpful since a vector of zeros would correspond to a
completely heterozygous genome, without preference over any allele in a given
position. We used Captum-0.6.0 python library [20] for calculating SM and IG.

Permutation-based feature importance (PM) is another model-agnostic technique
(in this case, the model does not even need to be differentiable) that
measures the impact of feature perturbations on model output [21]. For a
neural network model $f$, input vector $x$ and perturbed input vector
$x^{(j)}$ (where the $j^{th}$ feature is replaced by a random feature), the
permutation feature importance for the $j^{th}$ feature is defined as:

$PM_{j}(x)=\left|f(x)-f(x^{(j)})\right|$

Here, we opted to measure the rate of change in predictions for the feature
importance instead of a performance metric such as accuracy to increase
sensitivity.

Using these three methods, we obtained absolute feature importance scores for
a set of case samples (both in real and simulated scenarios) for all trained
models, averaged the scores over samples and applied min-max normalisation to
map the values to [0,1] range. We call this averaged value mean attribution
score (MAS). Moreover, the obtained case set consisted of genotypes with
neural network predicted probabilities above the 70th percentile. This
thresholding was performed mainly for computational ease but also for
utilising case samples with prediction values closer to 1. Different
percentile thresholds and using all true case genotypes provided similar
results in preliminary analyses.

### Feature importance to potentially associated loci

To account for the stochasticity of neural networks and reduce false positive
rate, we devised a streamlined method for obtaining potentially associated
loci (PAL) from the MAS computed from multiple models trained with different
seeds. More specifically, let $A$ be the MAS matrix sized $m\times L$ where
$m$ is the number of models (i.e., models trained with different seeds) and
$L$ is the length of the genotype. We obtain the average MAS vector $\mu$:

$\mu=\left(\mu_{1},\mu_{2},\ldots,\mu_{L}\right)$

where

$\mu_{j}=\frac{1}{m}\sum_{a=1}^{m}A_{aj}\quad\text{for }j=1,2,\ldots,L$

Then we set a significance threshold $\theta$ and obtain weights $w$ based on
the number of occurrences above the $\theta$ threshold among $m$ models
trained with different seeds:

$w_{j}=\frac{1}{m}\sum_{a=1}^{m}\mathbf{1}(A_{aj}>\theta),\quad\text{where
}\mathbf{1}\text{ is the indicator function}$

To obtain PAL, first we define the adjusted mean attribution score (AMAS) as:

$AMAS_{j}=w_{j}\cdot\mu_{j},\quad\text{for all }j\text{ such that
}\mu_{j}>\theta$

Then the PAL set contains all positions above the theta threshold based on
their AMAS:

$PAL_{AMAS}=\left\\{j\mid AMAS_{j}>\theta\right\\}$

We can also define PAL as positions above the theta threshold in all $m$
models:

$PAL_{Common}=\left\\{j\mid\forall a,A_{aj}>\theta\right\\}$

In practice, we chose $m=10$ for models trained with different seeds to be
able to perform a large number of trials with realistic computational
resources. Since there were multiple PAL in high LD detected by models trained
with different seeds, we also considered these positions ($r^{2}>0.5$) between
different models as common signals between models to increase detection power
both for $PAL_{Common}$ and $PAL_{AMAS}$. Furthermore, we set $\theta$ to
99.99 percentile (strict) and 99.95 percentile (relaxed) thresholds. Although
the choice of $\theta$ is essentially arbitrary, the strict threshold produces
30 PAL and the relaxed threshold 145 PAL, which are plausible numbers for
keeping the false positive rate low considering the expected high polygenicity
of schizophrenia [22]. An overview of our proposed approach is provided in Fig
1.

Figure 1: Overview of the approach for obtaining potentially associated loci
(PAL) from feature attribution scores obtained via SM, IG and PM methods.

### Logistic regression

To compare our method to conventional GWAS, we performed logistic regression
(LR) using all case (n = 1814) and control (n = 7325) samples (both in real
and simulated scenarios) using the top 3 principal components (PCs) from
principal component analysis (PCA) as covariates. We obtained the significance
threshold by dividing p-value = 0.05 by the number of tests (i.e., number of
SNPs: 290,522). We used python-3.9 and statsmodels 0.14.1 [23] for the code.
Furthermore, we performed an association analysis using REGENIE v2.2.4 [3] on
real data with the same case/control genotypes via the methodology described
by Pujol Gualdo et al. [24], except for the filtering steps to retain exactly
the same individuals and positions as LR.

Relevant code for the neural network model, interpretability methods and
logistic regression with mock simulated data can be found at
https://github.com/genodeco/Interpreting-ANNs-for-GWAS.

## Results

### Simulations

For 6 different phenotype simulations with different noise scaling factors
($k$ = 1, $k$ = 2, $k$ = 3) and different numbers of causal SNPs (100 and 1000
causal SNPs), we tested the performance of logistic regression (LR), saliency
map (SM), integrated gradient (IG) and permutation-based (PM) feature
importance approaches based on precision and ability to detect different types
of causal positions (interaction and dominant/recessive). An important note
here is that we defined true positive (TP) as a detected position (i.e., PAL)
which is in high correlation ($r^{2}>0.5$ in 100 SNP upstream and 100 SNP
downstream region, total $\sim$2Mb) with any true causal position (since high
LD, which is prevalent in proximal loci, makes it difficult to detect the
exact causal SNP). The results show that LR with Bonferroni correction seems
to produce the least number of false positives over all different simulations
but $PAL_{AMAS}$ and $PAL_{Common}$ obtained from the IG approach with strict
thresholding (99.99 percentile) was a close contender (S1 Fig). Interestingly
but perhaps expectedly, the LR approach predominantly detected the causal
positions with interactive genetic effects, whereas the highest percentage of
dominant/recessive positions were detected by gradient-based approaches (Table
1). We furthermore visualised the performance of IG on Manhattan plots
obtained from LR (S2 Fig). Even under high noise schemes, some high-effect
positions were detected both by LR and our method.

Table 1: Percentage of dominant/recessive positions (as opposed to interactive
positions) in correctly detected causal loci averaged over all 6 simulation
scenarios with varying parameters. Significance-based detection refers to
Bonferroni correction for LR (logistic regression) and strict threshold (99.99
percentile) for IG (integrated gradients), SM (saliency map) and PM
(permutation-based) approaches. Since significance thresholds are not directly
comparable between LR and other methods, percentages were also provided for
the top 5 and top 10 correctly detected causal loci. The expected frequency
for dominant/recessive positions based on the phenotype simulations is 0.5.
Values in parentheses demonstrate counts as
[dominant/recessive]/([dominant/recessive] + [interactive]).
Category | Top 5 | Top 10 | Significance-based
---|---|---|---
LR | 0.0000 (0/24) | 0.0213 (1/47) | 0.0902 (12/133)
IG | 0.1304 (3/23) | 0.2143 (9/42) | 0.1268 (9/71)
SM | 0.0000 (0/23) | 0.1500 (6/40) | 0.1549 (11/71)
PM | 0.0500 (1/20) | 0.1579 (6/38) | 0.1618 (11/68)

### Application to the EstBB SCZ cohort

After testing our approaches on simulations, we applied the IG method on the
EstBB SCZ cohort to detect PAL. Although it is difficult to compare the
methods reliably (see Discussion), we decided not to use SM and PM on real
data due to higher false positives under certain simulation scenarios (S1
Fig). We detected a single locus with the strict threshold and multiple loci
with the relaxed threshold in genomic regions with multiple protein-coding
genes (Fig 2). The highest confidence signal on chromosome 4 (PAL d) included
variants within TET2 and PPA2 genes and rs10010325 SNP, which has been
reported to have pleiotropic effects for attention deficit hyperactivity
disorder, autism spectrum disorder and intelligence [25]. Furthermore, in this
LD region resides rs7674220 (albeit not being present in our non-imputed array
framework), which is associated with SCZ and mapped to TET2 [26].
Interestingly, rs2647259 present in this region in our data seems to be
associated with PPA2 expression levels in brain tissues (based on variant-gene
association analysis, explained in the next paragraph). Other high-confidence
regions we detected in chromosomes 3 (PAL c) and 11 (PAL f) are gene-dense
regions with multiple low p-value SNPs in Trubetskoy et al. metadata (i.e.,
largest SCZ GWAS meta-analysis to date) [27]. The high signal LD region we
detect in chromosome 11 includes rs1892928 (not present in our non-imputed
array data), which was reported to be associated with SCZ and mapped to SYT12.
In this region, we also detected rs7122539 (present in our array framework),
which was shown to be associated with bipolar disorder and mapped to PC gene
[28, 29]. The region in chromosome 3 includes rs4955417 (associated with
depressive symptoms and neuroticism, mapped to IHO1 and C3orf62 [30, 31]) and
rs7617480 (associated with major depressive disorder and unipolar depression
[32, 33, 34]). Furthermore, QRICH1 gene reported in the Schizophrenia Exome
Sequencing Meta-Analysis (SCHEMA) with p-value$<$0.01 and reported to be
associated with developmental disorders and autism spectrum disorder resides
in this region [35, 36, 27, 37]. Other lower confidence but potentially
interesting regions include C1orf87 (PAL a) (with gene ontology annotation
calcium ion binding), and AOX1 (PAL b) (with gene ontology annotations
oxidoreductase and electron transfer activity) genes. SNPs in detected PAL are
presented in S1 Table, genes in detected PAL are presented in S2 Table, and
SNPs with low p-values ($<$0.001) from Trubetskoy et al. 2022 metadata in
detected high LD regions are presented in S3 Table.

Figure 2: PAL detected by integrated gradients (IG) approach. Red dashed lines
indicate significance thresholds (relaxed and strict) and blue markers
indicate PAL above threshold over all trained 10 models (i.e.,
$PAL_{Common}$). For all PAL (a-g), protein coding genes in those regions were
provided except for multi-gene regions ($>$2 genes) for which the most likely
mapped gene based on literature was provided inside the parentheses.

We further assessed whether variants in the PAL we detected were enriched in
genes with brain associations. For this purpose, we used the open access Adult
Genotype Tissue Expression (GTEx) Project [38] variant-gene association
dataset (GTEx_Analysis_v8_eQTL). Since we did not fine map the causal
variants, we assessed the significant expression levels (significance defined
as pval_nominal $<$ pval_nominal_threshold in significant variant-gene pairs
dataset from GTEx *.signif_variant_gene_pairs.txt.gz files) in all brain
tissues for all detected loci and close proximity loci (100 SNP upstream and
100 SNP downstream region of the PAL, total $\sim$2Mb) with high LD ($>$0.5).
We initially obtained the empirical cumulative distribution function (ECDF)
for significant gene-tissue (only brain tissues) expression change in all
290,522 positions. Mean value for all positions was 0.858 (i.e., on average, a
SNP had 0.858 gene-tissue combinations with significantly altered expression
in brain) with most positions (250,692 SNPs) having no variants significantly
affecting gene expression in brain tissues. Percentile ranks for all methods
detecting PAL was higher than 90, indicating potential enrichment in brain
expressed genes in these detected regions (S3 Fig). As an additional step, we
used the GENE2FUNC tool in FUMA [39] to assess functionality for the set of
genes present in the detected PAL. In terms of expression, there was
significant downregulation in brain tissues and in terms of GWAS enrichment,
lowest enrichment p-values were related to brain morphology and volume with
multiple involved genes (S4 Fig). It is important to note here that we
directly tested all the genes in PAL without any gene mapping attempt, with
most genes coming from gene-rich loci from chromosomes 3 and 11 (PAL c and f),
which might potentially be introducing bias to these results. Another
potential biasing caveat for both ECDF and functional analysis is that
detected PAL generally consist of multiple SNPs in high LD, and it is likely
that only a small portion of SNPs in these LD blocks indicate a true signal.

We furthermore compared these approaches with conventional GWAS using LR and
REGENIE association analysis on the real dataset. Although there was a high
correlation between MAS and -log(p) values (S5 Fig), there was no significant
signal detected using GWAS, and only a single PAL overlap was present in the
top 10 signals (S6 Fig).

## Discussion

In this work, we assessed post hoc neural network interpretability methods for
detecting disease-associated genomic loci and proposed a general framework for
using these methods within a genotype-phenotype context. Simulations
demonstrated that these approaches can provide low false positive rates with
strict thresholding even under relatively low genetic signal scenarios. An
important note here is that the noise factor of 3 in simulations approximately
corresponds to the SNP-based heritability estimate for SCZ, which is around
0.25 [27]. However, this estimate is based on linear additive models, whereas
our simulated traits were based on interactive and dominant/recessive effect
positions, possibly making detection of associated loci substantially more
difficult. This brings us to a fundamental challenge in method development for
causal detection, which is the lack of ground truth disease models. Without
knowing the exact disease models, it is not trivial to assess these methods
due to the likely discrepancy between simulated and real domains. In this
regard, a potential solution, knockoff variables, has been adopted for neural
networks in previous works [40, 13, 10]. Briefly, these are implanted “fake”
positions following the distribution of the real input features to assess the
false positive rates when detecting relevant features. However, sampling
knockoff variables is not straightforward because deviations from the
distribution of real input features would bias the false positive estimates.
Furthermore, implementing a sufficient number of knockoff variables would
potentially deteriorate the realness of data and significantly increase the
computational burden due to the fully connected architecture of our models.
Future parameter-light architectures (such as convolutions) could be more
suitable to explore this approach. Another challenge for method assessment
arises due to LD. Even through simulations with ground truth information,
evaluation metrics such as TP/FP counts can be substantially biased since
methods (including LR) tend to present highly correlated signals. With our
definition of TP (i.e., a detected position in high correlation with any true
causal position), a method might present higher TP (or FP) counts simply due
to multiple signals in a high LD block as opposed to another method with
signals consisting of isolated SNPs. Another LD-related challenge (but
specific to our application) is due to models trained with different seeds
detecting different signals in high LD. To increase our power, we counted
these high LD PAL as common signals between models, which boosted TP signals
in simulation analyses (S1 Fig). Although a substantial portion of this boost
was probably due to overlapping signals, the same TP/FP analysis where we only
considered exact detected SNPs as common between models demonstrated that our
initial approach (common signals based on LD) allows catching more TP signals
and FP signals do not change substantially between two approaches for the IG
method (S7 Fig). Considering all these issues, future research might benefit
from evaluation metrics better tailored for genomic data.

Another important subject related to neural networks is stochasticity. In our
approach, we employed both heavy dropout regularisation (including unusual
dropout masking before the input layer) and ensemble-like multiple seed
training to take advantage of stochasticity and limit false positives.
Interestingly, even with this heavy regularisation, we were able to train
models successfully (i.e., positive correlation between training and
validation loss), even though the predictive power of these models was
limited. Receiver operating characteristic (ROC) curves for phenotype
prediction in the real dataset were comparable between neural network and
logistic regression models, with logistic regression being marginally better
(S8 Fig). Initially, this might suggest that the real dataset primarily
exhibits linear effects. Yet, as a counter-evidence to the “mainly linear
effects” scheme, logistic regression could not detect any significant
positions when applied to the real cohort, and there was no correspondence
between the top detected positions of logistic regression and neural network
approaches (S6 Fig), which might suggest the existence of nonlinearities in
the real dataset. These aspects of effect types and how heavy regularisation
alters the capacity of neural networks for detecting nonlinearities are
interesting subjects for future research. Considering the relatively small
dataset for a case/cohort type analysis in this work, increasing the training
data and experimenting with more simulation scenarios (such as including more
nonlinearities, gene-environment interactions or covariate effects) could be
helpful to distinguish the differences between linear and nonlinear
approaches.

An intriguing finding was that all approaches mainly detected positions with
interactive effects while neural network methods detected more positions with
dominant/recessive effects compared to LR. This is potentially due to
dominant/recessive positions being more difficult to detect since the SNP
effect completely disappears for a proportion of genotypes (see Materials and
Methods). Another related point is that the outcome space of functions
defining the disease genomics is considerably smaller due to the ternary
nature of the genotype variables. Research on the topology of this space could
provide valuable insights on limitations of methods developed for in silico
detection of causal loci and provide better roadmaps during development.

Despite the state-of-the-art deep learning algorithms and the ever-increasing
amount of data in biobanks, there are still multiple challenges to the wide-
spread adoption of neural networks for GWAS, such as stochasticity stemming
from initial random states, lack of established uncertainty quantification and
lack of ground truth disease models. These challenges are also fundamentally
connected to very interesting questions, such as the contribution of nonlinear
effects (gene-gene or gene-environment interactions or other polynomial
disease models) or omnigenic vs polygenic models for complex disease
development [41]. In the near future, assumption-free interpretable neural
networks can become important tools to explore these questions and potentially
replace linear models as the conventional method of conducting GWAS and
obtaining interpretable polygenic risk scores.

## Acknowledgements

BY, MA, LM were supported by the European Union’s Horizon 2020 Research and
Innovation Programme under Grant agreements No 964874 Realment and 847776
CoMorMent. Thanks to the Agence Nationale de la Recherche
(ANR-20-CE45-0010-01) for funding. Thanks to the High Performance Computing
Center (HPC) of the University of Tartu for providing computational resources.
Thanks to the Estonian Biobank Research Team (Andres Metspalu, Tõnu Esko,
Reedik Mägi, Mari Nelis, Georgi Hudjashov) for the genotype and phenotype
datasets used in this work. Thanks to Liis Karo-Astover for administrative
support and proofreading. Data used for the gene-variant expression analyses
described in this manuscript were obtained from the Genotype-Tissue Expression
(GTEx) Portal on 04/04/2024.

## Contributions

Conceptualisation: BY. Data curation: BY, MA. Methodology: BY. Investigation:
BY, MA, FJ. Formal analysis: BY, MA. Visualisation: BY. Supervision: BY, LM.
Original draft: BY. Review & editing: BY, MA, FJ, LM. Funding acquisition: LM.

## Ethics statement

The activities of the EstBB are regulated by the Human Genes Research Act,
which was adopted in 2000 specifically for the operations of the EstBB.
Individual level data analysis was carried out under ethical approval
1.1-12/624 (issued 24.03.2020), 1.1-12/2618 (issued 04.08.2022) and
1.1-12/3454 (issued 20.10.2022) from the Estonian Committee on Bioethics and
Human Research (Estonian Ministry of Social Affairs).

## References

  * [1] A. Abdellaoui, L. Yengo, K. J. H. Verweij, and P. M. Visscher, “15 years of GWAS discovery: Realizing the promise,” The American Journal of Human Genetics, vol. 110, pp. 179–194, Feb. 2023.
  * [2] L. Jiang, Z. Zheng, H. Fang, and J. Yang, “A generalized linear mixed model association tool for biobank-scale data,” Nature Genetics, vol. 53, pp. 1616–1621, Nov. 2021.  Publisher: Nature Publishing Group.
  * [3] J. Mbatchou, L. Barnard, J. Backman, A. Marcketta, J. A. Kosmicki, A. Ziyatdinov, C. Benner, C. O’Dushlaine, M. Barber, B. Boutkov, L. Habegger, M. Ferreira, A. Baras, J. Reid, G. Abecasis, E. Maxwell, and J. Marchini, “Computationally efficient whole-genome regression for quantitative and binary traits,” Nature Genetics, vol. 53, pp. 1097–1103, July 2021.  Publisher: Nature Publishing Group.
  * [4] Y. Lu and J. Lu, “A universal approximation theorem of deep neural networks for expressing probability distributions,” in Proceedings of the 34th International Conference on Neural Information Processing Systems, NIPS ’20, (Red Hook, NY, USA), pp. 3094–3105, Curran Associates Inc., Dec. 2020.
  * [5] W. S. Alharbi and M. Rashid, “A review of deep learning applications in human genomics using next-generation sequencing data,” Human Genomics, vol. 16, p. 26, July 2022.
  * [6] X. Zhou, Y. Chen, F. C. F. Ip, Y. Jiang, H. Cao, G. Lv, H. Zhong, J. Chen, T. Ye, Y. Chen, Y. Zhang, S. Ma, R. M. N. Lo, E. P. S. Tong, V. C. T. Mok, T. C. Y. Kwok, Q. Guo, K. Y. Mok, M. Shoai, J. Hardy, L. Chen, A. K. Y. Fu, and N. Y. Ip, “Deep learning-based polygenic risk analysis for Alzheimer’s disease prediction,” Communications Medicine, vol. 3, pp. 1–20, Apr. 2023.  Publisher: Nature Publishing Group.
  * [7] A. van Hilten, S. A. Kushner, M. Kayser, M. A. Ikram, H. H. H. Adams, C. C. W. Klaver, W. J. Niessen, and G. V. Roshchupkin, “GenNet framework: interpretable deep learning for predicting phenotypes from genetic data,” Communications Biology, vol. 4, pp. 1–9, Sept. 2021.  Number: 1 Publisher: Nature Publishing Group.
  * [8] J. Xu, C. Mao, Y. Hou, Y. Luo, J. L. Binder, Y. Zhou, L. M. Bekris, J. Shin, M. Hu, F. Wang, C. Eng, T. I. Oprea, M. E. Flanagan, A. A. Pieper, J. Cummings, J. B. Leverenz, and F. Cheng, “Interpretable deep learning translation of GWAS and multi-omics findings to identify pathobiology and drug repurposing in Alzheimer’s disease,” Cell Reports, vol. 41, p. 111717, Nov. 2022.
  * [9] L. Liu, Q. Meng, C. Weng, Q. Lu, T. Wang, and Y. Wen, “Explainable deep transfer learning model for disease risk prediction using high-dimensional genomic data,” PLOS Computational Biology, vol. 18, p. e1010328, July 2022.  Publisher: Public Library of Science.
  * [10] A. Badré and C. Pan, “Explainable multi-task learning improves the parallel estimation of polygenic risk scores for many diseases through shared genetic basis,” PLOS Computational Biology, vol. 19, p. e1011211, July 2023.  Publisher: Public Library of Science.
  * [11] Y. Liu, D. Wang, F. He, J. Wang, T. Joshi, and D. Xu, “Phenotype Prediction and Genome-Wide Association Study Using Deep Convolutional Neural Network of Soybean,” Frontiers in Genetics, vol. 10, 2019.
  * [12] B. Mieth, A. Rozier, J. A. Rodriguez, M. M. C. Höhne, N. Görnitz, and K.-R. Müller, “DeepCOMBI: explainable artificial intelligence for the analysis and discovery in genome-wide association studies,” NAR Genomics and Bioinformatics, vol. 3, p. lqab065, Sept. 2021.
  * [13] P. H. Kassani, F. Lu, Y. Le Guen, M. E. Belloy, and Z. He, “Deep neural networks with controlled variable selection for the identification of putative causal genetic variants,” Nature Machine Intelligence, vol. 4, pp. 761–771, Sept. 2022.  Publisher: Nature Publishing Group.
  * [14] L. Leitsalu, T. Haller, T. Esko, M.-L. Tammesoo, H. Alavere, H. Snieder, M. Perola, P. C. Ng, R. Mägi, L. Milani, K. Fischer, and A. Metspalu, “Cohort Profile: Estonian Biobank of the Estonian Genome Center, University of Tartu,” International Journal of Epidemiology, vol. 44, pp. 1137–1147, Aug. 2015.
  * [15] D. Ho, K. Imai, G. King, and E. A. Stuart, “MatchIt: Nonparametric Preprocessing for Parametric Causal Inference,” Journal of Statistical Software, vol. 42, pp. 1–28, June 2011.
  * [16] M. Buda, A. Maki, and M. A. Mazurowski, “A systematic study of the class imbalance problem in convolutional neural networks,” Neural Networks, vol. 106, pp. 249–259, Oct. 2018.
  * [17] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Köpf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala, “PyTorch: An Imperative Style, High-Performance Deep Learning Library,” Dec. 2019.  arXiv:1912.01703 [cs, stat].
  * [18] K. Simonyan, A. Vedaldi, and A. Zisserman, “Deep Inside Convolutional Networks: Visualising Image Classification Models and Saliency Maps,” Apr. 2014.  arXiv:1312.6034 [cs].
  * [19] M. Sundararajan, A. Taly, and Q. Yan, “Axiomatic attribution for deep networks,” in Proceedings of the 34th International Conference on Machine Learning - Volume 70, ICML’17, (Sydney, NSW, Australia), pp. 3319–3328, JMLR.org, Aug. 2017.
  * [20] N. Kokhlikyan, V. Miglani, M. Martin, E. Wang, B. Alsallakh, J. Reynolds, A. Melnikov, N. Kliushkina, C. Araya, S. Yan, and O. Reblitz-Richardson, “Captum: A unified and generic model interpretability library for PyTorch,” 2020.  _eprint: 2009.07896.
  * [21] L. Breiman, “Random Forests,” Machine Learning, vol. 45, pp. 5–32, Oct. 2001.
  * [22] M. J. Owen, “Genomic insights into schizophrenia,” Royal Society Open Science, vol. 10, p. 230125, Feb. 2023.  Publisher: Royal Society.
  * [23] S. Seabold and J. Perktold, “statsmodels: Econometric and statistical modeling with python,” in 9th Python in Science Conference, 2010.
  * [24] N. Pujol Gualdo, Estonian Biobank Research Team, R. Mägi, and T. Laisk, “Genome-wide association study meta-analysis supports association between MUC1 and ectopic pregnancy,” Human Reproduction, vol. 38, pp. 2516–2525, Dec. 2023.
  * [25] S. Rao, A. Baranova, Y. Yao, J. Wang, and F. Zhang, “Genetic Relationships between Attention-Deficit/Hyperactivity Disorder, Autism Spectrum Disorder, and Intelligence,” Neuropsychobiology, vol. 81, pp. 484–496, June 2022.
  * [26] M. Lam, W. D. Hill, J. W. Trampush, J. Yu, E. Knowles, G. Davies, E. Stahl, L. Huckins, D. C. Liewald, S. Djurovic, I. Melle, K. Sundet, A. Christoforou, I. Reinvang, P. DeRosse, A. J. Lundervold, V. M. Steen, T. Espeseth, K. Räikkönen, E. Widen, A. Palotie, J. G. Eriksson, I. Giegling, B. Konte, A. M. Hartmann, P. Roussos, S. Giakoumaki, K. E. Burdick, A. Payton, W. Ollier, O. Chiba-Falek, D. K. Attix, A. C. Need, E. T. Cirulli, A. N. Voineskos, N. C. Stefanis, D. Avramopoulos, A. Hatzimanolis, D. E. Arking, N. Smyrnis, R. M. Bilder, N. A. Freimer, T. D. Cannon, E. London, R. A. Poldrack, F. W. Sabb, E. Congdon, E. D. Conley, M. A. Scult, D. Dickinson, R. E. Straub, G. Donohoe, D. Morris, A. Corvin, M. Gill, A. R. Hariri, D. R. Weinberger, N. Pendleton, P. Bitsios, D. Rujescu, J. Lahti, S. Le Hellard, M. C. Keller, O. A. Andreassen, I. J. Deary, D. C. Glahn, A. K. Malhotra, and T. Lencz, “Pleiotropic Meta-Analysis of Cognition, Education, and Schizophrenia Differentiates Roles of Early Neurodevelopmental and Adult Synaptic Pathways,” The American Journal of Human Genetics, vol. 105, pp. 334–350, Aug. 2019.
  * [27] V. Trubetskoy, A. F. Pardiñas, T. Qi, et al., “Mapping genomic loci implicates genes and synaptic biology in schizophrenia,” Nature, vol. 604, pp. 502–508, Apr. 2022.  Publisher: Nature Publishing Group.
  * [28] N. Mullins, A. J. Forstner, K. S. O’Connell, et al., “Genome-wide association study of more than 40,000 bipolar disorder cases provides new insights into the underlying biology,” Nature Genetics, vol. 53, pp. 817–829, June 2021.  Publisher: Nature Publishing Group.
  * [29] E. A. Stahl, G. Breen, A. J. Forstner, et al., “Genome-wide association study identifies 30 loci associated with bipolar disorder,” Nature Genetics, vol. 51, pp. 793–803, May 2019.  Publisher: Nature Publishing Group.
  * [30] P. Turley, R. K. Walters, O. Maghzian, A. Okbay, J. J. Lee, M. A. Fontana, T. A. Nguyen-Viet, R. Wedow, M. Zacher, N. A. Furlotte, P. Magnusson, S. Oskarsson, M. Johannesson, P. M. Visscher, D. Laibson, D. Cesarini, B. M. Neale, and D. J. Benjamin, “Multi-trait analysis of genome-wide association summary statistics using MTAG,” Nature Genetics, vol. 50, pp. 229–237, Feb. 2018.  Publisher: Nature Publishing Group.
  * [31] B. M. L. Baselmans, R. Jansen, H. F. Ip, J. van Dongen, A. Abdellaoui, M. P. van de Weijer, Y. Bao, M. Smart, M. Kumari, G. Willemsen, J.-J. Hottenga, D. I. Boomsma, E. J. C. de Geus, M. G. Nivard, and M. Bartels, “Multivariate genome-wide analyses of the well-being spectrum,” Nature Genetics, vol. 51, pp. 445–451, Mar. 2019.  Publisher: Nature Publishing Group.
  * [32] D. F. Levey, M. B. Stein, F. R. Wendt, G. A. Pathak, H. Zhou, M. Aslan, R. Quaden, K. M. Harrington, Y. Z. Nuñez, C. Overstreet, K. Radhakrishnan, G. Sanacora, A. M. McIntosh, J. Shi, S. S. Shringarpure, J. Concato, R. Polimanti, and J. Gelernter, “Bi-ancestral depression GWAS in the Million Veteran Program and meta-analysis in >1.2 million individuals highlight new therapeutic directions,” Nature Neuroscience, vol. 24, pp. 954–963, July 2021.  Publisher: Nature Publishing Group.
  * [33] X. Meng, G. Navoly, O. Giannakopoulou, D. F. Levey, D. Koller, G. A. Pathak, N. Koen, K. Lin, M. J. Adams, M. E. Rentería, Y. Feng, J. M. Gaziano, D. J. Stein, H. J. Zar, M. L. Campbell, D. A. van Heel, B. Trivedi, S. Finer, A. McQuillin, N. Bass, V. K. Chundru, H. C. Martin, Q. Q. Huang, M. Valkovskaya, C.-Y. Chu, S. Kanjira, P.-H. Kuo, H.-C. Chen, S.-J. Tsai, Y.-L. Liu, K. S. Kendler, R. E. Peterson, N. Cai, Y. Fang, S. Sen, L. J. Scott, M. Burmeister, R. J. F. Loos, M. H. Preuss, K. V. Actkins, L. K. Davis, M. Uddin, A. H. Wani, D. E. Wildman, A. E. Aiello, R. J. Ursano, R. C. Kessler, M. Kanai, Y. Okada, S. Sakaue, J. A. Rabinowitz, B. S. Maher, G. Uhl, W. Eaton, C. S. Cruz-Fuentes, G. A. Martinez-Levy, A. I. Campos, I. Y. Millwood, Z. Chen, L. Li, S. Wassertheil-Smoller, Y. Jiang, C. Tian, N. G. Martin, B. L. Mitchell, E. M. Byrne, S. Awasthi, J. R. I. Coleman, S. Ripke, T. Sofer, R. G. Walters, A. M. McIntosh, R. Polimanti, E. C. Dunn, M. B. Stein, J. Gelernter, C. M. Lewis, and K. Kuchenbaecker, “Multi-ancestry genome-wide association study of major depression aids locus discovery, fine mapping, gene prioritization and causal inference,” Nature Genetics, vol. 56, pp. 222–233, Feb. 2024.  Publisher: Nature Publishing Group.
  * [34] T. D. Als, M. I. Kurki, J. Grove, G. Voloudakis, K. Therrien, E. Tasanko, T. T. Nielsen, J. Naamanka, K. Veerapen, D. F. Levey, J. Bendl, J. Bybjerg-Grauholm, B. Zeng, D. Demontis, A. Rosengren, G. Athanasiadis, M. Bækved-Hansen, P. Qvist, G. Bragi Walters, T. Thorgeirsson, H. Stefánsson, K. L. Musliner, V. M. Rajagopal, L. Farajzadeh, J. Thirstrup, B. J. Vilhjálmsson, J. J. McGrath, M. Mattheisen, S. Meier, E. Agerbo, K. Stefánsson, M. Nordentoft, T. Werge, D. M. Hougaard, P. B. Mortensen, M. B. Stein, J. Gelernter, I. Hovatta, P. Roussos, M. J. Daly, O. Mors, A. Palotie, and A. D. Børglum, “Depression pathophysiology, risk prediction of recurrence and comorbid psychiatric disorders using genome-wide analyses,” Nature Medicine, vol. 29, pp. 1832–1844, July 2023.  Publisher: Nature Publishing Group.
  * [35] A. J. Schork, H. Won, V. Appadurai, R. Nudel, M. Gandal, O. Delaneau, M. Revsbech Christiansen, D. M. Hougaard, M. Bækved-Hansen, J. Bybjerg-Grauholm, M. Giørtz Pedersen, E. Agerbo, C. Bøcker Pedersen, B. M. Neale, M. J. Daly, N. R. Wray, M. Nordentoft, O. Mors, A. D. Børglum, P. Bo Mortensen, A. Buil, W. K. Thompson, D. H. Geschwind, and T. Werge, “A genome-wide association study of shared risk across psychiatric disorders implicates gene regulation during fetal neurodevelopment,” Nature Neuroscience, vol. 22, pp. 353–361, Mar. 2019.  Publisher: Nature Publishing Group.
  * [36] J. Kaplanis, N. Akawi, G. Gallone, J. F. McRae, E. Prigmore, C. F. Wright, D. R. Fitzpatrick, H. V. Firth, J. C. Barrett, M. E. Hurles, and o. b. o. t. D. D. D. Study, “Exome-wide assessment of the functional impact and pathogenicity of multinucleotide mutations,” Genome Research, vol. 29, pp. 1047–1056, July 2019.  Company: Cold Spring Harbor Laboratory Press Distributor: Cold Spring Harbor Laboratory Press Institution: Cold Spring Harbor Laboratory Press Label: Cold Spring Harbor Laboratory Press Publisher: Cold Spring Harbor Lab.
  * [37] Y. Chen, W. Li, L. Lv, and W. Yue, “Shared Genetic Determinants of Schizophrenia and Autism Spectrum Disorder Implicate Opposite Risk Patterns: A Genome-Wide Analysis of Common Variants,” Schizophrenia Bulletin, p. sbae044, Apr. 2024.
  * [38] F. Aguet, A. A. Brown, S. E. Castel, et al., “Genetic effects on gene expression across human tissues,” Nature, vol. 550, pp. 204–213, Oct. 2017.  Publisher: Nature Publishing Group.
  * [39] K. Watanabe, E. Taskesen, A. van Bochoven, and D. Posthuma, “Functional mapping and annotation of genetic associations with FUMA,” Nature Communications, vol. 8, p. 1826, Nov. 2017.  Publisher: Nature Publishing Group.
  * [40] E. Candès, Y. Fan, L. Janson, and J. Lv, “Panning for Gold: ‘Model-X’ Knockoffs for High Dimensional Controlled Variable Selection,” Journal of the Royal Statistical Society Series B: Statistical Methodology, vol. 80, pp. 551–577, June 2018.
  * [41] E. A. Boyle, Y. I. Li, and J. K. Pritchard, “An Expanded View of Complex Traits: From Polygenic to Omnigenic,” Cell, vol. 169, pp. 1177–1186, June 2017.