# Spatial Confidence Regions for Combinations of Excursion Sets in Image
Analysis

Thomas Maullin-Sapey1, Armin Schwartzman2,3 and Thomas E. Nichols1

1Big Data Institute, Li Ka Shing Centre for Health Information and Discovery,
Nuffield Department of Population Health, University of Oxford, Oxford, UK;
2Division of Biostatistics, University of California, San Diego, CA, USA;
3Halicioğlu Data Science Institute, University of California, San Diego, USA
TMS<EMAIL_ADDRESS>AS<EMAIL_ADDRESS>TEN:
<EMAIL_ADDRESS>
###### Abstract

The analysis of excursion sets in imaging data is essential to a wide range of
scientific disciplines such as neuroimaging, climatology and cosmology.
Despite growing literature, there is little published concerning the
comparison of processes that have been sampled across the same spatial region
but which reflect different study conditions. Given a set of asymptotically
Gaussian random fields, each corresponding to a sample acquired for a
different study condition, this work aims to provide confidence statements
about the intersection, or union, of the excursion sets across all fields.
Such spatial regions are of natural interest as they directly correspond to
the questions “Where do all random fields exceed a predetermined threshold?”,
or “Where does at least one random field exceed a predetermined threshold?”.
To assess the degree of spatial variability present, we develop a method that
provides, with a desired confidence, subsets and supersets of spatial regions
defined by logical conjunctions (i.e. set intersections) or disjunctions (i.e.
set unions), without any assumption on the dependence between the different
fields. The method is verified by extensive simulations and demonstrated using
a task-fMRI dataset to identify brain regions with activation common to four
variants of a working memory task.

Keywords: Spatial Statistics, Excursion Sets, Confidence Regions, Linear Model

##  1 Introduction

The collection and analysis of imaging data, modelled as a random field
sampled on a spatial domain, is central to a broad range of scientific
disciplines such as neuroimaging, climatology, and cosmology. Often, spatial
data drawn from $n$ observations of a random field are combined to obtain an
estimate, $\hat{\mu}_{n}$, of some spatially varying target function $\mu$.
The target function $\mu$ is defined on a closed spatial domain,
$S\subset\mathbb{R}^{N}$, and maps spatial locations to some variable of
interest in $\mathbb{R}$ (e.g. seismic activity, infrared heat, changes in
blood oxygenation in the brain, etc.). In such applications, it is typically
desirable to ask “At which locations does the target function exceed a certain
value, $c$?”. For example, “Where has a significant change in temperature
occurred?”; “Where do significant heat readings indicate the presence of
celestial bodies?”. Such questions are addressed by the study of excursion
sets, i.e. sets of the form $\\{s\in S:\mu(s)\geq c\\}$.

A wealth of literature has previously focused on documenting the geometric and
topological properties of excursion sets of random functions (c.f. Adler
(1981), Torres (1994), Cao (1999), Azas and Wschebor (2009), Adler and Taylor
(2009)). These properties include, for example, the Euler characteristic (a
measure of topological structure), the Hausdorff dimension (indicating how
fractal the set may be), Lipschitz Killing curvatures (describing high-
dimensional volumes and areas) and Betti Numbers (describing how many
stationary points appear above the threshold, $c$) (Worsley (1996), Adler
(1977), Adler et al. (2017), Pranav et al. (2019)). Much of this work, though,
is limited to homogeneous stationary processes, i.e. those with zero mean and
a spatial covariance structure which is dependent upon only distance. Only the
most recent literature, including the present paper, concerns processes with a
non-zero mean function and a potentially heterogeneous covariance function.

In addition, at present there is little published concerning how excursion
sets may be compared to one another. In many applications, it is typical for a
researcher to collect data from two or more study conditions and contrast the
results to draw some meaningful inference (e.g. temperature changes in winter
vs summer, heat measurements gathered using different imaging modalities,
brain activity in healthy subjects across a range of tasks, etc.). Such
settings are akin to having $M$ spatially aligned, potentially correlated,
estimates,
$\hat{\mu}_{n}^{1},\hat{\mu}_{n}^{2},...\hat{\mu}_{n}^{M}:S\rightarrow\mathbb{R}$,
of $M$ target functions, $\mu^{1},\mu^{2},...\mu^{M}$.

Often, pertinent questions that arise in such settings may be formally
expressed using logical statements which involve conjunctions (logical
‘and’s), negations (logical ‘not’s) and disjunctions (logical ‘or’s). For
example, a climatologist may ask where significant changes in temperature were
observed during either winter or summer (i.e. “Where does $\mu^{1}$ or
$\mu^{2}$ exceed $c$?’). Alternatively, a neuroscientist may ask “At which
locations in the brain did subjects exhibit signs of cognitive behaviour
during one task and not another?” (e.g. “At which locations does $\mu^{1}$,
and not $\mu^{2}$, exceed $c$?”).

One approach to answering such questions is to employ null-hypothesis testing.
A hypothesis test of a disjunction of nulls to assess evidence for a
conjunction of alternative hypotheses can be conducted with an ‘intersection-
union’ test. An intersection-union test rejects at level $\alpha$ when all of
the individual nulls are rejected at level $\alpha$. This approach has proven
particularly useful for tests of bioequivalence (Berger and Hsu (1996), Berger
(1997)), but does not consider the spatial aspect that is our focus here.

It is common practice for researchers to compare images corresponding to
different study conditions via rudimentary visual inspection. (In fMRI, just a
few recent examples of qualitative assessment of ‘overlap’ and ‘differences’
are Zhang et al. (2021), Dijkstra et al. (2021), Ferreira et al. (2021)). Such
a subjective practice may lead to biased or misleading results. In recent
years, much literature has been published on the theory of confidence regions,
which focuses on quantifying spatial uncertainty by providing, with a fixed
confidence, sub- and super-sets bounding the excursion set of a single target
function (Sommerfeld et al. (2018), Bowring et al. (2019), Bowring et al.
(2021)). However, the theory of confidence regions does not currently offer a
means for investigating logical conjunctions and disjunctions of exceedance
statements (i.e. intersections and unions of excursion sets). This work
addresses this issue by first proposing a theory of spatial confidence regions
for ‘conjunction inference’ (i.e. inference for logical statements involving
conjunctions) in image analysis. Following this, the proposed method is
extended to allow the investigation of statements containing disjunctions and
negations.

Figure 1:  Illustration of the set $\mathcal{F}_{c}$ for a setting in which
$M=2$. Shown in $(a)$ are two spatially-varying target functions, $\mu^{1}(s)$
(blue) and $\mu^{2}(s)$ (green), overlaid and thresholded at the level $c$. In
$(b)$, the regions at which $\mu^{1}(s)\geq c$ and $\mu^{2}(s)\geq c$ are
displayed as red and yellow circles, respectively, with the intersection set,
$\mathcal{F}_{c}$, illustrated in purple. In $(c)$, potential confidence
regions, $\hat{\mathcal{F}}_{c}^{+}$ (grey) and $\hat{\mathcal{F}}_{c}^{-}$
(orange), for $\mathcal{F}_{c}$ are illustrated.

Formally, our work primarily focuses on the intersection of $M$ excursion sets
(i.e. the region over which the conjunction statement “$\mu^{1}(s)\geq c$ and
$\mu^{2}(s)\geq c$ and … and $\mu^{M}(s)\geq c$” holds). We denote
$\mathcal{M}=\\{1,...,M\\}$, and for each $i\in\mathcal{M}$ (i.e. for each
study condition) define the $i^{th}$ true excursion set and $i^{th}$ estimated
excursion set as:

$\mathcal{A}^{i}_{c}=\\{s\in S:\mu^{i}(s)\geq c\\}\quad\text{ and
}\quad\hat{\mathcal{A}}^{i}_{c}=\\{s\in S:\hat{\mu}_{n}^{i}(s)\geq c\\},$

respectively. Next, we define $\mathcal{F}_{c}$ to be the intersection of the
sets $\\{\mathcal{A}_{c}^{i}\\}_{i\in\mathcal{M}}$, and
$\hat{\mathcal{F}}_{c}$ to be the intersection of the sets
$\\{\hat{\mathcal{A}}_{c}^{i}\\}_{i\in\mathcal{M}}$. In other words, the
spatial region we are interested in, and our estimate of that region, are
given by:

$\mathcal{F}_{c}=\bigg{\\{}s\in S:\min_{i\in\mathcal{M}}\mu^{i}(s)\geq
c\bigg{\\}}\quad\text{and}\quad\hat{\mathcal{F}}_{c}=\bigg{\\{}s\in
S:\min_{i\in\mathcal{M}}\hat{\mu}_{n}^{i}(s)\geq c\bigg{\\}},$

respectively. An illustration of $\mathcal{F}_{c}$ in a setting in which $M=2$
is provided by Fig. 1.

We note here that the above construction can be adjusted to allow $c$ to vary
across study conditions (i.e. “$\mu^{1}(s)\geq c_{1}$ and … and
$\mu^{M}(s)\geq c_{M}$” for $c_{1},...,c_{M}\in\mathbb{R}$). Such an
adjustment would involve simply translating each field upwards or downwards
(e.g. substituting $\\{\mu^{i}\\}_{i\in\mathcal{M}}$ for
$\\{\mu^{i}-c+c_{i}\\}_{i\in\mathcal{M}}$). For ease in the proceeding text,
we shall assume $c$ is equal for all study conditions.

Our goal is to obtain confidence regions $\hat{\mathcal{F}}^{+}_{c}$ and
$\hat{\mathcal{F}}^{-}_{c}$ such that the below probability holds
asymptotically;

$\mathbb{P}\bigg{[}\hat{\mathcal{F}}^{+}_{c}\subseteq\mathcal{F}_{c}\subseteq\hat{\mathcal{F}}^{-}_{c}\bigg{]}=1-\alpha$
(1)

for a predefined tolerance level $\alpha$ ($\alpha=0.05$, for example). To
generate such regions, we build on the theory of Sommerfeld et al. (2018), to
show that under an appropriate definition of $\hat{\mathcal{F}}^{+}_{c}$ and
$\hat{\mathcal{F}}^{-}_{c}$, the above probability may be approximated using
quantiles of a well-defined random variable. Using a wild $t-$bootstrapping
procedure, we will then demonstrate that the relationship between this random
variable and the above probability can be used to generate
$\hat{\mathcal{F}}^{+}_{c}$ and $\hat{\mathcal{F}}^{-}_{c}$ for any desired
value of $\alpha$. Unlike much of the previous literature on random excursion
sets, in this work no assumption is made on the processes’ means or spatial
covariance structure, and we do not make any assumption of between-‘study
condition’ independence (i.e. $\\{\hat{\mu}_{n}^{i}\\}_{i\in\mathcal{M}}$ may
be correlated with one another).

In the following sections, we first describe the notation and assumptions upon
which our theory relies. Following this, we provide the central theoretical
result of this work, which relates Equation (1) to an exceedance statement for
a well-defined random variable. Next, via the use of the wild $t-$bootstrap,
we detail how this result may be employed to obtain confidence regions for
conjunction inference in the setting of linear regression modelling. Finally,
we validate our results with simulations and a real data example based on fMRI
data taken from the Human Connectome Project (Van Essen et al. (2013)). Proofs
of the theory presented in this work, alongside further illustration and
extensive simulation results, are provided as supplementary material.

##  2 Confidence Regions for Excursion Set Combinations

###  2.1 Notation

In the following sections, we shall need notation to describe the numerous
possible low dimensional sub-manifolds which can arise from intersecting the
boundaries of $M$ excursion sets. To do so, we first denote the power set of a
finite set, $A$, as $\mathcal{P}(A)$, and define $\mathcal{P}^{+}(A)$ as
$\mathcal{P}^{+}(A)=\mathcal{P}(A)\setminus\\{\emptyset\\}$. For each
$i\in\mathcal{M}$, we define $\partial\mathcal{A}_{c}^{i}$ as the level set
$\\{s\in S:\mu^{i}(s)=c\\}$. Similarly, we define $\partial\mathcal{F}_{c}$ as
the level set $\\{s\in S:\min_{i\in\mathcal{M}}\mu^{i}(s)=c\\}$. For a spatial
set, $B\subseteq S$, its complement in $S$ is denoted $B^{c}=S\setminus B$,
its topological closure is denoted $\overline{B}$, its interior is denoted
$B^{\circ}$ and its boundary is denoted $\partial B$. In addition, for closed
$B$, the notation $B^{1}$ and $B^{-1}$ shall represent $B$ and
$\overline{B^{c}}$, respectively. We note here that the notation
$\partial\mathcal{A}^{i}_{c}$ and $\partial B$ may conflict with one another.
This conflict is resolved, however, by the assumptions of Section 2.2 which
ensure that the boundaries of $\\{\mathcal{A}_{c}^{i}\\}_{i\in\mathcal{M}}$
and $\mathcal{F}_{c}$ are equal to the level sets
$\\{\partial\mathcal{A}_{c}^{i}\\}_{i\in\mathcal{M}}$ and
$\partial\mathcal{F}_{c}$ (at least locally, in the vicinity of
$\partial\mathcal{F}_{c}$, which is sufficient for our purposes).

We can now partition the level set $\partial\mathcal{F}_{c}$ into sub-
manifolds using the set of possible intersections of the $M$ excursion set
boundaries. For all $\alpha\in\mathcal{P}^{+}(\mathcal{M})$, we define the
boundary segment $\partial^{\alpha}\mathcal{F}_{c}$ as follows;

$\partial^{\alpha}\mathcal{F}_{c}=\partial\mathcal{F}_{c}\cap\bigg{(}\bigcap_{i\in\alpha}\partial\mathcal{A}_{c}^{i}\bigg{)}\cap\bigg{(}\bigcap_{j\in\mathcal{M}\setminus\alpha}(\partial\mathcal{A}_{c}^{j})^{c}\bigg{)}.$

We note it is possible that, for some $\alpha\in\mathcal{P}^{+}(\mathcal{M})$,
the boundary segment $\partial^{\alpha}\mathcal{F}_{c}$ is empty (in fact,
this is often the case in practice). This notation is crucial to the statement
of Theorem 2.1 and is illustrated by Fig. 2.

Figure 2:  Annotated images of ‘overlapping’ excursion sets (top) with
corresponding illustrations of boundary partitions (bottom) for settings in
which the number of study conditions, $M$, is equal to $2$ (left) or $3$
(right).

###  2.2 Assumptions

In this section, we outline and discuss the assumptions upon which our theory
relies. Additional discussion of why each assumption is required, alongside
illustration, is given in Supplementary Theory Section S2.

######  Assumption 2.2.1.

For each $i\in\mathcal{M}$, there exists a bounded function
$\sigma^{i}:S\rightarrow\mathbb{R}^{+}$ and positive sequence
$\tau_{n}\rightarrow 0$, such that the below Central Limit Theorem (CLT)
holds:

$\bigg{\\{}\frac{\hat{\mu}^{i}_{n}(s)-\mu^{i}(s)}{\tau_{n}\sigma^{i}(s)}\bigg{\\}}_{s\in
S,i\in\mathcal{M}}\xrightarrow{d}\\{G^{i}(s)\\}_{s\in S,i\in\mathcal{M}}$

where $\\{G^{i}(s)\\}_{s\in S,i\in\mathcal{M}}$ is a well-defined multi-
variate random field with continuous sample paths in $S$ and $\xrightarrow{d}$
represents convergence in distribution.

######  Remark.

This assumption introduces the notation $\sigma^{i}(s)$ and $\tau_{n}$. In
typical applications, $\sigma^{i}(s)$ represents standard deviation across
observations (for study condition $i$, at spatial location $s$) and $\tau_{n}$
is a decreasing function of sample size. For conciseness in the following
text, we now define $\\{g^{i}\\}_{i\in\mathcal{M}}$ and
$\\{\hat{g}_{n}^{i}\\}_{i\in\mathcal{M}}$ as;

$g^{i}(s)=\frac{\mu^{i}(s)-c}{\sigma^{i}(s)}\quad\text{ and
}\quad\hat{g}_{n}^{i}(s)=\frac{\hat{\mu}_{n}^{i}(s)-c}{\sigma^{i}(s)}.$

######  Assumption 2.2.2.

We assume that:

  1. (a)

The functions $\\{g^{i}\\}_{i\in\mathcal{M}}$ are continuous on $S$.

  2. (b)

For $n$ large enough, the functions $\\{\hat{g}^{i}_{n}\\}_{i\in\mathcal{M}}$
are continuous on $S$.

######  Remark.

In practice, assumptions of this form are already common in the imaging
literature, as in many applications $\\{\mu^{i}\\}_{i\in\mathcal{M}}$,
$\\{\hat{\mu}_{n}^{i}\\}_{i\in\mathcal{M}}$ and
$\\{\sigma^{i}\\}_{i\in\mathcal{M}}$ are assumed to be continuous across
space. For further remarks, see Supplementary Theory Section S2.2.

######  Assumption 2.2.3.

We assume that:

  1. (a)

Every open ball around a point on the boundary $\partial\mathcal{F}_{c}$ has a
non-empty intersection with $(\mathcal{F}_{c})^{\circ}$. In addition, for all
$\alpha\in\mathcal{P}^{+}(\mathcal{M})$ every open ball around a point in
$\partial^{\alpha}\mathcal{F}_{c}$ has a non-empty intersection with the set:

$\mathcal{J}_{c}^{\alpha}:=\bigg{(}\bigcap_{i\in\alpha}(\mathcal{A}^{i}_{c})^{c}\bigg{)}\cap\bigg{(}\bigcap_{j\in\mathcal{M}\setminus\alpha}(\mathcal{A}^{j}_{c})^{\circ}\bigg{)}$

where, if $\mathcal{M}\setminus\alpha$ is empty, the last term is defined to
be equal to $S$.

  2. (b)

There is an open neighbourhood of $\partial\mathcal{F}_{c}$ over which the
functions $\\{g^{i}\\}_{i\in\mathcal{M}}$ and, for $n$ large enough,
$\\{\hat{g}^{i}_{n}\\}_{i\in\mathcal{M}}$ are $C^{1}$ with finite, non-zero,
gradients.

  3. (c)

For every point $s\in\partial\mathcal{F}_{c}$ and
$\alpha\in\mathcal{P}^{+}(\mathcal{M})$, if
$s\in\partial\big{(}\bigcap_{i\in\alpha}\mathcal{A}_{c}^{i}\big{)}$ then the
set $\partial\big{(}\bigcap_{i\in\alpha}\mathcal{A}_{c}^{i}\big{)}$ partitions
every sufficiently small open ball around $s$ into exactly two components,
each of which is path connected.

######  Remark.

The above statements ensure that $\mathcal{F}_{c}$ is non-empty and has a well
defined boundary which is equal to the level set
$\partial\mathcal{F}_{c}=\\{s\in S:\min_{i\in\mathcal{M}}g^{i}(s)=0\\}$. All
three statements ensure that no changes in topology occur at the level $c$ by
requiring that $\partial\mathcal{F}_{c}$ does not contain plateaus (spatial
regions of non-zero measure over which $\min_{i\in\mathcal{M}}g^{i}(s)=0$
exactly), local minima or maxima. In addition, each statement handles
pathological cases which the other statements do not. For brevity, we do not
list such cases here but instead provide further discussion in Supplementary
Theory Section S2.3.

Each of the examples presented in Fig. 2 satisfy Assumptions 2.2.2 and 2.2.3
and highlight several features worthy of note. Firstly, we do not assume
$\partial\mathcal{F}_{c}$ is a $C^{1}$ curve, but rather piece-wise $C^{1}$
(c.f. examples $(d)$ and $(e)$). Secondly, the assumptions allow for the
possibility that the excursion sets could be nested within one another (c.f.
examples $(c)$ and $(f)$). Thirdly, the assumptions permit the excursion sets
to share common boundaries (see examples $(b)$, $(c)$ and $(f)$). The last two
observations are of practical relevance as, in many applications, it may be
expected that the target functions for different study conditions exhibit a
strong similarity. For instance, in the neuroimaging example, it may be
expected that some anatomical regions of the brain display evidence of
activation for all of the study conditions. We note here that there is no
requirement that $\mathcal{F}_{c}$ be (globally) path-connected;
$\mathcal{F}_{c}$ may consist of several disconnected components without
violating the assumptions of this section.

###  2.3 Theory

In this section, we present the central result of this work, Theorem 2.1,
alongside two corollaries, Corollary 2.1 and Corollary 2.2. To do so, we now
define the nested sets $\hat{\mathcal{F}}^{-}_{c}$ and
$\hat{\mathcal{F}}^{+}_{c}$ by thresholding the statistic
$\tau_{n}^{-1}\min_{i\in\mathcal{M}}\hat{g}^{i}_{n}(s)$:

$\hat{\mathcal{F}}^{-}_{c}:=\hat{\mathcal{F}}^{-}_{c}(a)=\\{s\in
S:\tau_{n}^{-1}\min_{i\in\mathcal{M}}\hat{g}^{i}_{n}(s)\geq-a\\},$
$\hat{\mathcal{F}}^{+}_{c}:=\hat{\mathcal{F}}^{+}_{c}(a)=\\{s\in
S:\tau_{n}^{-1}\min_{i\in\mathcal{M}}\hat{g}^{i}_{n}(s)\geq+a\\},$

using some constant $a\in\mathbb{R}^{+}$. To find an appropriate value of $a$,
such that Equation (1) holds for a desired tolerance level (e.g.
$\alpha=0.05$), Theorem 2.1 is required.

######  Theorem 2.1.

Under the assumptions of Section 2.2, the below holds:

$\lim_{n\rightarrow\infty}\mathbb{P}[\hat{\mathcal{F}}^{+}_{c}\subseteq\mathcal{F}_{c}\subseteq\hat{\mathcal{F}}^{-}_{c}]=\mathbb{P}[H\leq
a],$ (2)

where the variable $H$ is defined as follows;

$H=\max_{\alpha\in\mathcal{P}^{+}(\mathcal{M})}\bigg{(}\sup_{s\in\partial^{\alpha}\mathcal{F}_{c}}\big{|}\min_{i\in\alpha}(G^{i}(s))\big{|}\bigg{)}.$

Conceptually, the event $\\{H>a\\}$ may be thought of as the situation in
which, at some point, $s$, inside some boundary segment,
$\partial^{\alpha}\mathcal{F}_{c}$, a large value was observed for the
statistic $|\min_{i\in\alpha}(G^{i}(s))|$. Therefore, Theorem 2.1 tells us the
probability that the statement
$\hat{\mathcal{F}}^{+}_{c}\subseteq\mathcal{F}_{c}\subseteq\hat{\mathcal{F}}^{-}_{c}$
is violated is directly determined by such points. If we can find a value of
$a$ such that $\mathbb{P}[H\leq a]=1-\alpha$ then asymptotically, the
inclusion probability,
$\mathbb{P}[\hat{\mathcal{F}}^{+}_{c}\subseteq\mathcal{F}_{c}\subseteq\hat{\mathcal{F}}^{-}_{c}]$,
will also equal $1-\alpha$. In other words, if the $(1-\alpha)^{th}$ quantile
of the distribution of $H$ is known, then confidence regions
$\hat{\mathcal{F}}^{-}_{c}$ and $\hat{\mathcal{F}}^{+}_{c}$ may be constructed
which satisfy Equation (1). We shall take up the task of evaluating the
quantiles of $H$ in Section 3.

Theorem 2.1 is key to generating confidence regions for conjunction inference
(intersections of excursion sets). However, Theorem 2.1 may also be extended
to allow investigation of other logical statements involving negations or
disjunctions. Below we provide two corollaries, each of which is illustrated
via a worked example. Corollary 2.1 extends Theorem 2.1 to account for logical
negations, whilst Corollary 2.2 provides an equivalent statement for inference
performed upon logical disjunctions (unions of excursion sets).

######  Example 2.1.

Suppose, given two target functions, $\mu^{1}$ and $\mu^{2}$, we were
interested in the region defined by the logical conjunction ‘$\mu^{1}\geq c$
and $\mu^{2}\leq c$’ (i.e. “Where did a variable of interest exceed a
threshold under one study condition and not under another?”). In the notation
of Section 2.1, this region is readily seen to be given by
$(\mathcal{A}_{c}^{1})^{1}\cap(\mathcal{A}_{c}^{2})^{-1}$.

To apply the result of Theorem 2.1, we first define $\tilde{\mu}^{1}=\mu^{1}$
and $\tilde{\mu}^{2}=2c-\mu^{2}$ and note that the logical statement
‘$\mu^{1}\geq c$ and $\mu^{2}\leq c$’ is identical to the statement
‘$\tilde{\mu}^{1}\geq c$ and $\tilde{\mu}^{2}\geq c$’. By using the latter
statement, Theorem 2.1 may now be applied to obtain the desired confidence
regions (assuming that an appropriate mechanism exists for evaluating the
quantiles of $H$). However, during this process, the limiting variables
$\\{G^{i}\\}_{i\in\\{1,2\\}}$, which correspond to the target functions
$\\{\mu^{i}\\}_{i\in\\{1,2\\}}$, must be replaced by variables corresponding
to the functions $\\{\tilde{\mu}^{i}\\}_{i\in\\{1,2\\}}$,
$\\{\tilde{G}^{i}\\}_{i\in\\{1,2\\}}$. By noting the definitions of
$\\{\tilde{\mu}^{i}\\}_{i\in\\{1,2\\}}$ and
$\\{\tilde{G}^{i}\\}_{i\in\\{1,2\\}}$, it can be seen that
$\tilde{G}^{1}=G^{1}$ and $\tilde{G}^{2}=-G^{2}$.

In summary, a procedure for generating confidence regions corresponding to the
logical conjunction ‘$\mu^{1}\geq c$ and $\mu^{2}\leq c$’ would be identical
to that previously described, except that $G^{2}$ would be substituted for
$-G^{2}$ and $\partial\mathcal{F}_{c}$ would be substituted for
$\partial((\mathcal{A}^{1}_{c})^{1}\cap(\mathcal{A}^{2}_{c})^{-1})$. In
general, this example may be extended to allow for inference to be performed
upon arbitrary conjunctions of statements and negated statements. To achieve
this, the definitions of $\mathcal{F}_{c}$, $\hat{\mathcal{F}}^{-}_{c}$ and
$\hat{\mathcal{F}}^{+}_{c}$ must be extended as follows.

######  Corollary 2.1.

Let $\\{\delta_{i}\\}_{i\in\mathcal{M}}$ be an arbitrary sequence of integers
in $\\{-1,1\\}$ and define $\mathcal{F}_{c}$ as
$\mathcal{F}_{c}=\bigcap_{i\in\mathcal{M}}(\mathcal{A}_{c}^{i})^{\delta_{i}}$.
Similarly, extend the definition of the sets $\hat{\mathcal{F}}^{+}_{c}$ and
$\hat{\mathcal{F}}^{-}_{c}$ to:

$\hat{\mathcal{F}}^{+}_{c}=\bigcap_{i\in\mathcal{M}}\\{s\in
S:\tau_{n}^{-1}\hat{g}^{i}_{n}(s)\geq+a\\}^{\delta_{i}}\quad\text{and}\quad\hat{\mathcal{F}}^{-}_{c}=\bigcap_{i\in\mathcal{M}}\\{s\in
S:\tau_{n}^{-1}\hat{g}^{i}_{n}(s)\geq-a\\}^{\delta_{i}},$

respectively and, for each $\alpha\in\mathcal{P}^{+}(\mathcal{M})$, define
$\partial^{\alpha}\mathcal{F}_{c}$ as before. If $\mathcal{F}_{c}$ satisfies
the assumptions of Section 2.2, then the below statement holds asymptotically:

$\lim_{n\rightarrow\infty}\mathbb{P}[\hat{\mathcal{F}}^{+}_{c}\subseteq\mathcal{F}_{c}\subseteq\hat{\mathcal{F}}^{-}_{c}]=\mathbb{P}\bigg{[}\max_{\alpha\in\mathcal{P}^{+}(\mathcal{M})}\bigg{(}\sup_{s\in\partial^{\alpha}\mathcal{F}_{c}}\big{|}\min_{i\in\alpha}(\delta_{i}G^{i}(s))\big{|}\bigg{)}\leq
a\bigg{]}$

######  Example 2.2.

Suppose, given two target functions, $\mu^{1}$ and $\mu^{2}$, interest lay in
the logical disjunction of ‘$\mu^{1}\geq c$ or $\mu^{2}\geq c$’ (i.e. “Where
did at least one of the variables of interest exceed the threshold?”). The
region of space over which this statement holds is given by
$\mathcal{A}_{c}^{1}\cup\mathcal{A}_{c}^{2}$, and shall be denoted
$\mathcal{G}_{c}$.

To generate confidence regions for $\mathcal{G}_{c}$, we first assume that
$\mathcal{G}_{c}$ satisfies Assumptions 2.2.1-2.2.3. By Assumptions 2.2.2 and
2.2.3 and De Morgan’s law, it can be seen that
$\mathcal{G}_{c}=((\mathcal{A}_{c}^{1})^{-1}\cap(\mathcal{A}_{c}^{2})^{-1})^{-1}$.
It must be noted that Assumptions 2.2.2 and 2.2.3 are necessary for this
statement to hold due to the fact that
$(\mathcal{A}_{c}^{i})^{-1}=\overline{(\mathcal{A}_{c}^{i})^{c}}$ and not
$(\mathcal{A}_{c}^{i})^{c}$. By employing Corollary 2.1, for a fixed value of
$\alpha$, confidence regions, $\hat{\mathcal{F}}^{-}_{c}$ and
$\hat{\mathcal{F}}^{+}_{c}$, may be obtained which satisfy the following;

$\lim_{n\rightarrow\infty}\mathbb{P}[(\hat{\mathcal{F}}^{-}_{c})^{c}\subseteq((\mathcal{A}_{c}^{1})^{-1}\cap(\mathcal{A}_{c}^{2})^{-1})^{c}\subseteq(\hat{\mathcal{F}}^{+}_{c})^{c}]=1-\alpha.$

By taking the closure of each of the sets in the above probability, and again
via careful consideration of Assumptions 2.2.2 and 2.2.3, the below is
obtained;

$\lim_{n\rightarrow\infty}\mathbb{P}[(\hat{\mathcal{F}}^{-}_{c})^{-1}\subseteq\mathcal{G}_{c}\subseteq(\hat{\mathcal{F}}^{+}_{c})^{-1}]=1-\alpha.$

In other words, the sets $(\hat{\mathcal{F}}^{-}_{c})^{-1}$ and
$(\hat{\mathcal{F}}^{+}_{c})^{-1}$ serve as confidence regions for disjunction
inference on $\mu^{1}$ and $\mu^{2}$. Note that, as in the previous example,
great attention must paid to the signs which appear in front of the variables
$\\{G^{i}\\}_{i\in\mathcal{M}}$ in the definition of $H$. As this example
began by generating confidence regions for
$(\mathcal{A}_{c}^{1})^{-1}\cap(\mathcal{A}_{c}^{2})^{-1}$, it can be seen
that the variables $G^{1}$ and $G^{2}$ must be appended with negative signs
(i.e. when applying the result of Corollary 2.1, both $\delta_{1}$ and
$\delta_{2}$ were set to $-1$). This example may be expanded upon to obtain
similar results for larger values of $M$, as shown by Corollary 2.2.

######  Corollary 2.2.

Let $\\{\delta_{i}\\}_{i\in\mathcal{M}}$ be an arbitrary sequence of integers
in $\\{-1,1\\}$ and define $\mathcal{G}_{c}$ as
$\mathcal{G}_{c}=\bigcup_{i\in\mathcal{M}}(\mathcal{A}_{c}^{i})^{\delta_{i}}$.
Define the sets $\hat{\mathcal{G}}^{+}_{c}$ and $\hat{\mathcal{G}}^{-}_{c}$
as:

$\hat{\mathcal{G}}^{+}_{c}=\bigg{(}\bigcap_{i\in\mathcal{M}}\\{s\in
S:\tau_{n}^{-1}\hat{g}^{i}_{n}(s)\geq-a\\}^{\delta_{i}}\bigg{)}^{-1}\quad\text{and}\quad\hat{\mathcal{G}}^{-}_{c}=\bigg{(}\bigcap_{i\in\mathcal{M}}\\{s\in
S:\tau_{n}^{-1}\hat{g}^{i}_{n}(s)\geq+a\\}^{\delta_{i}}\bigg{)}^{-1},$

respectively and, for each $\alpha\in\mathcal{P}^{+}(\mathcal{M})$, define
$\partial^{\alpha}\mathcal{G}_{c}$ analogously to
$\partial^{\alpha}\mathcal{F}_{c}$ in the previous sections. If
$\mathcal{G}_{c}$ satisfies the assumptions of Section 2.2, then the below
statement holds asymptotically:

$\lim_{n\rightarrow\infty}\mathbb{P}[\hat{\mathcal{G}}^{+}_{c}\subseteq\mathcal{G}_{c}\subseteq\hat{\mathcal{G}}^{-}_{c}]=\mathbb{P}\bigg{[}\max_{\alpha\in\mathcal{P}^{+}(\mathcal{M})}\bigg{(}\sup_{s\in\partial^{\alpha}\mathcal{G}_{c}}\big{|}\min_{i\in\alpha}(-\delta_{i}G^{i}(s))\big{|}\bigg{)}\leq
a\bigg{]}.$

##  3 Spatial Conjunction Inference for the Linear Model

In this section, we build upon the work of Sommerfeld et al. (2018) and
Bowring et al. (2019), in which methods for generating confidence regions were
presented for a single target function, $\mu(s)$, derived from a linear
regression model. In Section 3.1 we formally describe the spatially-varying
Linear Model (LM) for $M$ study conditions and, in Section 3.2, we explain how
the wild $t$-bootstrap may be applied to estimate quantiles of the variable
$H$. Further implementation details, for when data are sampled from a discrete
lattice rather than across a continuous space, alongside pseudocode, are
provided in Supplementary Results Section S$2$. This section focuses solely on
conjunction inference. However, the methods described here may equally be
applied to perform inference on other logical statements via the use of the
corollaries and examples presented in Section 2.3.

###  3.1 Model Specification

For each $i\in\mathcal{M}$ (i.e. for each study condition), we assume that the
data follow a spatially-varying Linear Model (LM) defined at location $s\in S$
as:

$Y^{i}(s)=X^{i}\beta^{i}(s)+\epsilon^{i}(s),\quad\epsilon^{i}(s)\sim
N(0,\Sigma^{i}(s)).$ (3)

The known quantities in the model are; the $(n\times 1)$ vector of responses,
$Y^{i}(s)$, and the $(n\times p)$ design matrix, $X^{i}$. The unknown model
parameters are; the $(p\times 1)$ vector of regression coefficients,
$\beta^{i}(s)$, and the $(n\times n)$ random error covariance matrix,
$\Sigma^{i}(s)$. For each spatial location, $s\in S$, the Generalized Least
Squares (GLS) estimator of the parameter vector $\beta^{i}(s)$,
$\hat{\beta}^{i}(s)$, is given by:

$\hat{\beta}^{i}(s)=(X^{i^{\prime}}\Sigma^{i}(s)^{-1}X^{i})^{-1}X^{i^{\prime}}\Sigma^{i}(s)^{-1}Y^{i}(s).$

In this context, under the $i^{th}$ study condition, at each spatial location
$s$, our interest lies in assessing whether linear relationships hold between
the elements of the parameter vector $\beta^{i}(s)$. Such relationships may be
expressed using a contrast vector $L^{i}$ of dimension $(p\times 1)$. In the
notation of the previous sections, the target and estimator functions for the
spatially-varying LM are given as $\mu^{i}(s)=L^{i^{\prime}}\beta^{i}(s)$ and
$\hat{\mu}^{i}_{n}(s)=L^{i^{\prime}}\hat{\beta}^{i}(s)$ respectively, and
$\tau_{n}^{-1}\sigma^{i}(s)$ is the contrast standard error, given as:

$\tau_{n}^{-1}\sigma^{i}(s)=\big{[}L^{i^{\prime}}(X^{i^{\prime}}\Sigma^{i}(s)^{-1}X^{i})^{-1}L^{i}\big{]}^{\frac{1}{2}}.$

Whilst in the above notation we have presented
$\\{\hat{\beta}^{i}(s)\\}_{i\in\mathcal{M}}$ as being estimated separately
using $M$ different models, we note that nothing in our theory prevents the
same model being used across all $M$ study conditions. For example, it is
possible for the model matrices, $\\{Y^{i}(s)\\}_{i\in\mathcal{M}}$ and
$\\{X^{i}\\}_{i\in\mathcal{M}}$, and parameter vector,
$\\{\beta^{i}(s)\\}_{i\in\mathcal{M}}$, to be equal for all $i\in\mathcal{M}$,
with the only distinction between study conditions being represented by the
contrast vector $L^{i}$. This observation is noteworthy as, in many
applications, it is common for researchers to compile all study conditions
into a single LM to obtain a heightened statistical power.

To apply Theorem 2.1 to the above LM, we now assume that Assumptions
2.2.1-2.2.3 hold. Such assumptions are commonly satisfied by standard
conditions placed on the continuity and boundedness of the increments and
moments of $\\{\beta^{i}(s)\\}_{s\in S}$ and $\\{\epsilon^{i}(s)\\}_{s\in S}$
(see Sommerfeld et al. (2018) for further detail). We note here that the
covariance matrix, $\Sigma^{i}(s)$, is usually unknown in practice, meaning
the function $\tau_{n}^{-1}\sigma^{i}(s)$ cannot be calculated. As
demonstrated in Sommerfeld et al. (2018), however, $\Sigma^{i}(s)$ can be
replaced by any consistent estimator $\hat{\Sigma}^{i}(s)$ and the assumptions
given in Section 2.2 are still satisfied. Such estimation is common in the
statistics literature and is typically achieved by assuming that some fixed
correlation structure applies to $\Sigma^{i}(s)$ (e.g. diagonal independence,
auto-regressive, etc.) so that the number of independent variance parameters
which must be estimated is small.

###  3.2 The Wild $t-$bootstrap

In practice, the distribution of $H$ is unknown and must be estimated. In this
section, for the model specification described in Section 3.1, we employ a
wild $t-$bootstrap resampling scheme to obtain the quantile, $a$, of $H$ which
satisfies $\mathbb{P}[H\leq a]=1-\alpha$.

To do so, we first define the $(n\times 1)$-dimensional residual vector,
$R^{i}$, for the LM of the $i^{th}$ study condition ($i\in\mathcal{M}$) as:

$[R^{i}_{1}(s),...,R^{i}_{n}(s)]^{\prime}=R^{i}(s)=[\hat{\Sigma}^{i}(s)]^{-\frac{1}{2}}(Y^{i}(s)-X^{i}\hat{\beta}^{i}(s))$

The wild $t-$bootstrap proceeds by, in each bootstrap instance, generating $n$
i.i.d. Rademacher random variables, $\\{r_{1},...,r_{n}\\}$, (random variables
which take the values $-1$ and $+1$ with equal probability) independently of
the data. For the $i^{th}$ study condition, at spatial location $s$, a new
bootstrap sample is then obtained by multiplying the elements of the residual
vector by the Rademacher variables (i.e. the new sample is given by
$\\{r_{1}R^{i}_{1}(s),...,r_{n}R^{i}_{n}(s)\\}$). Once the boostrap sample has
been generated, the standard deviation of the bootstrap sample,
$\hat{\sigma}^{*,i}(s)$, is calculated. A bootstrap instance of the variable
$G^{i}$, $\tilde{G}^{i}$, is now generated as follows;

$\tilde{G}^{i}(s)=n^{-\frac{1}{2}}\sum_{l=1}^{n}\frac{r_{l}R^{i}_{l}(s)}{\hat{\sigma}^{*,i}(s)}.$

A bootstrap instance of the variable $H$, $\tilde{H}$, may then be computed
using the bootstrap variables $\\{\tilde{G}^{i}(s)\\}_{s\in
S,i\in\mathcal{M}}$ as follows:

$\tilde{H}=\max_{\alpha\in\mathcal{P}^{+}(\mathcal{M})}\bigg{(}\sup_{s\in\partial^{\alpha}\hat{\mathcal{F}}_{c}}\big{|}\min_{i\in\alpha}(\tilde{G}^{i}(s))\big{|}\bigg{)}$
(4)

Quantiles of the distribution of $H$ may now be estimated empirically from the
observed sampling distribution of $\tilde{H}$. Discussion of how the above
supremum, taken across continuous space, is evaluated in practice is deferred
to Supplementary Results Section S$2.1$. However, we do briefly note here
that, as the location of the true boundary segment
$\partial^{\alpha}\mathcal{F}_{c}$ is in practice unknown, the boundary
segment $\partial^{\alpha}\mathcal{F}_{c}$ has been replaced by an estimate
$\partial^{\alpha}\hat{\mathcal{F}}_{c}$. It should be noted that, as a result
of this substitution, our proposed method may not be applied when
$\hat{\mathcal{F}}_{c}$ is empty.

Several strong references exist that argue and verify the correctness of using
the wild $t-$bootstrap for estimating quantiles of maxima distributions in the
above manner. Of particular note are the works of Chang et al. (2017) and
Bowring et al. (2019), in which the wild $t-$bootstrap is proposed for the
estimation of Gaussian maxima distributions and the generation of confidence
regions, respectively. Discussion of the wild $t-$bootstrap in an imaging
context may also be found in Telschow and Schwartzman (2021). Other
significant references which serve as precursors to this work include
Chernozhukov et al. (2013) and Sommerfeld et al. (2018), in which the Gaussian
multiplier bootstrap was investigated for estimating Gaussian maxima
distributions and generating confidence regions, respectively. More general
introductions to such bootstrapping procedures may be found in Wu (1986),
Efron and Tibshirani (1994) and Hesterberg (2014).

In the form that it has been presented here, the wild $t-$bootstrap requires
that $\\{G^{i}\\}_{i\in\mathcal{M}}$ be symmetric. Further, as noted above,
the wild $t-$bootstrap has been extensively verified for estimating maxima
distributions only in settings in which $\\{G^{i}\\}_{i\in\mathcal{M}}$ are
Gaussian. Here we stress that, whilst the theory presented in Section 2 makes
no assumptions on the distribution of $\\{G^{i}\\}_{i\in\mathcal{M}}$, the
bootstrap theory presented throughout Section 3 requires
$\\{G^{i}\\}_{i\in\mathcal{M}}$ to be Gaussian. In order to utilise the
results of Section 2 when the random variables $\\{G^{i}\\}_{i\in\mathcal{M}}$
are not Gaussian, alternative methods must be employed for estimating the
quantiles of $H$.

##  4 Simulations

###  4.1 Simulation Settings

To assess the correctness and performance of the method, a range of
simulations were conducted using synthetic data. Each simulation was designed
to investigate how the empirical coverage (i.e. the proportion of simulation
instances in which the inclusion
$\hat{\mathcal{F}}^{+}_{c}\subseteq\mathcal{F}_{c}\subseteq\hat{\mathcal{F}}^{-}_{c}$
held) was affected by various secondary factors of practical interest. In this
section, we describe three simulations designed to investigate how the
methods’ performance was influenced by; $(1)$ the degree of overlap between
excursion sets, $(2)$ the presence of correlation between the noise fields for
different study conditions, and $(3)$ the magnitude of the spatial rate of
change of the signal.

In all simulations, for $i\in\mathcal{M}$, the model used to generate the
synthetic data took the form;

$Y^{i}(s)=\mu^{i}(s)+\epsilon^{i}(s),\quad\epsilon^{i}(s)\sim
N(0,\sigma^{i}(s)I_{n}).$ (5)

where $n$ was allowed to vary from $40$ to $500$ (in increments of $20$). The
aim, across all simulations, was to assess the performance of the method when
employed for conjunction inference on the true signal
$\\{\mu^{i}\\}_{i\in\mathcal{M}}$ (i.e. when asked the question “Where do all
of the $\\{\mu^{i}\\}_{i\in\mathcal{M}}$ exceed a predefined threshold?”).

The mechanism used to generate the true signal,
$\\{\mu^{i}\\}_{i\in\mathcal{M}}$, varied across simulations. In Simulations 1
and 2, $\\{\mu^{i}\\}_{i\in\mathcal{M}}$ were generated using two binary
images of circles and squares, respectively, positioned in a ‘Venn diagram’
arrangement. Each binary image was scaled by a predefined amount and then
smoothed using an isotropic Gaussian filter with a Full-Width Half Maximum
(FWHM) of $5$ pixels. In Simulation $3$, $\mu^{1}$ and $\mu^{2}$ were
simulated as a horizontal and vertical linear ramp, respectively, with
predefined gradients.

Each simulation was performed twice, once using noisier ‘low-Signal-to-Noise
Ratio (SNR)’ synthetic data and again using less noisy ‘high-SNR’ synthetic
data. To generate the high-SNR synthetic data for Simulations $1$ and $2$, all
simulated $\\{\mu^{i}\\}_{i\in\mathcal{M}}$ were scaled by a factor of $3$
prior to smoothing. This data generation process is identical to that employed
in Bowring et al. (2019), in which it is noted that synthetic data generated
using these parameters strongly resembles that of an fMRI analysis when voxels
are of dimension $2$mm3. In Simulation $3$, for the high-SNR data, the linear
ramps were initially simulated with a gradient of $8$ per $50$ pixels. Across
all simulations, the low-SNR data was generated by reducing the signal
magnitude of the high-SNR data by a factor of $4$. In all simulations,
thresholds of $c=1/2$ and $c=2$ were employed for the low-SNR data and high-
SNR data, respectively.

As Simulation $1$ aimed to investigate the method’s performance as the overlap
between excursion sets increased, in this simulation, the distance between the
centre of the circles was varied, ranging from $0$ pixels (full overlap) to
$50$ pixels (barely overlapping) in increments of $2$ pixels. Similarly, as
Simulation $2$ aimed to assess the effect of correlation between
$\epsilon^{1}$ and $\epsilon^{2}$, in this simulation, this correlation was
varied between $-1$ and $1$ in increments of $0.1$ (in all other simulations,
the noise fields were generated independently). As Simulation $3$ served to
assess how sensitive the empirical coverage was to the rate of change of
$\mu^{1}$ and $\mu^{2}$ over space, in this simulation, the initial ramp
gradients were multiplied by a factor of $k$, where $k$ was varied from $0.25$
to $1.75$ in $0.05$ increments.

All simulation results were obtained as averages taken across $2500$
simulation instances and compared to the nominal coverage using $95\%$
binomial confidence intervals. In all simulation instances, the image
dimensions were $(100\times 100)$ pixels, confidence regions were computed for
a tolerance level of $\alpha=0.05$ using $5000$ bootstrap realizations,
$\sigma^{i}(s)$ was computed using the ordinary least squares estimator and
$\tau_{n}$ was computed as $\tau_{n}=n^{-0.5}$. In each simulation instance,
the assessment of whether the inclusion statement,
$\hat{\mathcal{F}}^{+}_{c}\subseteq\mathcal{F}_{c}\subseteq\hat{\mathcal{F}}^{-}_{c}$,
held was verified using the interpolation-based methods of Bowring et al.
(2019). In this approach, the inclusion statement was deemed to hold if, and
only if, the sets of pixels which were identified as
$\hat{\mathcal{F}}_{c}^{-}$, $\mathcal{F}_{c}$ and $\hat{\mathcal{F}}_{c}^{+}$
were appropriately nested and, in addition, interpolation along the boundary
$\partial\mathcal{F}_{c}$ confirmed that no violations of the inclusion
statement had occurred.

A further six simulations which investigated the influence of other secondary
factors of interest were also conducted. The secondary factors considered by
these simulations include; the presence of spatial structure in the noise, the
effect of $\partial^{\\{1\\}}\mathcal{F}_{c}$ and
$\partial^{\\{2\\}}\mathcal{F}_{c}$ having differing lengths, the effect of
$\epsilon^{1}$ having a much larger variance than $\epsilon^{2}$ and the
impact of varying the value of $M$. A full discussion of these simulations is
deferred to Supplementary Results Section S$4$. In addition, in keeping with
the previous literature, tolerance levels of $\alpha=0.1$ and $\alpha=0.2$, as
well as the substitution of $\partial\mathcal{F}_{c}$ for
$\partial\hat{\mathcal{F}}_{c}$ in Equation (4), were also considered for
simulation (c.f. Sommerfeld et al. (2018), Bowring et al. (2019), Bowring et
al. (2021)). Again, the results of these simulations may be found in the
Supplementary Results document in Sections S$5$ and S$6$.

###  4.2 Simulation Results

For a majority of the simulations conducted, the empirical coverage estimates
were tightly distributed around the expected nominal coverage. In particular,
for the high-SNR data, across all simulations, the $95\%$ binomial confidence
intervals consistently captured the nominal coverage at the predicted rate.
However, when the simulations employed the low-SNR data, it can be seen that
the confidence regions produced by the method were conservative when the data
were generated under certain unfavourable conditions.

Figure 3:  Simulation results generated using the low-SNR synthetic data.
Nominal coverage is shown as a dashed line, with a corresponding binomial
confidence interval also shaded. All data points are averages taken across
$2500$ simulation instances.

For example, in Fig. 3, it can be seen that some over-coverage was observed
for Simulation $1$ when the data was noisy, and the number of subjects was
low. This observation is corroborated by the results of the remaining low-SNR
synthetic data simulations (c.f. Supplementary Results Section S$5$). However,
in all cases, the degree of agreement between the empirical and nominal
coverage improved as the sample size increased. Furthermore, no such over-
coverage was observed for the equivalent high-SNR simulations. This
observation matches expectation, as Theorem 2.1 holds only asymptotically, and
is supported by the previous literature on confidence regions, in which
similar results were observed for the single-‘study condition’ setting (c.f.
Sommerfeld et al. (2018), Bowring et al. (2019)).

The presence of strong negative correlation between the study conditions and
the rate at which the signal varied across space both also appeared to
influence the methods’ performance (c.f. Fig 3, Simulations 2 and 3). In both
instances, it is likely that the observed over-coverage is due to difficulties
in estimating the true location of the boundary $\partial\mathcal{F}_{c}$, as
each of these factors make it harder to identify the locations at which small
changes in $\min_{i\in\mathcal{M}}\mu^{i}$ occur. When high-SNR data was used
in the place of low-SNR data, both of these factors had a substantially
reduced impact on the empirical coverage.

Despite the above observations, the simulation results overwhelmingly
supported the claim that the proposed method is robust to a range of secondary
factors. Included in the list of such factors are; variation in the shape and
size of $\mathcal{F}_{c}$, spatial correlation in the noise, variation in the
lengths of individual boundary segments, between-‘study condition’
correlation, realistic levels of variation in the magnitude of the noise,
large numbers of study conditions and situations in which boundary segments
are shared by many excursion sets. For further detail, see Supplementary
Results Sections S$5$ and S$6$.

In terms of time efficiency, the wild $t-$bootstrap provided extremely fast
computational performance. For instance, using an Intel(R) Xeon(R) Gold 6126
2.60GHz processor with 16GB RAM, and averaging over the $2500$ simulation
instances conducted for $n=500$ high-SNR observations, in Simulation $1$, the
time taken to generate confidence regions for a circle separation of $20$
pixels was $5.98$ seconds. For a comprehensive summary of computation times,
see Supplementary Results Section S$7$.

##  5 Real Data Application

To demonstrate the method in practice, we apply it to fMRI data drawn from the
Human Connectome Project (HCP) dataset (Van Essen et al. (2013)). In this
example, task fMRI data were collected as part of a block design from $80$
healthy, unrelated, young adults as part of a working memory task that had
four distinct components, each using a different stimuli type: pictures of
places, tools, faces and body parts. In Supplementary Results Section S$3$, we
provide a brief overview of the imaging acquisition protocol, task paradigm,
preprocessing stages and first-level analysis employed for generating this
dataset (for exhaustive detail, see Barch et al. (2013) and Glasser et al.
(2013)). Following first-level analysis, the data consisted of $4$ images for
each of the $80$ subjects, measuring the %BOLD response to each of the $4$
stimuli types. In this example, our interest lies in identifying, at the
group-level, which regions of the brain are associated with working memory
regardless of stimuli type (i.e. “Which regions are active in response to all
four working memory stimuli types?”).

As this work has predominantly focused on two-dimensional excursion sets, a
single slice of the brain was chosen for analysis ($z=46$mm, covering portions
of the frontal gyrus involved in working memory). For each stimuli type, an
$n=80$ group-level linear model was constructed. Using the proposed method,
confidence regions were then generated at the $5\%$ confidence level to assess
where the group-level percentage BOLD change exceeded $1\%$ for all four
stimuli types. To compare the proposed method to standard fMRI inference
procedures, a group-level contrast was also generated using the Big Linear
Model toolbox for each of the four stimuli types. In addition, single-‘study
condition’ confidence regions were also generated using the methods of
Sommerfeld et al. (2018) for each of the four stimuli types.

The results are shown in Fig. 4. The conjunction inference identified several
localized regions, including the Superior Frontal Gyrus, which is well known
for its involvement in working memory (c.f. Boisgueheneuc et al. (2006), Vogel
et al. (2016), Alagapan et al. (2019)), and the Angular Gyri, which are well
known to be involved in a range of tasks including memory retrieval, attention
and spatial cognition (c.f. Seghier (2013), Bréchet et al. (2018)).

Figure 4:  $\%$BOLD images (top) and 95% confidence regions (bottom left) for
the HCP working memory task for each of the four visual stimuli types,
alongside conjunction confidence regions assessing the overlap of the four
thresholded $\%$BOLD images (bottom right). Displayed is axial slice $z=46$mm.
The thresholds employed were $c=1\%$ BOLD change, for all images, and
$\alpha=0.05$, for the confidence regions. Upper and lower confidence regions
are displayed in red and blue, respectively. Across the bottom row, the yellow
sets are the point estimates
$\hat{\mathcal{A}}_{c}^{1},...\hat{\mathcal{A}}_{c}^{4}$ and
$\hat{\mathcal{F}}_{c}$ respectively. The red conjunction set has localized
regions in the Superior Frontal Gyrus (top) and the Angular Gyri (left/right),
for which we can assert with $95\%$ confidence that, for all four study
conditions, there was (at least) a $1.0\%$ change in BOLD response.

Whilst the confidence regions for conjunction inference illustrated in Fig. 4
exhibit a clear resemblance to the four contrast images, we note that the red
set, $\hat{\mathcal{F}}_{c}^{+}$, is very small in comparison to the blue set,
$\hat{\mathcal{F}}_{c}^{-}$, which is large and diffuse. This observation
conveys important information about the spatial variability of the regions
identified. In particular, it can be seen that, in the regions immediately
surrounding the Superior Frontal Gyrus and Angular Gyri, there is a much
higher resemblance between the blue ($\hat{\mathcal{F}}_{c}^{-}$) and yellow
($\hat{\mathcal{F}}_{c}$) sets than is seen across the rest of the brain. This
resemblance provides an insight into the degree of spatial variation present
surrounding these regions. In this case, the strength and localization of this
resemblance suggest that the spatial variability of $\hat{\mathcal{F}}_{c}$ is
less severe in and around the aforementioned anatomical regions than it is
across the rest of the brain. It may be concluded that the estimated yellow
‘blobs’ corresponding to the Superior Frontal Gyrus and the Angular Gyri have
been more reliably localized than the other yellow ‘blobs’ appearing in the
image.

In general, the single-‘study condition’ confidence regions can be seen to be
much larger than those observed for the conjunction inference. This
observation matches expectation as the overlap of the four excursion sets is
smaller than each set individually. In addition, in this example, the
conjunction method employed the entire experimental design for analysis (as
opposed to the single-‘study condition’ method which employed only the data
concerning the stimuli of interest) and is therefore expected to have a higher
statistical power and, thus, tighter confidence bounds.

##  6 Conclusion and Discussion

In this work, we have produced a method for generating confidence regions for
intersections and unions of multiple excursion sets. The confidence regions
generated by the method generalise the notion of confidence intervals to
arbitrary spatial dimensions and possess the usual frequentist interpretation
that, were the $0.95$ procedure repeated on numerous samples, the proportion
of confidence regions, $\hat{\mathcal{F}}^{+}_{c}$ and
$\hat{\mathcal{F}}^{-}_{c}$, that correctly enclose $\mathcal{F}_{c}$ would
tend to $0.95$. Such confidence regions serve as an indicator of the
reliability of $\hat{\mathcal{F}}_{c}$ as an estimate of $\mathcal{F}_{c}$.
Confidence regions which closely resemble the set $\hat{\mathcal{F}}_{c}$ may
be interpreted as signifying that $\hat{\mathcal{F}}_{c}$ is a reliable
estimate of $\mathcal{F}_{c}$, whilst little resemblance may suggest that
there is a high degree of spatial variability present in the data and that the
estimated shape, size and locale of $\hat{\mathcal{F}}_{c}$ are not
particularly reliable. Such statements are of particular value in imaging
applications where the aim of a statistical analysis is typically to assess
how reliably some form of activity may be localized to a particular spatial
region.

We stress here that, whilst $\hat{\mathcal{F}}_{c}$ is equal to the
intersection of $\\{\hat{\mathcal{A}}_{c}^{i}\\}_{i\in\mathcal{M}}$, the
confidence regions $\hat{\mathcal{F}}^{\pm}_{c}$ are not the intersections of
the single-‘study condition’ confidence regions
$\\{\hat{\mathcal{A}}^{\pm,i}_{c}\\}_{i\in\mathcal{M}}$ which would be
obtained using the methods of Sommerfeld et al. (2018). This can be seen by
noting that $\\{\hat{\mathcal{A}}^{\pm,i}_{c}\\}_{i\in\mathcal{M}}$ are
defined using separate thresholds from different bootstraps and can be heavily
influenced by the behavior of
$\\{\partial\mathcal{A}_{c}^{i}\\}_{i\in\mathcal{M}}$ in spatial regions far
from $\mathcal{F}_{c}$. In fact, treating the naive intersection of the
($1-\alpha$) confidence regions
$\\{\hat{\mathcal{A}}^{\pm,i}_{c}\\}_{i\in\mathcal{M}}$ as confidence regions
for $\mathcal{F}_{c}$ is not a valid method in general and can result in
undesirable asymptotic coverage lying anywhere inside the range
$[1-M\alpha,1]$. This claim is further discussed in Supplementary Theory
Section S4. To our knowledge, the only valid approach for generating
confidence regions for conjunction inference is that outlined in this work.

One potential limitation of the specific methods of Section 3 is that they
allow for confidence regions to be generated only when the target functions,
$\\{\mu^{i}\\}_{i\in\mathcal{M}}$, are linear combinations of regression
parameter estimates. As noted by Bowring et al. (2021), it is often preferable
for an analysis to consider standardized effect sizes, such as Cohen’s $d$
(Cohen (2013)) or Hedges’ $G$ (Hedges (1981)), instead of regression parameter
estimates, as standardized effect sizes allow for information about
statistical power to be incorporated into the analysis. To this end, Bowring
et al. (2021) outlined an approach for generating confidence regions (in the
single-‘study condition’ setting) for images of Cohen’s $d$ effect sizes.
Here, we suggest that a potential avenue for future work is to combine the
methods of Bowring et al. (2021) and those proposed in this work to allow for
conjunction, or disjunction, inference to be performed on standardized effect
sizes rather than regression parameter estimates.

Another limitation noted here is that, whilst our proposed method
theoretically works for data of arbitrary dimensions, the simulations of
Section 4 verify the performance of our method only for two-dimensional data.
Whilst it is typical for many practical applications to focus on two-
dimensional data (e.g. climate maps, surface images), higher-dimensional
datasets are abundant in a wealth of disciplines (e.g. brain images,
astrological maps). For this reason, we intend to investigate further the
performance of the method in higher dimensions in future work.

## Acknowledgment(s)

This work was partially supported by the NIH under Grant [R01EB026859] (A.S.,
T.M., T.N.).

Data were provided by the Human Connectome Project, WU-Minn Consortium
(Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657)
funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for
Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at
Washington University. The code used to obtain the results of Sections 4 and 5
is freely available and may be found at:
https://github.com/TomMaullin/ConfSets

The authors would like to thank Dr. Fabian J.E. Telschow for his early
involvement in formulating the project, and Professor Gary R. Lawlor for his
email correspondence regarding the differentiability requirements for
L’hospital’s rule for multivariable functions.

## Disclosure statement

The authors report there are no competing interests to declare.

## References

  * Adler [1981] R.J. Adler.  _The geometry of random fields_.  Wiley series in probability and mathematical statistics. Probability and mathematical statistics. J. Wiley, 1981.  ISBN 9780471278443.
  * Adler and Taylor [2009] R.J. Adler and J.E. Taylor.  _Random Fields and Geometry_.  Springer Monographs in Mathematics. Springer New York, 2009.  ISBN 9780387481166.
  * Adler [1977] Robert J. Adler.  Hausdorff Dimension and Gaussian Fields.  _The Annals of Probability_ , 5(1):145 – 151, 1977.  doi: 10.1214/aop/1176995900.
  * Adler et al. [2017] Robert J. Adler, Kevin Bartz, Sam C. Kou, and Anthea Monod.  Estimating thresholding levels for random fields via euler characteristics, 2017.
  * Alagapan et al. [2019] Sankaraleengam Alagapan, Eldad Lustenberger, Caroline; Hadar, Hae Won Shin, and Flavio Fröhlich.  Low-frequency direct cortical stimulation of left superior frontal gyrus enhances working memory performance.  _NeuroImage_ , 184, 01 2019.  doi: 10.1016/j.neuroimage.2018.09.064.
  * Azas and Wschebor [2009] Jean-Marc Azas and Mario Wschebor.  Level sets and extrema of random processes and fields.  _Level Sets and Extrema of Random Processes and Fields_ , 01 2009\.  doi: 10.1002/9780470434642.
  * Barch et al. [2013] Deanna M. Barch, Gregory C. Burgess, Michael P. Harms, Steven E. Petersen, Bradley L. Schlaggar, Maurizio Corbetta, Matthew F. Glasser, Sandra Curtiss, Sachin Dixit, Cindy Feldt, Dan Nolan, Edward Bryant, Tucker Hartley, Owen Footer, James M. Bjork, Russ Poldrack, Steve Smith, Heidi Johansen-Berg, Abraham Z. Snyder, and David C. Van Essen.  Function in the human connectome: Task-fmri and individual differences in behavior.  _NeuroImage_ , 80:169–189, 2013.  ISSN 1053-8119.  doi: https://doi.org/10.1016/j.neuroimage.2013.05.033.  Mapping the Connectome.
  * Berger [1997] Roger L. Berger.  _Likelihood Ratio Tests and Intersection-Union Tests_ , pages 225–237.  Birkhäuser Boston, Boston, MA, 1997.  ISBN 978-1-4612-2308-5.  doi: 10.1007/978-1-4612-2308-5˙15.
  * Berger and Hsu [1996] Roger L. Berger and Jason C. Hsu.  Bioequivalence trials, intersection-union tests and equivalence confidence sets.  _Statistical Science_ , 11(4):283 – 319, 1996\.  doi: 10.1214/ss/1032280304.
  * Boisgueheneuc et al. [2006] Foucaud du Boisgueheneuc, Richard Levy, Emmanuelle Volle, Magali Seassau, Hughes Duffau, Serge Kinkingnehun, Yves Samson, Sandy Zhang, and Bruno Dubois.  Functions of the left superior frontal gyrus in humans: a lesion study.  _Brain_ , 129(12):3315–3328, 09 2006.  ISSN 0006-8950.  doi: 10.1093/brain/awl244.
  * Bowring et al. [2019] Alexander Bowring, Fabian Telschow, Armin Schwartzman, and Thomas E. Nichols.  Spatial confidence sets for raw effect size images.  _NeuroImage_ , 203:116187, 2019.  ISSN 1053-8119.  doi: https://doi.org/10.1016/j.neuroimage.2019.116187.
  * Bowring et al. [2021] Alexander Bowring, Fabian J.E. Telschow, Armin Schwartzman, and Thomas E. Nichols.  Confidence sets for cohen’s d effect size images.  _NeuroImage_ , 226:117477, 2021.  ISSN 1053-8119.  doi: https://doi.org/10.1016/j.neuroimage.2020.117477.
  * Bréchet et al. [2018] Lucie Bréchet, Petr Grivaz, Baptiste Gauthier, and Olaf Blanke.  Common recruitment of angular gyrus in episodic autobiographical memory and bodily self-consciousness.  _Frontiers in Behavioral Neuroscience_ , 12, 11 2018.  doi: 10.3389/fnbeh.2018.00270.
  * Cao [1999] J. Cao.  The size of the connected components of excursion sets of $\chi$2, t and f fields.  _Advances in Applied Probability_ , 31(3):579–595, 1999.  doi: 10.1239/aap/1029955192.
  * Chang et al. [2017] Chung Chang, Xuejing Lin, and R. Todd Ogden.  Simultaneous confidence bands for functional regression models.  _Journal of Statistical Planning and Inference_ , 188:67–81, 2017.  ISSN 0378-3758.  doi: https://doi.org/10.1016/j.jspi.2017.03.002.
  * Chernozhukov et al. [2013] Victor Chernozhukov, Denis Chetverikov, and Kengo Kato.  Gaussian approximations and multiplier bootstrap for maxima of sums of high-dimensional random vectors.  _The Annals of Statistics_ , 41(6), Dec 2013.  ISSN 0090-5364.  doi: 10.1214/13-aos1161.
  * Cohen [2013] J. Cohen.  _Statistical Power Analysis for the Behavioral Sciences_.  Elsevier Science, 2013.  ISBN 9781483276489.  URL https://books.google.co.uk/books?id=rEe0BQAAQBAJ.
  * Dijkstra et al. [2021] Nadine Dijkstra, Simon van Gaal, Linda Geerligs, Sander E Bosch, and Marcel van Gerven.  No overlap between unconscious and imagined representations, May 2021\.
  * Efron and Tibshirani [1994] Bradley Efron and R.J. Tibshirani.  _An Introduction to the Bootstrap_.  Chapman and Hall/CRC, 1 edition, 1994.
  * Ferreira et al. [2021] Roberto A. Ferreira, David Vinson, Ton Dijkstra, and Gabriella Vigliocco.  Word learning in two languages: Neural overlap and representational differences.  _Neuropsychologia_ , 150:107703, 2021.  ISSN 0028-3932.  doi: https://doi.org/10.1016/j.neuropsychologia.2020.107703.
  * Glasser et al. [2013] Matthew F. Glasser, Stamatios N. Sotiropoulos, J. Anthony Wilson, Timothy S. Coalson, Bruce Fischl, Jesper L. Andersson, Junqian Xu, Saad Jbabdi, Matthew Webster, Jonathan R. Polimeni, David C. Van Essen, and Mark Jenkinson.  The minimal preprocessing pipelines for the human connectome project.  _NeuroImage_ , 80:105–124, 2013.  ISSN 1053-8119.  doi: https://doi.org/10.1016/j.neuroimage.2013.04.127.  Mapping the Connectome.
  * Hedges [1981] Larry V. Hedges.  Distribution theory for glass’s estimator of effect size and related estimators.  _Journal of Educational Statistics_ , 6(2):107–128, 1981.  ISSN 03629791.  URL http://www.jstor.org/stable/1164588.
  * Hesterberg [2014] Tim Hesterberg.  What teachers should know about the bootstrap: Resampling in the undergraduate statistics curriculum, 2014.
  * Pranav et al. [2019] Pratyush Pranav, Rien van de Weygaert, Gert Vegter, Bernard J T Jones, Robert J Adler, Job Feldbrugge, Changbom Park, Thomas Buchert, and Michael Kerber.  Topology and geometry of gaussian random fields i: on betti numbers, euler characteristic, and minkowski functionals.  _Monthly Notices of the Royal Astronomical Society_ , 485(3):4167–4208, Feb 2019.  ISSN 1365-2966.  doi: 10.1093/mnras/stz541.
  * Seghier [2013] Mohamed L. Seghier.  The angular gyrus: multiple functions and multiple subdivisions.  _The Neuroscientist : a review journal bringing neurobiology, neurology and psychiatry_ , 19(1):43–61, Feb 2013.  ISSN 1089-4098.  doi: 10.1177/1073858412440596.
  * Sommerfeld et al. [2018] Max Sommerfeld, Stephan Sain, and Armin Schwartzman.  Confidence regions for spatial excursion sets from repeated random field observations, with an application to climate.  _Journal of the American Statistical Association_ , 113(523):1327–1340, 2018.  doi: 10.1080/01621459.2017.1341838.  PMID: 31452557.
  * Telschow and Schwartzman [2021] Fabian Telschow and Armin Schwartzman.  Simultaneous confidence bands for functional data using the gaussian kinematic formula.  _Journal of Statistical Planning and Inference_ , 216, 06 2021.  doi: 10.1016/j.jspi.2021.05.008.
  * Torres [1994] Sergio Torres.  Topological Analysis of COBE-DMR Cosmic Microwave Background Maps.  _apjl_ , 423:L9, March 1994.  doi: 10.1086/187223.
  * Van Essen et al. [2013] David Van Essen, Stephen Smith, Deanna Barch, Timothy Behrens, Essa Yacoub, and Kamil Ugurbil.  The wu-minn human connectome project: an overview.  _NeuroImage_ , 80, 05 2013.  doi: 10.1016/j.neuroimage.2013.05.041.
  * Vogel et al. [2016] Tobias Vogel, Renata Smieskova, Anna Schmidt, André; Walter, Fabienne Harrisberger, Anne Eckert, Undine E. Lang, Anita Riecher-Rössler, Marc Graf, and Stefan Borgwardt.  Increased superior frontal gyrus activation during working memory processing in psychosis: Significant relation to cumulative antipsychotic medication and to negative symptoms.  _Schizophrenia Research_ , 4 2016.  doi: 10.1016/j.schres.2016.03.033.
  * Worsley [1996] Keith J. Worsley.  The geometry of random images.  _CHANCE_ , 9(1):27–40, 1996.  doi: 10.1080/09332480.1996.10542483.
  * Wu [1986] C. F. J. Wu.  Jackknife, Bootstrap and Other Resampling Methods in Regression Analysis.  _The Annals of Statistics_ , 14(4):1261 – 1295, 1986.  doi: 10.1214/aos/1176350142.
  * Zhang et al. [2021] Xiaoxia Zhang, Linling Li, Gan Huang, Li Zhang, Zhen Liang, Li Shi, and Zhiguo Zhang.  A multisensory fmri investigation of nociceptive-preferential cortical regions and responses.  _Frontiers in neuroscience_ , Apr 2021.


#  A slime mold inspired local adaptive mechanism for flow networks

Vidyesh Rao Anisetti,1, Ananth Kandala,2 J. M. Schwarz1,3 1 Physics
Department, Syracuse University, Syracuse, NY 13244 USA
2Department of Physics, University of Florida, FL 32611-8440, USA
3 Indian Creek Farm, Ithaca, NY 14850 USA

###### Abstract

In the realm of biological flow networks, the ability to dynamically adjust to
varying demands is paramount. Drawing inspiration from the remarkable
adaptability of Physarum polycephalum, we present a novel physical mechanism
tailored to optimize flow networks. Central to our approach is the principle
that each network component—specifically, the tubes— harnesses locally
available information to collectively minimize a global cost function. Our
findings underscore the scalability of this mechanism, making it feasible for
larger, more complex networks. We construct a comprehensive phase diagram,
pinpointing the specific network parameters under which successful adaptation,
or tuning, is realized. There exists a phase boundary in the phase diagram,
revealing a distinct satisfiability-unsatisfiability (SAT-UNSAT) phase
transition delineating successful and unsuccessful adaptation.

##  I Introduction

Biological circulatory systems, despite their intricacy, exhibit remarkable
adaptability. They adeptly modify various attributes—such as vessel diameter,
wall thickness, and the count of micro-vessels cater to the evolving metabolic
needs of tissues hu2012blood ; Humphrey2021 . This adaptability can be
perceived as an outcome of an optimization mechanism, where a global cost
function is optimized. Intriguingly, this optimization is not typically
directed by a central authority, but emerges from local interactions between
different aspects of the system. For example, prior research illustrates how
vascular flow efficiency ronellenfitsch2016global  can be enhanced by altering
attributes, e.g. tube thickness, based on local data such as flow within a
tube.

Given the circulatory system example, as well as others Tero2010 ; Eytan2003 ,
a fundamental inquiry centers on how local interaction rules give rise to
adaptive behaviors, and how these behaviors subsequently manifest as
optimization algorithms in biological systems Holland1992 . In this
manuscript, we introduce a straightforward physical mechanism that potentially
allows biological systems to implement such optimization strategies. Moreover,
limitations on such optimization algorithms do indeed exist. For instance, one
type of optimization algorithm for a large swath of some parameter space may
not be feasible. Here, we explore the limitations on our physical mechanism
for adaptation in the context of what is known as a Satisfiability/SAT-
Unsatisfiability/UNSAT phase transition Mitchell1992 ; Kirkpatrick1994 ;
Monasson1999 . A similar transition was found in the study of limits on the
multifunctionality of tunable flow and mechanical networks Rocks2019 .

The specific question we are interested in is: How does one tune the node
pressures of a flow network, by altering its tube thickness, via a physical
mechanism that uses only local information? To answer this question, we take
inspiration from the adaptive behaviour observed in Physarum polycephalum, or
slime mold. This organism, in spite of not having a brain or nervous system,
is capable of optimizing its network structure to determine the most efficient
path between food sources Tero2010 ; Tero2008 . Prior research Alim2017
indicates that slime mold possibly utilizes a form of chemical signaling for
this purpose. Upon encountering food, it releases a signaling chemical at the
location, which is then propagated throughout the organism via advection. This
chemical, in turn, induces an expansion in the tubes through which it is
flowing. As a result, the tubes that connect food sources through the shortest
path undergo more pronounced expansion compared to those on longer routes,
leading to an efficient connection between food sources. This behavior
exemplifies how biological systems can harness physical processes to optimize
for a specific task such as finding food.

Figure 1: Schematic for an adaptative flow mechanism. [Arrow 1] To increase
the pressure at node A, a specialized chemical (depicted in red) is
introduced. This chemical is advectively transported by the fluid flow within
the network. The fluid flow direction is indicated by the grey arrows. [Arrow
2] As this chemical traverses the tubes, it interacts with the tube’s
structure, causing it to constrict and thereby increasing the flow resistance.
This results in an increase in pressure at node A. [Arrows 3 & 4] Conversely,
to decrease the pressure at node A, a different chemical is released that
dilates the tubes. This dilation reduces flow resistance, leading to a
decrease in pressure at node A.

Our system is a flow network - a set of nodes connected by tubes. Fluid flows
in and out of the network through some of those nodes. This flow creates
pressure values at each node, depending upon the architecture of the network
and the conductance of the pipes. The task here is to modify the conductance
of these pipes such that the pressure values at some ‘output’ nodes matches
the desired pressures.

The first thing that we observe is that pressure at a node depends on the
resistance to flow downstream. Therefore, to alter the pressure at a node, we
must change the conductance of pipes downstream (Fig. 1). Consider an output
node where we want the pressure to decrease. To do so a chemical at that node
gets released which gets carried by the fluid flow. This chemical interacts
with the tube such that, when it is flowing through a tube, it increases the
conductance of the pipe by making it thicker. This increase in conductance
decreases the resistance to flow, which in turn decreases the pressure at the
output node. Similarly, when we wish to increase the output node pressure, we
must release a different kind of chemical which decreases the conductance of
the pipes by making it thinner. Through this mechanism, the entire network
undergoes adjustments. Each tube, relying on the locally available
information—namely, the chemical concentration within it—fine-tunes the
pressures at the output nodes to align with the target values. This localized
adjustment is facilitated by our introduced method, where the discrepancy at
the output nodes is conveyed into the network through a chemical signal,
subsequently influencing the network structure.

In what follows, we have assessed the performance of our tuning mechanism
across a range of network sizes. The scaling relations observed suggest that
our tuning approach remains effective even as the network size increases.
Notably, our results suggest a SAT-UNSAT phase transition
monasson2007introduction  with a discontinuous jump in the fraction of
networks that can be successfully adapted at the transition point along with
universal scaling features.

##  II The Tuning Process

###  II.1 The System

We create networks consisting of N nodes and M edges. We first create a
Barabási-Albert network with connection parameter 1 using the python package,
NetworkX. This creates a minimally connected network with N nodes and N-1
edges. To create a network with M edges we add M-(N+1) unique edges.

We select $j$ boundary nodes, denoted as ${q_{1},q_{2},...,q_{j}}$, and apply
an external potential $\mathbf{u}=[u(q_{1}),u(q_{2}),...,u(q_{j})]^{T}$. The
resulting response potentials at the remaining nodes, termed interior nodes,
are calculated by solving the discrete Laplace equation using the Schur
compliment method Curtis . From these interior nodes, we identify $k$ output
nodes ${p_{1},p_{2},..,p_{k}}$. The potentials at these nodes are represented
as ${v(p_{1}),v(p_{2}),...,v(p_{k})}$. The objective is to adjust this network
by altering the conductance of its pipes, aiming to align the output potential
vector $\mathbf{v}=[v(p_{1}),v(p_{2}),...,v(p_{k})]^{T}$ with a target output
potential vector
$\mathbf{v_{des}}=[v_{des}(p_{1}),v_{des}(p_{2}),v_{des}(p_{3}),..,v_{des}(p_{k})]^{T}$.
For context, the target output potential vector $\mathbf{v_{des}}$ is derived
by introducing a perturbation to the potential at each node, with the
perturbation value sourced from a uniform distribution spanning
$[-\Delta,+\Delta]$.

###  II.2 Implementation of the Tuning Process :

  1. 1.

The input potential vector $\mathbf{u}$ is applied at boundary nodes. A
supervisor checks the output potentials $\mathbf{v}$ at output nodes and
compares them to the vector of desired output potentials $\mathbf{{v}_{des}}$
.

  2. 2.

There are two kinds of chemicals, $s_{+}$ and $s_{-}$. $s_{+}$ increases the
conductance of the pipe when it is passing through it, and vice versa for
$s_{-}$. We assume that the output nodes release a chemical whose amount is
proportional to the difference between the present output potentials and
desired output potentials. At $t=0$ for some $a\in O,~{}v(a)\neq{v}_{des}(a)$,
then

$\displaystyle v(a)>{v}_{des}(a)$ $\displaystyle\Rightarrow
s_{+}(a)=\alpha(v(a)-{v}_{des}(a))$ (1) $\displaystyle v(a)<{v}_{des}(a)$
$\displaystyle\Rightarrow s_{-}(a)=\alpha({v}_{des}(a)-v(a))$ (2)

Where $\alpha$ is the factor that controls the chemical response given by the
node to the difference in potentials. Moreover, $s_{+}(a)$ denotes amount of
chemical (e.g. no of molecules) at node ”a”, and $\mathbf{s_{+}}(t)$ denotes
the array of chemical amount at each node at time $t$ .

  3. 3.

This chemical is carried by the current in the network. Therefore, in the next
time step the chemical flows to the neighbouring nodes of $a$ that are
downstream to $a$111For simplicity we assume that the chemical has negligible
diffusion and the only way it can spread is via the network currents. We call
all such downstream neighbours of $a$ as $\mathcal{D}(a)$. Then for all
$b\in\mathcal{D}(a):$

$\displaystyle
s_{+}(b,t+1)=s_{+}(a,t)\times\dfrac{i(b,a)}{\sum_{x\in\mathcal{D}(a)}i(x,a)}$
(3) $\displaystyle+{\rm(incoming~{}chemical~{}from~{}other~{}nodes)}$ (4)

where $i(x,a)$ represents the current from $a$ to $x$. This is how the entire
array $\mathbf{s_{+}}$ and $\mathbf{s_{-}}$ is modified at each time step.
Note that all the chemical initially present at $a$ flows downstream after one
time step.

  4. 4.

Using the above equation, an N $\times$ N array $\hat{S}_{+}$ is generated,
where each entry $i,j$ denotes the amount of chemical passing through the pipe
$\\{i,j\\}$ at step $t$ $\xrightarrow{}$ $t+1$. Let $\hat{W}$ denote the
conductance matrix of the graph, where each entry $\\{i,j\\}$ denotes the
conductance of that pipe. Then

$\displaystyle\hat{W}(t+1)=\hat{W}(t)+\beta(\hat{S}_{+}-\hat{S}_{-}),$ (5)

$\beta$ controls the response of the pipe to the passing chemical.

  5. 5.

The new potentials are calculated on the interior vertices using
$\hat{W}(t+1)$. Again, the supervisor checks if $v(a)=v_{des}(a)$. The
chemical takes some time to reach the boundary nodes, where it drains out of
the network 222The current flows into and out of the network through boundary
nodes. Therefore, the total change in potential due to the chemical released
at the output nodes is observed after some amount of time. Therefore, we
introduce a time delay $\tau$ before releasing the chemical once again 333This
helps the potentials at the output node to converge nicely at the desired
potentials. If this delay is not introduced, the output potential oscillates
about the desired potential. Also, this helps in avoiding the buildup of
excess chemical in the network.

  6. 6.

This process is repeated iteratively.

##  III Results

We implemented the aforementioned procedure on a network comprising 50 nodes
and 118 edges, as depicted in Fig. 2. The conductance values for the flow
network were uniformly sampled from the range $[10^{-5},1]$. External
pressures were applied to the input nodes (highlighted in green), with each
node’s external pressure being uniformly sampled from $[-10,10]$. The
objective was to adjust the resulting pressure at the output nodes (indicated
in red) to align with specified target values. These target pressures were
determined by setting $\Delta=1$. Chemicals were introduced at these nodes at
intervals of $\tau=5$ units, influencing the conductance of downstream pipes.
Successful tuning was evident as the error, represented as
$||\mathbf{v}-\mathbf{v_{des}}||$, decreased significantly, and the pressure
at the output nodes converged to the desired values.

Figure 2: The tuning process. [a,b] The network undergoes modifications due to
the tuning process. The color bar illustrates the conductance values both pre
and post-training. [c] A plot of error versus time. The error plateaus since
tuning a particular node ‘$a$’ is stopped when $|v(a)-v_{des}(a)|<10^{-6}$).
[d] Pressures at the output nodes converge to their target values, represented
by similarly colored translucent horizontal lines.

Fig. 3 presents the $P_{SAT}$ values across varying parameters: number of
edges (E), output nodes (M), total nodes (N), and with a fixed $\Delta=0.1$.
(While simulations were conducted for N = 50, 100, 150, 200, and 250, Fig.3
specifically displays results for N=50, 150, and 250). The other training
parameters remain consistent with those used in Fig. 2. We define $P_{SAT}$ as
the proportion of networks that achieve successful tuning. A tuning process is
deemed successful if the ratio of initial error to final error is less than
$10^{-2}$. This ratio was determined for each pixel through 100 training
repetitions. A notable surge in $P_{SAT}$ is observed, transitioning from a
‘hard phase’—where tuning is ineffective—to an ‘easy phase’ where it
predominantly succeeds. This phenomenon is further elaborated as a SAT-USAT
phase transition in Fig. 4. The relevant tuning parameter in such transitions
is the clause density $\alpha$ eli_sat_unsat ; monasson2007introduction ,
which represents the ratio of clauses to variables. In our context, it is the
ratio of nodes to tune to the number of edges ($M/E$). The right column of
Fig. 3 illustrates the decline of $P_{SAT}$ with an increasing $\alpha$. In
these plots, curves for constant $E$ are shown. Increasing $\alpha$—achieved
by increasing $M$ while maintaining $E$—signifies with an increasing problem
hardness. We fit these curves to a sigmoid-like function given by:

$f(x,a,b)=\frac{1}{2}\left(1-\frac{{e^{b\cdot(x-a)}-e^{-b\cdot(x-a)}}}{{e^{b\cdot(x-a)}+e^{-b\cdot(x-a)}}}\right)$
(6)

From these fits, we deduce the critical point $\alpha_{c}$ and the transition
width $w$. $\alpha_{c}$ is the $\alpha$ value at which $P_{SAT}=0.5$, and $w$
represents the horizontal span between $P_{SAT}=0.25$ and $P_{SAT}=0.75$.

Figure 3: Adaptation phase diagram. (a-c) Depictions of $P_{SAT}$ for networks
comprising 50, 150, and 250 nodes, respectively. The color bar indicates
$P_{SAT}$ values. The x-axis represents the number of output nodes (M), while
the y-axis denotes the number of edges (E). (d-f) Illustrations of the decline
in $P_{SAT}$ in relation to an increase in $\alpha=M/E$. Each curve is
representative of a horizontal slice of the 2D phase diagram (displayed on the
left) at a constant $E$ value. The data is fitted to a sigmoid like curve as
outlined in Eq. 6.

By analyzing how $\alpha_{c}$ and $w$ vary with $E$, we deduce the critical
number of nodes $M_{c}=\alpha_{c}E$ that can be successfully tuned for varying
system sizes, specifically for N=50, 100, 150, 200, and 250. The width around
this critical value, in terms of the number of nodes for varying system sizes,
is given by $\Delta M_{c}=wE$ (refer to Fig. 4). We observe that the critical
number of nodes $\alpha_{c}E$ scales as a power law with respect to $E$ with
an exponent greater than 1 ($M_{c}\sim E^{1.03}$). This indicates that our
tuning process is scalable and will remain effective for larger system sizes.
To demonstrate that this is a SAT-UNSAT phase transition and not just a
crossover, we present two supporting arguments: 1) We observe that $wE$ scales
with $E$ with an exponent less than 1 ($wE\sim E^{0.76}$). This implies that
$w$ scales with $E$ with a negative exponent, indicating that in the
thermodynamic limit as $E\rightarrow\infty$, the transition width $w$
vanishes. 2) We note that this transition is universal, as all the transition
plots in Fig. 3 (d-f) (and for sizes 100 and 150) collapse onto a universal
function upon rescaling by $\alpha\rightarrow(\alpha-\alpha_{c})/w$.

Figure 4: Universal scaling near the transition. a) Illustrates the scaling of
the critical number of nodes that can be tuned with the varying number of
edges $E$ for different system sizes. (b) Depicts the scaling of the width
around this critical number of nodes with respect to the number of edges $E$
for varying system sizes. (c) Shows that the transition plots collapse onto
one another upon rescaling for $N=50,100..250.$

##  IV Discussion

Inspired by Physarum polycephalum, we introduce a physical mechanism to tune
flow networks, where each network component uses locally available information
to optimise a global cost function. We show that this mechanism is scalable to
larger networks and present a phase diagram demonstrating for what values of
network parameters successful tuning is observed. Additionally, we showed that
this optimization strategy exhibits a SAT-UNSAT phase transition.

In previous work, we employed a similar concept to train physical networks
Anisetti2023 , where the gradient information was encoded in a chemical signal
diffused throughout the network. While this approach is equivalent to
executing gradient descent, its biological plausibility is questionable due to
the slow diffusion-based transport of chemical signals. In contrast, our
current methodology involves advective transport of these signals so that the
adaptation occurs over faster time scales. However, with this approach, we do
not anticipate a gradient descent of the MSE, as the chemical signals only
modify the conductances downstream.

Given the parallels with particle jamming, which encompasses a SAT-USAT
transition ruiz-garcia2019tuning , we propose that an analogous tuning
mechanism can be established for hard sphere systems. This has the potential
to lay the groundwork for methods to embed memory keim2019memory  within these
systems and as well learning within these systems given recent work reframing
the learning process as feedback-based aging in a glassy landscape
AnisettiAging . Moreover, the potential applications of our findings extend
further. One immediate and practical application lies in the domain of urban
water supply systems. As cities grow and evolve, their water supply networks
face the challenge of dynamically adjusting to meet varying demands. The
traditional methods of managing these pressures might not be efficient or
responsive enough to cater to the rapid changes in demand. Our proposed
mechanism offers a potential solution to this challenge, enabling water supply
systems to self-optimize and adapt in real-time. Furthermore, the principles
elucidated in our study could inspire innovative methodologies for other
human-made transport systems. For instance, power grids, which require a
delicate balance of supply and demand across vast networks, could benefit from
a similar approach. By integrating mechanisms that allow for localized
adjustments based on real-time feedback, these systems could become more
resilient, efficient, and adaptive.

The authors acknowledge Benjamin Scellier and Karen Alim for helpful
discussions. JMS acknowledges financial support from NSF-DMR-2204312.

## References

  * (1) D. Hu, D. Cai, and A. V. Rangan, “Blood vessel adaptation with fluctuations in capillary flow distribution,” PLOS ONE, 2012.  Published: September 27, 2012.
  * (2) J. D. Humphrey and M. A. Schwartz, “Vascular mechanobiology: homeostasis, adaptation, and disease,” Annual review of biomedical engineering, vol. 23, pp. 1–27, 2021.
  * (3) H. Ronellenfitsch and E. Katifori, “Global optimization, local adaptation, and the role of growth in distribution networks,” Phys. Rev. Lett., vol. 117, p. 138301, 2016.
  * (4) A. Tero, S. Takagi, T. Saigusa, K. Ito, D. P. Bebber, M. D. Fricker, K. Yumiki, R. Kobayashi, and T. Nakagaki, “Rules for biologically inspired adaptive network design,” Science, vol. 327, no. 5964, pp. 439–442, 2010.
  * (5) D. Eytan, N. Brenner, and S. Marom, “Selective adaptation in networks of cortical neurons,” Journal of Neuroscience, vol. 23, no. 28, pp. 9349–9356, 2003.
  * (6) J. H. Holland, Adaptation in natural and artificial systems: an introductory analysis with applications to biology, control, and artificial intelligence.  MIT press, 1992.
  * (7) D. Mitchell, B. Selman, H. Levesque, et al., “Hard and easy distributions of sat problems,” in Aaai, vol. 92, pp. 459–465, 1992.
  * (8) S. Kirkpatrick and B. Selman, “Critical behavior in the satisfiability of random boolean expressions,” Science, vol. 264, no. 5163, pp. 1297–1301, 1994.
  * (9) R. Monasson, R. Zecchina, S. Kirkpatrick, B. Selman, and L. Troyansky, “Determining computational complexity from characteristic ‘phase transitions’,” Nature, vol. 400, no. 6740, pp. 133–137, 1999.
  * (10) J. W. Rocks, H. Ronellenfitsch, A. J. Liu, S. R. Nagel, and E. Katifori, “Limits of multifunctionality in tunable networks,” Proceedings of the National Academy of Sciences, vol. 116, no. 7, pp. 2506–2511, 2019.
  * (11) A. Tero, K. Yumiki, R. Kobayashi, T. Saigusa, and T. Nakagaki, “Flow-network adaptation in Physarum amoebae,” Theory in Biosciences, vol. 127, no. 2, pp. 89–94, 2008.
  * (12) K. Alim, N. Andrew, A. Pringle, and M. P. Brenner, “Mechanism of signal propagation in Physarum polycephalum,” Proceedings of the National Academy of Sciences of the United States of America, vol. 114, no. 20, pp. 5136–5141, 2017.
  * (13) R. Monasson, “Introduction to phase transitions in random optimization problems,” arXiv preprint arXiv:0704.2536, 2007.  Lectures delivered in Les Houches during the Summer School on Complex Systems, July 2006.
  * (14) E. Curtis and J. Morrow, Inverse Problems for Electrical Networks.  Series on applied mathematics, World Scientific, 2000.  pg.40-43, sec. 3.5.
  * (15) For simplicity we assume that the chemical has negligible diffusion and the only way it can spread is via the network currents.
  * (16) The current flows into and out of the network through boundary nodes.
  * (17) This helps the potentials at the output node to converge nicely at the desired potentials. If this delay is not introduced, the output potential oscillates about the desired potential. Also, this helps in avoiding the buildup of excess chemical in the network.
  * (18) E. Chertkov, Phase transitions in random satisfiability problems.  2017\.  Lecture notes available at : https://guava.physics.uiuc.edu/~nigel/courses/563/Essays_2017/PDF/chertkov.pdf.
  * (19) V. R. Anisetti, B. Scellier, and J. M. Schwarz, “Learning by non-interfering feedback chemical signaling in physical networks,” Phys. Rev. Research, vol. 5, p. 023024, 2023.
  * (20) M. Ruiz-García, A. J. Liu, and E. Katifori, “Tuning and jamming reduced to their minima,” Phys. Rev. E, vol. 100, p. 052608, 2019.
  * (21) N. C. Keim, J. D. Paulsen, Z. Zeravcic, S. Sastry, and S. R. Nagel, “Memory formation in matter,” Rev. Mod. Phys., vol. 91, p. 035002, 2019.
  * (22) V. R. Anisetti, A. Kandala, and J. Schwarz, “Emergent learning in physical systems as feedback-based aging in a glassy landscape,” arXiv preprint arXiv:2309.04382, 2023.