# Learning Invariant Stabilizing Controllers for Frequency Regulation under
Variable Inertia ††thanks: This work was supported by NSF Award ECCS-1947050.

Priyank Srivastava Patricia Hidalgo-Gonzalez Jorge Cortés P. Srivastava is
with the Department of Mechanical Engineering, Massachusetts Institute of
Technology<EMAIL_ADDRESS>P. Hidalgo-Gonzalez and J. Cortés are with the
Department of Mechanical and Aerospace Engineering, UC San Diego,
{phidalgogonzalez<EMAIL_ADDRESS>
###### Abstract

Declines in cost and concerns about the environmental impact of traditional
generation have boosted the penetration of renewables and non-conventional
distributed energy resources into the power grid. The intermittent
availability of these resources causes the inertia of the power system to vary
over time. As a result, there is a need to go beyond traditional controllers
designed to regulate frequency under the assumption of invariant dynamics.
This paper presents a learning-based framework for the design of stable
controllers based on imitating datasets obtained from linear-quadratic
regulator (LQR) formulations for different switching sequences of inertia
modes. The proposed controller is linear and invariant, thereby interpretable,
does not require the knowledge of the current operating mode, and is
guaranteed to stabilize the switching power dynamics. We also show that it is
always possible to stabilize the switched system using a communication-free
local controller, whose implementation only requires each node to use its own
state. Simulations on a 12-bus 3-region network illustrate our results.

##  I Introduction

In power networks, any mismatch between electricity generation and consumption
leads to the deviation of the frequency from its nominal value. The increasing
penetration of renewable energy resources (RES), along with their intermittent
availability, has made ensuring frequency regulation more relevant than ever.
The presence of RES reduces the inertia of the system and makes it time-
varying. As such, traditional controllers designed for invariant systems are
no longer guaranteed to be stabilizing. Motivated by these considerations,
this paper addresses the problem of optimally stabilizing the frequency of
power networks with time-varying inertia.

_Literature Review:_ In the traditional paradigm of power systems, there
exists a number of mechanisms to prevent frequency excursions, cf. [1, 2].
Inertial response is the first (automatic) response when any power imbalances
occur. It originates from the kinetic energy stored in synchronous generators
and determines the instantaneous frequency when power imbalances arise. More
inertia in the system translates into a slower rate of change of frequency. As
the frequency starts deviating, some generators respond proportionally to this
deviation through the governor response or droop control [3]. After droop
control starts actuating, slower mechanisms (e.g., spinning reserves)
participate to restore the frequency. RES, such as wind and solar, are usually
connected through inverters, decoupling their rotational inertia (if existing)
from the grid. As a result, the system inertia is an inverse function of the
number of RES. In fact, since different distributed energy resources make
autonomous decisions when connecting to the grid, the inertia of the system
becomes time-varying [4]. This can provoke abrupt variation in the grid
frequency under mismatches of generation and demand. Without appropriate
measures, this can make the standard frequency control schemes too slow to
mitigate arising contingencies. The impact of low inertia in the future grid
is captured by system operators in various reports [5, 6, 7].

A growing body of work addresses this need by analyzing the effect of inertia
variations on frequency [8], designing robust controllers [9], and identifying
conditions on the power supply dynamics and rate of change of inertia that
ensure stability [10]. The work [11] uses a switched affine hybrid system
framework to model inertia variations and proposes a learning-based invariant
controller stabilizing each inertia mode of the closed-loop dynamics. This
formulation is extended in [12] to enhance sparsity, albeit each node needs to
communicate with a certain minimum threshold number of nodes, and take into
account the stability of the switched system, albeit there is no guarantee
that a feasible solution exists. Both [11, 12] assume that all the nodes have
equal inertia in each mode and stability is considered a posteriori once the
training is complete. Ideally, as pointed out in [13], stability guarantees
should be encoded in the training phase itself. In fact, the lack of
guarantees on stability is a shortcoming in many works employing machine
learning techniques for power systems, cf. [14, 15, 16, 17]. For example, the
work in [17] presents an overview on reinforcement learning (RL) techniques
for power systems, but does not touch upon the stability aspects. The work
[18] discusses the importance of stability when using RL in power systems and
how key RL assumptions may not hold in some power systems applications.
Recently, [19, 20] have developed RL approaches for frequency control with
stability guarantees, but the designed controllers do not consider time-
varying frequency dynamics due to the changing inertia.

_Statement of Contributions:_ We consider the problem of designing an
invariant controller to stabilize the frequency of a power network with time-
varying inertia. The fact that the controller is invariant makes it oblivious
to changes in inertia, hence facilitating its implementability by power system
operators. Our starting point is a formulation of the frequency dynamics of
the power network as a switched affine system, where each mode corresponds to
a different value of the inertia. To address the fact that changes in the
operating mode are not known a priori, we consider a candidate set of
switching sequences and, for each of them, solve a finite-horizon LQR problem
to generate optimal state-input trajectories to be used as data. We then
formulate the controller design problem as a constrained least-squares
optimization, where the objective function measures the fit of the
trajectories generated with the controller to the data, and the constraints
encode the stabilization requirement for the switched system. Our first result
considers the formulation where constraints correspond to the stabilization of
the individual modes and its proof is constructive, providing an explicit
stabilizing controller which is distributed over the power network. Our second
result generalizes our treatment to guarantee the system stability under
arbitrary switching, and establishes that regardless of the inertia of the
operating mode, stabilization is always possible using an invariant
controller. Our last result shows that, in fact, there always exists a
stabilizing controller which is local, meaning that its implementation only
requires each node to use its own state. Simulations on a 12-bus 3-region
network with variable inertia demonstrate the stabilizing performance of the
learned controllers with and without sparsity constraint.

##  II Problem Formulation

Consider111Throughout the paper, we use the following notation. Let
$,\mathbb{R}_{\geq 0},\mathbb{R}_{>0},\mathbb{Z},\mathbb{Z}_{>0}$ denote the
set of reals, non-negative reals, positive reals, integers, and positive
integers, respectively. We denote by $|\operatorname{\mathcal{X}}|$ the
cardinality of a set $\operatorname{\mathcal{X}}$. The symbol
$\operatorname{\mathbf{0}}$ represents the matrix of all zeros and $I$ denotes
the identity matrix, with appropriate dimensions. For a matrix $A$, $A_{ij}$
denotes its $ij$th element, $A^{\top}$ denotes its transpose, and $A^{-1}$ its
inverse. $A\succ\operatorname{\mathbf{0}}$
$(\succeq\operatorname{\mathbf{0}})$ and $A\prec\operatorname{\mathbf{0}}$
$(\preceq\operatorname{\mathbf{0}})$ denote, respectively, that $A$ is
positive definite (semidefinite) and negative definite (semidefinite).
$A\otimes B$ denotes the Kronecker product of matrices $A$ and $B$.
$\operatorname{diag}(a_{i})$ is the matrix with entries
$\\{a_{i}\\}_{i=1}^{m}$ in its main diagonal. With a slight abuse of notation,
we let $(x;y)\in^{m+n}$ denote the concatenated vector obtained by putting
together the entries of vectors $x\in^{m}$ and $y\in^{n}$. We also employ
basic concepts from graph theory following [21]. We denote a weighted
undirected graph by
$\operatorname{\mathcal{G}}=(\operatorname{\mathcal{V}},\operatorname{\mathcal{E}},\operatorname{\mathsf{A}})$,
with $\operatorname{\mathcal{V}}$ as the set of nodes and
$\operatorname{\mathcal{E}}\subseteq\operatorname{\mathcal{V}}\times\operatorname{\mathcal{V}}$
as the set of edges. $(i,j)\in\operatorname{\mathcal{E}}$ iff
$(j,i)\in\operatorname{\mathcal{E}}$ iff there is an edge from node $i$ to
$j$. A node $j\in\operatorname{\mathcal{V}}$ is a neighbor of $i$ if
$(i,j)\in\operatorname{\mathcal{E}}$. We denote the set of neighbors of node
$i$ by $\operatorname{\mathcal{N}}_{i}$. With
$|\operatorname{\mathcal{V}}|=n$, the _adjacency matrix_
$\operatorname{\mathsf{A}}\in^{n\times n}$ of $\operatorname{\mathcal{G}}$ is
such that $\operatorname{\mathsf{A}}_{ij}>0$ if
$(i,j)\in\operatorname{\mathcal{E}}$ and $\operatorname{\mathsf{A}}_{ij}=0$,
otherwise. The weighted degree of node $i$ is
$d(i)=\sum_{\operatorname{\mathcal{N}}_{i}}\operatorname{\mathsf{A}}_{ij}$.
Finally, the _Laplacian matrix_ $\operatorname{\mathsf{L}}\in^{n\times n}$ is
$\operatorname{\mathsf{L}}=\operatorname{diag}(d(i))-\operatorname{\mathsf{A}}$.
a power network with $n\in\mathbb{Z}_{>0}$ nodes, whose interconnection is
described by an undirected graph $\operatorname{\mathcal{G}}$. Following [22],
we consider a DC approximation of the power flow. The frequency and phase
angle dynamics for each node $i\in\\{1,\dots,n\\}$ are approximated as follows

$\displaystyle
m_{i}\ddot{\theta}_{i}+d_{i}\dot{\theta}_{i}=u_{i}-\sum\limits_{\operatorname{\mathcal{N}}_{i}}b_{ij}(\theta_{i}-\theta_{j}),$

where $u_{i}$ is the power input at node $i$ and $b_{ij}\in\mathbb{R}_{\geq
0}$ is the susceptance between lines $i$ and $j$. If node $i$ is a synchronous
generator, then $\theta_{i}\in$ denotes the rotor angle,
$m_{i}\in\mathbb{R}_{>0}$ the rotational inertia of the generator $i$ and
$d_{i}\in\mathbb{R}_{>0}$ the primary speed droop control at node $i$. If node
$i$ corresponds to a renewable or battery interfaced via a power electronics
converter, then $\theta_{i}$ is the voltage phase angle, $m_{i}$ is the power
measurement time constant or the virtual inertia through a controlled device,
and $d_{i}$ is the droop control coefficient. The joint state-space
representation of the network is

$\displaystyle\begin{bmatrix}\dot{\theta}\\\
\dot{\omega}\end{bmatrix}=\underbrace{\begin{bmatrix}\operatorname{\mathbf{0}}&I\\\
-M^{-1}\operatorname{\mathsf{L}}&-M^{-1}D\end{bmatrix}}_{A}\begin{bmatrix}\theta\\\
\omega\end{bmatrix}+\underbrace{\begin{bmatrix}\operatorname{\mathbf{0}}\\\
M^{-1}\end{bmatrix}}_{B}u,$ (1)

where $x=(\theta;\omega)\in^{2n}$ corresponds to the stacked vector of angle
and frequency deviations at each node,
$M=\operatorname{diag}(m_{i})\in^{n\times n}$ is the diagonal matrix with
inertia coefficients, $D=\operatorname{diag}(d_{i})\in^{n\times n}$ is the
diagonal matrix with droop control coefficients, and
$\operatorname{\mathsf{L}}$ is the Laplacian of the weighted version of
$\operatorname{\mathcal{G}}$ whose adjacency matrix is
$\operatorname{\mathsf{A}}_{ij}=b_{ij}$, $i,j\in\\{1,\dots,n\\}$. One can
verify that $(A,B)$ is stabilizable.

The formulation (1) assumes that the inertia of the system remains constant
and makes sense in the traditional paradigm of power systems. However, in
scenarios with increasing penetration of renewables, the inertia of the
network may change over time. Hence, it is reasonable to incorporate the time
dependence in the inertia at each node. We do this by considering a switched-
affine system representation as in [23], where each mode corresponds to a
different value of the inertia. If $m\in\mathbb{Z}_{>0}$ is the number of
modes, the frequency dynamics are then given by

$\displaystyle\begin{bmatrix}\dot{\theta}\\\
\dot{\omega}\end{bmatrix}=\underbrace{\begin{bmatrix}\operatorname{\mathbf{0}}&I\\\
-M_{q(t)}^{-1}\operatorname{\mathsf{L}}&-M_{q(t)}^{-1}D\end{bmatrix}}_{A_{q(t)}}\begin{bmatrix}\theta\\\
\omega\end{bmatrix}+\underbrace{\begin{bmatrix}\operatorname{\mathbf{0}}\\\
M_{q(t)}^{-1}\end{bmatrix}}_{B_{q(t)}}u.$ (2)

Here, at time $t$, the system is in mode $q(t)\in\\{1,\dots,m\\}$ and
$M_{q}(t)$ denotes the inertia of the network in mode $q(t)$. The inertia at
time $t$ depends on the online generators and the connected power electronics
converters at that time. When convenient, we drop the argument $t$ and refer
to $q(t)$ as $q$.

Our goal is to design an optimal controller that brings the system (2) to the
origin from any initialization. Since we might not have knowledge of the
current operating mode at all times, our aim is to design a time-invariant
controller of the form

$\displaystyle u=Kx,$

that stabilizes (2), minimizes the state deviation, and optimizes the control
input required. For a fixed linear system, this is achievable using the
solution to the linear-quadratic control (LQR) problem. However, for the
switched system, this cannot be done unless the switching sequence is known
beforehand. Optimizing instead for all possible switching sequences quickly
becomes computationally intractable. Therefore, we follow an offline, data-
driven, imitation-based approach that balances the goals of optimality and
stability: the basic idea is to consider a set of candidate switching
sequences, solve a finite-horizon LQR problem for each of them, and finally
use the resulting trajectories as a training set to design a stabilizing
controller imitating the observed behavior.

##  III Data-Driven Controller Design

In this section, we carry out our approach to design a common stabilizing
time-invariant controller using training data generated for system (2) for a
variety of scenarios. We start by describing in Section III-A how the data is
generated via a finite-horizon LQR formulation. Then we provide in Section
III-B a least-squares formulation to learn the controller while guaranteeing
the stability of each mode $q\in\\{1,\dots,m\\}$. Since the stability of all
the modes is not sufficient to guarantee the stability of the switched system,
we generalize in Section III-C our treatment to the stabilization of the
switched system via a common Lyapunov function.

###  III-A Training Data from Optimal Input Trajectories

In order to generate the training data which would later be used to learn the
controller gain $K$, we solve $\operatorname{\mathcal{S}}\in\mathbb{Z}_{>0}$
instances of the finite-horizon LQR problem

$\displaystyle\min_{\operatorname{\mathbf{x}},\mathbf{u}}$
$\displaystyle\int\limits_{0}^{T}\big{(}x(t)^{\top}Qx(t)+u(t)^{\top}Ru(t)\big{)}\;dt$
(3) s.t. $\displaystyle x(0)=x_{0}$
$\displaystyle\dot{x}(t)=A_{q(t)}x(t)+B_{q(t)}u(t),\quad t\in[0,T],$

where $Q\succeq\operatorname{\mathbf{0}}\in^{2n\times 2n}$ penalizes state
deviations, $R\succ\operatorname{\mathbf{0}}\in^{n\times n}$ represents a cost
associated to the control action, $T>0$ is the time horizon, $x_{0}\in^{2n}$
is the initial state, and $\operatorname{\mathbf{x}}(t)\in^{2n}$ and
$\mathbf{u}(t)\in^{n}$ are the variables describing the optimal state and
input trajectories, respectively.

We generate $\operatorname{\mathcal{S}}$ scenarios by selecting different
initial conditions $x_{0}$ and switching sequences $q(t)$, with the pair
$(\operatorname{\mathbf{x}}^{k}(t),\mathbf{u}^{k}(t))$ denoting the training
data for scenario $k\in\\{1,\dots,\operatorname{\mathcal{S}}\\}$. The
scenarios provide data in the form of desirable trajectories for the
controller to imitate. The amount of information available to capture
optimality grows with the number of scenarios considered, at the cost of an
increasing computational effort to handle them. Also, the number of
trajectories by itself does not guarantee that the resulting controller is
stable. Instead, in our design formulations below, we make sure the stability
of the controller is guaranteed independently of the number of scenarios
considered. Regarding the selection of initial conditions for the scenarios,
since the frequency deviation is usually bounded for real systems, from a
practical viewpoint, rather than taking them to be uniformly distributed
throughout the state space, it makes sense to consider initial conditions
close to the origin.

###  III-B Simultaneous Stabilization of All Switching Modes

Here, we are interested in designing a learned time-invariant controller which
guarantees stability for each mode $q\in\\{1,\dots,m\\}$. Let
$\operatorname{\mathcal{H}}$ denote the set of Hurwitz matrices. Then the
controller design problem described above can be cast as an optimization of
the form

$\displaystyle\min_{K}$
$\displaystyle\sum\limits_{k=1}^{\operatorname{\mathcal{S}}}\int\limits_{0}^{T}\|\mathbf{u}^{k}(t)-K\operatorname{\mathbf{x}}^{k}(t)\|_{2}^{2}\;dt$
(4) s.t. $\displaystyle A_{q}+B_{q}K\in\operatorname{\mathcal{H}},\quad\forall
q.$

Since the set of Hurwitz matrices is not convex, (4) is non-convex. In fact,
finding a feasible solution of (4), also referred to as the _simultaneous
stabilization_ problem, is NP-hard for general system and input matrices, cf.
[24]. However, the matrices $\\{A_{q}\\}_{q=1}^{m}$ and
$\\{B_{q}\\}_{q=1}^{m}$ in our setup are not arbitrary, and indeed have a
well-defined structure. Specifically, the only quantity that specifies the
operating mode $q\in\\{1,\dots,m\\}$ is the inertia matrix $M_{q}$. Building
on this insight, we prove that the simultaneous stabilization problem (4) is
always feasible. Our proof is constructive and relies on identifying a
controller that stabilizes all the modes.

######  Proposition III.1.

_(Feasibility of the simultaneous stabilization data-driven problem for
individual modes):_ Problem (4) is always feasible.

###### Proof:

Let $K=[K_{1}\quad K_{2}]$, where $K_{1},K_{2}\in^{n\times n}$. Then from
equation (2), the closed-loop system matrix for mode $q\in\\{1,\dots,m\\}$ is
given by

$\displaystyle A_{q}+B_{q}K$
$\displaystyle=\begin{bmatrix}\operatorname{\mathbf{0}}&I\\\
-M_{q}^{-1}\operatorname{\mathsf{L}}&-M_{q}^{-1}D\end{bmatrix}+\begin{bmatrix}\operatorname{\mathbf{0}}\\\
M_{q}^{-1}\end{bmatrix}\begin{bmatrix}K_{1}&K_{2}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}\operatorname{\mathbf{0}}&I\\\
-M_{q}^{-1}(\operatorname{\mathsf{L}}-K_{1})&-M_{q}^{-1}(D-K_{2})\end{bmatrix}.$
(5)

Let us first consider the case when, in a given mode $q$, the inertia
coefficient of all the nodes is the same, and is given by
$m_{q}\in\mathbb{R}_{>0}$. Then we have $M_{q}=m_{q}I$. Choosing

$\displaystyle K_{1}=\operatorname{\mathsf{L}}-I\text{ and }K_{2}=D-I,$ (6)

the closed-loop system matrix (5) becomes

$\displaystyle A_{q}+B_{q}K$
$\displaystyle=\begin{bmatrix}\operatorname{\mathbf{0}}&I\\\
-1/m_{q}I&-1/m_{q}I\end{bmatrix}$
$\displaystyle=\underbrace{\begin{bmatrix}0&1\\\
-1/m_{q}&-1/m_{q}\end{bmatrix}}_{S_{q}}\otimes I.$ (7)

The eigenvalues of the $2\times 2$ matrix $S_{q}$ are negative for all
$m_{q}>0$. Hence, $A_{q}+B_{q}K\in\operatorname{\mathcal{H}}$ for all
$q\in\\{1,\dots,m\\}$ .

Next we consider the general case where each node $i\in\\{1,\dots,n\\}$ might
have a different inertia coefficient. Once again, choose $K_{1}$ and $K_{2}$
according to (6). The closed-loop system matrix (5) now takes the form

$\displaystyle A_{q}+B_{q}K=\begin{bmatrix}\operatorname{\mathbf{0}}&I\\\
-M_{q}^{-1}&-M_{q}^{-1}\end{bmatrix}$

For each mode $q\in\\{1,\dots,m\\}$, consider the Lyapunov function candidate
$V_{q}:^{2n}\to$

$\displaystyle
V_{q}=x^{\top}\underbrace{\begin{bmatrix}I&\operatorname{\mathbf{0}}\\\
\operatorname{\mathbf{0}}&M_{q}\end{bmatrix}}_{P_{q}}x.$

The Lie derivative of $V_{q}$ is given by

$\displaystyle\operatorname{\mathcal{L}}_{f}V(x)$
$\displaystyle=x^{\top}\big{(}(A_{q}+B_{q}K)^{\top}P_{q}+P_{q}(A_{q}+B_{q}K)\big{)}x$
$\displaystyle=x^{\top}\begin{bmatrix}\operatorname{\mathbf{0}}&\operatorname{\mathbf{0}}\\\
\operatorname{\mathbf{0}}&-2I\end{bmatrix}x\leq 0.$

This means that each mode $q\in\\{1,\dots,m\\}$ is stable, and the result
follows. ∎

The proof of Proposition III.1 considers first the case of equal inertia at
each node, and then generalizes the argument to the case of different inertia
at the nodes. Although establishing the feasibility of the simultaneous
stabilization problem (4) in the former case is a special case of the latter,
it is interesting to consider it separately as the eigenvalues of the closed-
loop system can be explicitly characterized.

######  Remark 1.

_(Distributed learned controller stabilizing all the modes):_ The proof of
Proposition III.1 is constructive and relies on identifying a (not necessarily
optimal) controller stabilizing all the modes. It is interesting to note that
the controller identified in (6) is distributed over
$\operatorname{\mathcal{G}}$, meaning that to implement it, each node
$i\in\\{1,\dots,n\\}$ needs to know just its angle and frequency, and the
angle of the nodes to which it is electrically connected. $\bullet$

###  III-C Simultaneous Stabilization of the Switched System

The controller resulting from the simultaneous stabilization problem (4) in
Section III-B guarantees the stability of each individual mode, but does not
guarantee the stability of the overall switched system (2) in general, cf.
[25]. To address this, here we reformulate the synthesis of the learned time-
invariant controller by specifying a common Lyapunov function as a certificate
of its correctness. Formally, the controller design problem takes now the form

$\displaystyle\min_{K,P}$
$\displaystyle\sum\limits_{k=1}^{\operatorname{\mathcal{S}}}\int\limits_{0}^{T}\|\mathbf{u}^{k}(t)-K\operatorname{\mathbf{x}}^{k}(t)\|_{2}^{2}\;dt$
(8) s.t.
$\displaystyle(A_{q}+B_{q}K)^{\top}\\!P\\!+\\!P(A_{q}+B_{q}K)\prec\operatorname{\mathbf{0}},\;\forall
q$ $\displaystyle P\succ\operatorname{\mathbf{0}}.$

In this formulation, we aim to find a common quadratic Lyapunov function given
by $V(x)=x^{\top}Px$. The first constraint in (8) ensures that the Lie
derivative of the Lyapunov function along the evolution of (2) remains
negative for each mode $q\in\\{1,\dots,m\\}$, thereby guaranteeing the
stability of the switched system. Note that the problem (8) is bilinear in the
decision variables $K$ and $P$ and, hence, nonconvex. The next result
establishes the feasibility of problem (8).

######  Theorem III.2.

_(Feasibility of the simultaneous stabilization data-driven problem for the
switched system):_ Problem (8) is always feasible.

###### Proof:

Using $X=P^{-1}$ and $K=YX^{-1}$, cf. [26], the constraints in problem (8) can
be equivalently written as

$\displaystyle
A_{q}X\\!+\\!XA_{q}^{\top}\\!+\\!B_{q}Y\\!+\\!Y^{\top}B_{q}^{\top}\prec\operatorname{\mathbf{0}},\;\forall
q$ (9a) $\displaystyle X\succ\operatorname{\mathbf{0}}.$ (9b)

With $X_{1},X_{2},X_{3},Y_{1},Y_{2}\in^{n\times n}$, let
$X=\begin{bmatrix}X_{1}&X_{2}\\\ X_{2}^{\top}&X_{3}\end{bmatrix}$ and
$Y=\begin{bmatrix}Y_{1}&Y_{2}\end{bmatrix}$. Then using the structure of
$\\{A_{q}\\}_{q=1}^{m}$ and $\\{B_{q}\\}_{q=1}^{m}$, constraint (9a) can be
rewritten as

$\displaystyle\begin{bmatrix}\operatorname{\mathbf{0}}&I\\\
-M_{q}^{-1}\operatorname{\mathsf{L}}&-M_{q}^{-1}D\end{bmatrix}\begin{bmatrix}X_{1}&X_{2}\\\
X_{2}^{\top}&X_{3}\end{bmatrix}+\begin{bmatrix}\operatorname{\mathbf{0}}\\\
M_{q}^{-1}\end{bmatrix}\begin{bmatrix}Y_{1}&Y_{2}\end{bmatrix}+$
$\displaystyle\begin{bmatrix}X_{1}&X_{2}\\\
X_{2}^{\top}&X_{3}\end{bmatrix}\begin{bmatrix}\operatorname{\mathbf{0}}&-\operatorname{\mathsf{L}}M_{q}^{-1}\\\
I&-DM_{q}^{-1}\end{bmatrix}+\begin{bmatrix}Y_{1}^{\top}\\\
Y_{2}^{\top}\end{bmatrix}\begin{bmatrix}\operatorname{\mathbf{0}}&M_{q}^{-1}\end{bmatrix}\prec\operatorname{\mathbf{0}},$

for all $q\in\\{1,\dots,m\\}$. Performing the matrix multiplications and using
the abbreviated notation

$\displaystyle Z_{q}$
$\displaystyle=-M_{q}^{-1}\operatorname{\mathsf{L}}X_{1}-M_{q}^{-1}DX_{2}^{\top}+M_{q}^{-1}Y_{1}$
$\displaystyle W_{q}$
$\displaystyle=M_{q}^{-1}\operatorname{\mathsf{L}}X_{2}+M_{q}^{-1}DX_{3},$

the inequality can be further rewritten as

$\displaystyle\begin{bmatrix}X_{2}^{\top}&X_{3}\\\
Z_{q}&-W_{q}+M_{q}^{-1}Y_{2}\end{bmatrix}+\begin{bmatrix}X_{2}&Z_{q}^{\top}\\\
X_{3}^{\top}&-W_{q}^{\top}+Y_{2}^{\top}M_{q}^{-1}\end{bmatrix}\prec\operatorname{\mathbf{0}},$

for all $q\in\\{1,\dots,m\\}$. Hence, (9a) is satisfied if the matrix

$\displaystyle\begin{bmatrix}-X_{2}-X_{2}^{\top}&-X_{3}-Z_{q}^{\top}\\\
-X_{3}^{\top}-Z_{q}&W_{q}+W_{q}^{\top}-M_{q}^{-1}Y_{2}-Y_{2}^{\top}M_{q}^{-1}\end{bmatrix}$

is positive definite for all $q\in\\{1,\dots,m\\}$. Using the Schur
complement, cf. [27], the positive definiteness condition (and hence (9a)) is
equivalent to

$\displaystyle-X_{2}-X_{2}^{\top}\succ$
$\displaystyle\operatorname{\mathbf{0}},$ (10a) $\displaystyle
W_{q}+W_{q}^{\top}-M_{q}^{-1}Y_{2}-Y_{2}^{\top}M_{q}^{-1}\succ$
$\displaystyle\operatorname{\mathbf{0}},$ (10b) $\displaystyle
W_{q}+W_{q}^{\top}-M_{q}^{-1}Y_{2}-Y_{2}^{\top}M_{q}^{-1}+\qquad$
$\displaystyle(X_{3}^{\top}+Z_{q})(X_{2}+X_{2}^{\top})^{-1}(X_{3}+Z_{q}^{\top})\succ$
$\displaystyle\operatorname{\mathbf{0}},$ (10c)

for all $q\in\\{1,\dots,m\\}$. Now, choose $X\succ\operatorname{\mathbf{0}}$
satisfying (10a). Then, since $W_{q}$ and $Z_{q}$ are independent of $Y_{2}$,
there exists $N\prec\operatorname{\mathbf{0}}\in^{n\times n}$, independent of
$Y_{2}$ too, such that for all $q\in\\{1,\dots,m\\}$,

$\displaystyle W_{q}+W_{q}^{\top}-N\succ$
$\displaystyle\operatorname{\mathbf{0}},$ (11a)
$\displaystyle(X_{3}^{\top}+Z_{q})(X_{2}+X_{2}^{\top})^{-1}(X_{3}+Z_{q}^{\top})-N\succ$
$\displaystyle\operatorname{\mathbf{0}}.$ (11b)

Finally, using (10b)-(10c) and (11) along with the fact that
$N\prec\operatorname{\mathbf{0}}$, it suffices to show that there exists
$Y_{2}$ such that
$2N-M_{q}^{-1}Y_{2}-Y_{2}^{\top}M_{q}^{-1}\succ\operatorname{\mathbf{0}}$, for
all $q\in\\{1,\dots,m\\}$. Let $\overline{M}$ denote the matrix obtained after
taking the entry wise maximum of the inertia coefficient matrix at all nodes.
Then the above inequality is satisfied if $Y_{2}$ is chosen such that

$\displaystyle Y_{2}\prec N\overline{M},$ (12)

completing the proof. ∎

Proposition III.1 can be considered as a special case of Theorem III.2.
However, the results differ in their proof methodologies. The proof of
Proposition III.1 provides an explicit expression for a feasible controller,
which in addition is distributed over $\operatorname{\mathcal{G}}$. This,
however, does not mean that the optimizer of (4) is distributed (although it
does imply that one can look for solutions of (4) among controllers that are
distributed over $\operatorname{\mathcal{G}}$). Instead, the proof of Theorem
III.2 identifies an ordered sequence of steps that lead to the identification
of a controller stabilizing the switched system. In principle, there is no
guarantee that the resulting controller will be distributed. The following
result shows that a distributed controller does in fact exist. More precisely,
it shows that there exists a controller that does not need communication even
with neighboring nodes (we term this special form of distributed controller as
_local_). The proof methodology leverages the freedom in choosing various
parameters in the proof of Theorem III.2.

######  Corollary III.3.

_(Local controller stabilizing the switched system):_ There exists a
controller of the form $u=D_{1}\theta+D_{2}\omega$, where
$D_{1},D_{2}\in^{n\times n}$ are diagonal matrices, satisfying the constraints
in problem (8).

###### Proof:

Following the proof of Theorem III.2, we are interested in identifying $X$ and
$Y$ satisfying (9). Let us choose $X_{1}=I,~{}X_{2}=-I$, and
$Y_{1}=\operatorname{\mathbf{0}}$. Then using the Schur complement, (9b) is
satisfied iff

$\displaystyle X_{3}-I\succ\operatorname{\mathbf{0}}.$ (13)

To satisfy (9a), once again, following the same steps as in the proof of
Theorem III.2, choose $N$ satisfying (11) and then, $Y_{2}$ as a diagonal
matrix satisfying (12). The controller $K$ is then given by
$K=\begin{bmatrix}\operatorname{\mathbf{0}}&Y_{2}\end{bmatrix}\begin{bmatrix}I&-I\\\
-I&X_{3}\end{bmatrix}^{-1}$. Using the formula for the inverse of a
partitioned matrix [28, Section 0.7.3],

$\displaystyle K=$
$\displaystyle\begin{bmatrix}\operatorname{\mathbf{0}}&Y_{2}\end{bmatrix}\begin{bmatrix}(I-X_{3}^{-1})^{-1}&(X_{3}-I)^{-1}\\\
(X_{3}-I)^{-1}&(X_{3}-I)^{-1}\end{bmatrix}$ $\displaystyle=$
$\displaystyle\begin{bmatrix}Y_{2}(X_{3}-I)^{-1}&Y_{2}(X_{3}-I)^{-1}\end{bmatrix}.$

Now if one chooses $X_{3}$ to be a diagonal matrix (making the controller
local) satisfying (12) and (13), then the resulting controller stabilizes the
switched system. ∎

Although this result guarantees the existence of a local stabilizing
controller, restricting the feasible set of (8) to controllers of that form
could significantly affect the optimal value of the objective function.
Motivated by Corollary III.3 and this observation, we propose a middle ground
that reformulates the optimization problem to promote sparsity in the learned
controller. Formally, following [29], let $\beta>0$ be a design parameter that
specifies the importance of promoting sparsity as compared to the original
objective function of matching the data provided by the sampled optimal
trajectories. Let
$\operatorname{\mathcal{E}}^{c}=\\{(i,j)\;|\;(i,j)\notin\operatorname{\mathcal{E}}\\}$,
denote the set of indices whose corresponding vertices are not neighbors in
$\operatorname{\mathcal{G}}$. The sparse-promotion controller design problem
takes the form

$\displaystyle\min_{K,P}$
$\displaystyle\\!\sum\limits_{k=1}^{\operatorname{\mathcal{S}}}\int\limits_{0}^{T}\|\mathbf{u}^{k}(t)-K\operatorname{\mathbf{x}}^{k}(t)\|_{2}^{2}\;dt+\beta\sum\limits_{(i,j)\in\operatorname{\mathcal{E}}^{c}}|K_{ij}|$
(14) s.t.
$\displaystyle(A_{q}+B_{q}K)^{\top}\\!P\\!+\\!P(A_{q}+B_{q}K)\prec\operatorname{\mathbf{0}},\;\forall
q$ $\displaystyle P\succ\operatorname{\mathbf{0}}.$

Since (8) is feasible by Theorem III.2, problem (14) is feasible too. To find
a local controller, one could also consider a modified version of (14) where
all the non-diagonal entries of $K$ are penalized.

##  IV Simulations

In this section, we demonstrate the effectiveness of the proposed approach via
numerical experiments. We use the standard 12-bus 3-region network, shown in
Figure 1, that has also been used in [23, 1, 22].

Figure 1: The 12-bus 3-region network used in simulations.

We take $m=10$ and assume that at a given time $t$, the rotational inertia for
each node $i\in\\{1,\dots,n\\}$ is same. Hence, each mode
$q\in\\{1,\dots,m\\}$ of the hybrid system is given by one value of inertia in
the set {0.2, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 5, 9}. To generate the training
data-set, we use $Q=\begin{bmatrix}I&\operatorname{\mathbf{0}}\\\
\operatorname{\mathbf{0}}&10^{5}I\end{bmatrix}$ and $R=10I$. To implement (3),
we use its discrete-time counterpart

$\displaystyle\min_{\operatorname{\mathbf{x}},\mathbf{u}}$
$\displaystyle\sum\limits_{t=0}^{T}x_{t}^{\top}Qx_{t}+u_{t}^{\top}Ru_{t}$ (15)
s.t. $\displaystyle x_{t+1}=A^{d}_{q(t)}x_{t}+B^{d}_{q(t)}u_{t},\quad
t\in\\{0,\ldots,T-1\\},$

where $A^{d}_{q(t)}{}^{2n\times 2n}$ and $B^{d}_{q(t)}\in^{2n\times n}$ are
respectively, the state and input matrices of the discretized system using a
zero-order hold. We use a stepsize of $10^{-2}$ seconds and simulate 50
scenarios of (15), each for 50 time steps, using `cvx` [30]. The initial
conditions for all the scenarios (for both the angles and the frequencies) are
different, and drawn from a normal distribution with 0 mean and 0.1 variance.
Each scenario starts in mode 7 (3 seconds of inertia), and from there, based
on a uniform distribution draw, the inertia of the system can remain the same,
increase, or decrease every 2 time steps.

We design three sets of controllers: (a) _Optimal_ : To design the first
optimal and stable learned controller, we solve (8) using the BMI algorithm in
[31]. Since the algorithm requires a feasible initialization, we solve the
feasibility problem associated with the LMI constraints (9) using `cvx` to
find an initial point. (b) _Distributed_ : To design the second learned
controller, which is stable and sparse, we solve (14) for various values of
$\beta$, again using the algorithm in [31]. The controller turns out to be
distributed over $\operatorname{\mathcal{G}}$ for $\beta=100$. We observe that
instead of using the same initialization as in (a), taking the Optimal
controller as the initial point reduces the number of iterations to converge.
(c) _Unconstrained:_ The third learned controller that we design is based on
optimizing the objective function of fitting the controller to the sampled
data without any consideration of stability.

To compare the performance of the designed controllers, we display their
dynamical response for the same switching sequences. For each simulation, we
assume that the system starts in mode 10 (9 seconds of inertia), with an
initial frequency deviation of 0.05 Hz at each node, and can switch to any
other mode every 0.01 seconds. In Figure 2, we plot the frequency deviation at
node 1 for different switching (inertia) sequences.

Figure 2: Frequency deviation at node 1 for different switching sequences
using the learned controllers. Dashed vertical lines represent the switching
instances (every $10^{-2}$ second). Line styles correspond to different
switching sequences.

Frequency evolution with the unconstrained controller emphasizes the
importance of including the stability constraints for the switched system in
formulations (8) and (14). It is interesting to note that even though the
Optimal controller has an higher overshoot for all the switching sequences,
convergence is also faster. To further compare the Optimal and Distributed
controllers, we simulate the dynamics for 1 second in each mode, from an
initial frequency deviation of 0.15 Hz at every node. Table I provides the
total absolute value of the control input and the total absolute value of
frequency deviation for 3 fixed inertia modes $(q=1,5,10)$.

TABLE I: Performance metrics for the stable learned controllers under different inertia modes. Mode | Learned Controller | $\int\limits_{0}^{1}\sum\limits_{i=1}^{n}|u_{i}(t)|dt$ | $\int\limits_{0}^{1}\sum\limits_{i=1}^{n}|\omega_{i}(t)|dt$
---|---|---|---
1 | Optimal | $17.38$ | $0.004$
Distributed | $73.14$ | $0.027$
5 | Optimal | $163.66$ | $0.033$
Distributed | $627.97$ | $0.204$
10 | Optimal | $726.83$ | $0.140$
Distributed | $1308.10$ | $0.624$

As expected, the Optimal controller, which requires state information from all
the nodes, outperforms the Distributed controller. The mean of performance
differences taken over the 10 nodes is 62% for the cumulative control action,
and 79% for the cumulative frequency deviation. This trade-off in performance
comes with a saving of 90% in communication without compromising the system
stability.

##  V Conclusions and Future Work

We have presented a framework to synthesize data-driven controllers to
regulate the frequency of power networks under time-varying inertia. The
proposed learning-based design seeks to imitate, under suitable stability
constraints, optimal trajectories for different scenarios of changes in
inertia generated by finite-horizon LQR formulations. We establish that,
regardless of the inertia values, stabilizing learned controllers are
guaranteed to exist and are amenable to distributed implementation. Future
work will explore the design of efficient algorithms to identify distributed
controllers which take optimality with respect to the training data into
account and the extension of our approach to nonlinear AC power dynamics.

## References

  * [1] P. Kundur, _Power System Stability and Control_. McGraw-Hill, 1994.
  * [2] F. Dörfler, J. W. Simpson-Porco, and F. Bullo, “Breaking the hierarchy: Distributed control & economic optimality in microgrids,” _IEEE Transactions on Control of Network Systems_ , vol. 3, no. 3, pp. 241–253, 2016\.
  * [3] E. Ela, M. Milligan, and B. Kirby, “Operating reserves and variable generation,” National Renewable Energy Laboratory, Tech. Rep., Aug 2011.
  * [4] A. Ulbig, T. S. Borsche, and G. Andersson, “Impact of low rotational inertia on power system stability and operation,” _IFAC Proceedings Volumes_ , vol. 47, no. 3, pp. 7290–7297, 2014.
  * [5] “Update report– Black system event in South Australia on 28 September 2016,” Australian Energy Market Operator, Tech. Rep., October 2016.
  * [6] “Frequency stability evaluation criteria for the synchronous zone of continental Europe,” European Network of Transmission System Operators for Electricity, Tech. Rep., March 2016.
  * [7] “Future ancillary services in ERCOT,” Electricity Reliability Council of Texas, Tech. Rep., 2013.
  * [8] T. S. Borsche, T. Liu, and D. J. Hill, “Effects of rotational inertia on power system damping and frequency transients,” in _IEEE Conf. on Decision and Control_ , Osaka, Japan, 2015, pp. 5940–5946.
  * [9] G. S. Misyris, S. Chatzivasileiadis, and T. Weckesser, “Robust frequency control for varying inertia power systems,” in _IEEE PES Innovative Smart Grid Technologies Conference Europe_ , Sarajevo, Bosnia and Herzegovina, 2018, pp. 1–6.
  * [10] A. Kasis, S. Timotheou, and M. Polycarpou, “Stability of power networks with time-varying inertia,” in _IEEE Conf. on Decision and Control_ , Austin, TX, 2021, pp. 2788–2793.
  * [11] P. Hidalgo-Gonzalez, R. Henriquez-Auba, D. S. Callaway, and C. J. Tomlin, “Frequency regulation using data-driven controllers in power grids with variable inertia due to renewable energy,” in _IEEE PES General Meeting_ , Atlanta, GA, August 2019, pp. 1–5.
  * [12] ——, “Frequency regulation using sparse learned controllers in power grids with variable inertia due to renewable energy,” in _IEEE Conf. on Decision and Control_ , Nice, France, Dec 2019, pp. 3253–3259.
  * [13] R. Dobbe, P. Hidalgo-Gonzalez, S. Karagiannopoulos, R. Henriquez-Auba, G. Hug, D. S. Callaway, and C. J. Tomlin, “Learning to control in power systems: Design and analysis guidelines for concrete safety problems,” _Electric Power Systems Research_ , vol. 189, p. 106615, 2020.
  * [14] W. Wang, N. Yu, J. Shi, and Y. Gao, “Volt-var control in power distribution systems with deep reinforcement learning,” in _IEEE International Conference on Communications, Control, and Computing Technologies for Smart Grids_ , Beijing, China, October 2019, pp. 1–7.
  * [15] D. Ye, M. Zhang, and D. Sutanto, “A hybrid multiagent framework with Q-learning for power grid systems restoration,” _IEEE Transactions on Power Systems_ , vol. 26, no. 4, pp. 2434–2441, Nov 2011.
  * [16] M. H. Khooban and M. Gheisarnejad, “A novel deep reinforcement learning controller based type-II fuzzy system: Frequency regulation in microgrids,” _IEEE Transactions on Emerging Topics in Computational Intelligence_ , vol. 5, no. 4, pp. 689–699, 2021.
  * [17] M. Glavic, R. Fonteneau, and D. Ernst, “Reinforcement learning for electric power system decision and control: Past considerations and perspectives,” _IFAC-PapersOnLine_ , vol. 50, no. 1, pp. 6918 – 6927, 2017.
  * [18] D. Ernst, M. Glavic, and L. Wehenkel, “Power systems stability control: Reinforcement learning framework,” _IEEE Transactions on Power Systems_ , vol. 19, no. 1, pp. 427–435, Feb 2004.
  * [19] W. Cui and B. Zhang, “Reinforcement learning for optimal primary frequency control from inverter-based resources: A Lyapunov approach,” _Arxiv preprint arXiv:2009.05654_ , 2020.
  * [20] M. Jin and J. Lavaei, “Stability-certified reinforcement learning: A control-theoretic perspective,” _IEEE Access_ , vol. 8, pp. 229 086–229 100, 2020.
  * [21] C. D. Godsil and G. F. Royle, _Algebraic Graph Theory_ , ser. Graduate Texts in Mathematics. Springer, 2001, vol. 207.
  * [22] B. K. Poolla, S. Bolognani, and F. Dorfler, “Optimal placement of virtual inertia in power grids,” _IEEE Transactions on Automatic Control_ , vol. 62, no. 12, pp. 6209–6220, 2017.
  * [23] P. Hidalgo-Gonzalez, D. S. Callaway, R. Dobbe, R. Henriquez-Auba, and C. J. Tomlin, “Frequency regulation in hybrid power dynamics with variable and low inertia due to renewable energy,” in _IEEE Conf. on Decision and Control_ , Miami Beach, FL, 2018, pp. 1592–1597.
  * [24] V. Blondel and J. N. Tsitsiklis, “NP-hardness of some linear control design problems,” _SIAM Journal on Control and Optimization_ , vol. 35, no. 6, pp. 2118–2127, 1997.
  * [25] M. S. Branicky, “Multiple Lyapunov functions and other analysis tools for switched and hybrid systems,” _IEEE Transactions on Automatic Control_ , vol. 43, no. 4, pp. 475–482, 1998.
  * [26] G. E. Dullerud and F. Paganini, _A Course in Robust Control Theory_ , ser. Texts in Applied Mathematics. Springer, 2000, no. 36.
  * [27] S. Boyd and L. Vandenberghe, _Convex Optimization_. Cambridge University Press, 2004.
  * [28] R. A. Horn and C. R. Johnson, _Matrix Analysis_. Cambridge University Press, 1985.
  * [29] N. K. Dhingra, M. Jovanović, and Z.-Q. Luo, “An ADMM algorithm for optimal sensor and actuator selection,” in _IEEE Conf. on Decision and Control_ , 2014, pp. 4039–4044.
  * [30] M. Grant and S. Boyd, “CVX: Matlab software for disciplined convex programming, version 2.1,” Mar. 2014, available at http://cvxr.com/cvx.
  * [31] Q. Tran Dinh, S. Gumussoy, W. Michiels, and M. Diehl, “Combining convex–concave decompositions and linearization approaches for solving BMIs, with application to static output feedback,” _IEEE Transactions on Automatic Control_ , vol. 57, no. 6, pp. 1377–1390, 2012.


– Under Review 2023 Full Paper – MIDL 2023 submission Ricardo Mokhtari1
<EMAIL_ADDRESS>
Azam Hamidinekoo1<EMAIL_ADDRESS>
Daniel Sutton1<EMAIL_ADDRESS>
Arthur Lewis1<EMAIL_ADDRESS>
Bastian R. Angermann2<EMAIL_ADDRESS>
Ulf Gehrmann2<EMAIL_ADDRESS>
Pål Lundin3<EMAIL_ADDRESS>
Hibret Adissu4<EMAIL_ADDRESS>
Junmei Cairns2<EMAIL_ADDRESS>
Jessica Neisen2<EMAIL_ADDRESS>
Emon Khan5<EMAIL_ADDRESS>
Daniel Marks3<EMAIL_ADDRESS>
Nia Khachapuridze2,3<EMAIL_ADDRESS>
Talha Qaiser1<EMAIL_ADDRESS>
Nikolay Burlutskiy1<EMAIL_ADDRESS>
1 Imaging and Data Analytics, AstraZeneca, Cambridge, UK
2 Translational Science and Experimental Medicine, AstraZeneca, Gothenburg,
Sweden
3 Early Clinical Development, AstraZeneca, Cambridge, UK
4 Respiratory and Immunology Safety Pathology, AstraZeneca, Gaithersburg, USA
5 Late Respiratory and Immunology, AstraZeneca, Cambridge, UK

# Interpretable histopathology-based prediction of disease relevant features
in Inflammatory Bowel Disease biopsies using weakly-supervised deep learning

###### Abstract

Crohn’s Disease (CD) and Ulcerative Colitis (UC) are the two main Inflammatory
Bowel Disease (IBD) types. We developed interpretable deep learning models to
identify histological disease features for both CD and UC using only
endoscopic labels. We explored fine-tuning and end-to-end training of two
state-of-the-art self-supervised models for predicting three different
endoscopic categories (i) CD vs UC (AUC=0.87), (ii) normal vs lesional
(AUC=0.81), (iii) low vs high disease severity score (AUC=0.80). With the
support of a pathologist, we explored the relationship between endoscopic
labels, model predictions and histological evaluations qualitatively and
quantitatively and identified cases where the pathologist’s descriptions of
inflammation were consistent with regions of high attention. In parallel, we
used a model trained on the Colon Nuclei Identification and Counting (CoNIC)
dataset to predict and explore 6 cell populations. We observed consistency
between areas enriched with the predicted immune cells in biopsies and the
pathologist’s feedback on the attention maps. Finally, we identified several
cell level features indicative of disease severity in CD and UC. These models
can enhance our understanding about the pathology behind IBD and can shape our
strategies for patient stratification in clinical trials.

###### keywords:

Weakly supervised learning, self-supervised learning, attention maps, IBD.

††editors: Under Review for MIDL 2023

##  1 Introduction

Inflammatory bowel diseases (IBD) are chronic, relapsing-remitting
inflammatory disorders of the gastrointestinal (GI) tract. Crohn’s disease
(CD) and Ulcerative colitis (UC) represent the two main types of IBD. Both
result from a complex interplay of several factors, including abnormal immune
responses, genetics, microbiome and environmental triggers [Baumgart and
Carding (2007)]. CD can affect any region of the gut, and is characterised by
segmental mucosal ulceration, transmural inflammation, fissures, fibrosis, and
stricture formation. In contrast, UC principally involves the colon, and
manifests as continuous mucosal inflammation extending from the rectum
proximally, with variable extent, and a more superficial inflammatory
infiltrate. CD and UC share similar symptoms, but are pathophysiologically
distinct diseases [Shanahan (2001)].

IBD is diagnosed based on a combination of clinical presentation,
radiographic, endoscopic and histopathological findings [DeRoche et al.
(2014)]. Defining the extent and severity of inflammation in IBD can influence
treatment decisions and support prediction of a patient’s prognosis. While
endoscopic evaluation assesses the macroscopic tissue, histopathological
evaluation assesses the microscopic tissue and is typically carried out
through a trained pathologist’s visual inspection of a Haematoxylin & Eosin
(H&E)-stained tissue, digitised into a whole slide image (WSI). While there is
a strong correlation between endoscopic and histopathological assessment
[Irani et al. (2018)], the relationship between these two data modalities
(especially the relationship of endoscopic scores to histopathology) is not
completely understood [Lemmens et al. (2013)] and can potentially be improved
through the development of interpretable machine learning models that predict
endoscopic categories from H&E stained WSIs coupled with the interpretation of
model predictions by a pathologist.

In this paper we demonstrate that applying self-supervised learning coupled
with weakly-supervised learning to H&E-stained IBD biopsies can accurately
distinguish disease type, macroscopic tissue appearance, and endoscopic
scores. Specifically, we trained two recent state-of-the-art architectures,
including Dual-Stream Multiple Instance Learning (DSMIL) [Li et al. (2021a)]
and Hierarchical Image Pyramid Transformer (HIPT) [Chen et al. (2022)] on the
large SPARC IBD dataset containing 1394 WSIs, and explored two training
strategies - fine-tuning and end-to-end (E2E) training. We explored the
relationship between endoscopic prediction and histology through pathologist
collaboration, where we found that the high attention regions identified by
the models were confirmed, qualitatively, to contain epithelial and stromal
morphological/structural features consistent with inflammation. We further
validated these models by leveraging a model trained on the publicly available
Colon Nuclei Identification and Counting (CoNIC) dataset [Graham et al.
(2021)], which can segment and classify 6 types of cells. These predictions
were compared with the attention maps produced by the weakly-supervised
models.

These models have applications in clinical trials since they can speed up the
workflow of pathologists by ranking WSIs by disease severity so that
pathologists can prioritise the most severe tissues first as well as drawing
their attention to the most diseased regions in the WSI. To our knowledge,
this is the first work on exploring and validating weakly-supervised methods
to associate endoscopic appearance with H&E morphology for large-scale IBD
biopsies. Our code is available online at https://github.com/AstraZeneca/ibd-
interprethttps://github.com/AstraZeneca/ibd-interpret.

##  2 Related work

WSI classification using only slide-level labels is commonly framed as a
multiple instance learning (MIL) problem [Augustine (2022)] in which image
patches are mapped to fixed-length embeddings followed by aggregation by an
operator, for example max-pooling. However, these MIL methods suffer from
several limitations that have been recently addressed by the research
community. For example, MIL methods are prone to misclassifications in the
case of an unbalanced number of positive instances when using a simple
aggregation operation such as max-pooling [Li et al. (2021a)]. This can be
addressed by utilising attention-based aggregation where each instance is
given a specific attention weight during training [Tomita et al. (2019)].
Additionally, MIL methods can ignore the correlation among different instances
across the WSI, which can be mitigated by using transformers to consider
morphological and spatial information [Shao et al. (2021)] as well as
generating attention maps for improved interpretability. Attention map
interpretability was further improved by incorporating non-local attention
into a dual-stream MIL (DSMIL) architecture to calculate an attention weight
for each patch [Li et al. (2021a)]. MIL models often use fixed patch-based
features extracted by a CNN or only features extracted by a fine-tuned model
since end-to-end training is expensive and time-demanding for large slides.
These sub-optimal extracted features can be improved by incorporating a multi-
resolution feature fusion mechanism to leverage varied-size tissue features
such as glands vs cells [Li et al. (2021b)]. Recently, HIPT [Chen et al.
(2022)] extended multi-resolution feature fusion by introducing hierarchical
pre-training, enabling a bottom-up aggregation of tissue features from cells
to tissue morphology.

Following these advancements, we leverage DSMIL [Li et al. (2021a)] and HIPT
[Chen et al. (2022)] due to their superior performance over the aforementioned
methods on large public H&E classification datasets such as The Cancer Genome
Atlas (TCGA). DSMIL achieves strong performance by combining self-supervised
pre-training using SimCLR [Chen et al. (2020)] with a dual-stream attention-
based MIL aggregator. Since both local and global contexts are important for
WSI classification, DSMIL takes patches extracted at multiple magnifications
with the same resolution as input. This is disadvantageous since fine details
are lost at low magnification. By contrast, HIPT extracts patches at a fixed
high magnification but with a relatively large resolution, such that fine and
coarse-grained details are captured at the same resolution.

In this paper, we predict clinically relevant IBD endoscopic categories
including disease types, macroscopic appearances, and endoscopic scores from
H&E biopsies which, to our knowledge, has not been addressed previously. Prior
work mostly focused on electronic health records [Reddy et al. (2019)],
genomics, metagenomics [Abbas et al. (2019)], biological and clinical
parameters for predicting endoscopic scores in UC patients [Popa et al.
(2020)]. Recently a fully supervised approach using H&E images to predict
remission and the Nancy Histological Index scores for UC was investigated by
Najdawi et al. (2022). In contrast, we consider both CD and UC and use weakly
supervised methods with macroscopic patient level labels, without any detailed
pathologist annotations at the microscopic level, which are inexpensive
relative to the supervised method.

##  3 Methodology

Figure 1: Overview of the implemented pipeline.

The proposed methodology is summarised in Figure 1. First, we used an image
quality control (QC) algorithm and extracted patches for training DSMIL and
HIPT to predict 3 clinically relevant endoscopic categories and generate
visual attention maps. A cell classification model was also applied to all
slides and the individual cell predictions were aggregated into heatmaps for
comparison with the weakly-supervised attention maps. Finally, attention maps
and cell heatmaps were visualised in HALO Link for discussion with the
pathologist. More details on these steps are in the subsections below.

###  3.1 Patch Extraction and Quality Control (QC)

To process gigapixel WSIs using deep learning models, we first split each WSI
into many small, non-overlapping image patches using the Histolab Marcolini et
al. (2022) and CLAM Lu et al. (2021) Python packages. We trained a QC model
based on DenseNet Huang et al. (2017) to automatically identify and remove
regions from the WSIs with imaging and/or tissue artefacts such as
overstaining, tissue folds, debris, out-of-focus areas and variations in
contrast and hue markings. The QC model provided readouts at the slide-level,
identifying how much of the slide was rejected as a proportion of the tissue.
These readouts were used to exclude slides containing large areas of rejected
tissue (e.g. $>$50%). The QC model was also integrated into the patch
extraction pipeline to save individual patches containing $>$50% accepted
tissue (see Appendix A).

###  3.2 Self-Supervised Pre-training and Weakly-Supervised Classification

In the following sections, we briefly outline how we distinguish self-
supervised pre-training and weakly-supervised classification in DSMIL and HIPT
- further details are provided in Section 4 and Appendix A.

Self-supervised Pre-training. We only consider the self-supervised
representation learning components of DSMIL and HIPT. For DSMIL, this involves
the SimCLR component, while for HIPT, this involves the ViT${}_{16}-256$ and
ViT${}_{256}-4096$ components, which learn to output embeddings at the patch
level. These self-supervised components, once trained, are used to generate
patch embeddings for all WSIs, which are subsequently used for training the
weakly-supervised classification components.

Weakly-supervised classification. We only consider the weakly-supervised
components of DSMIL and HIPT for predicting 3 endoscopic categories. For DSMIL
the attention-based MIL aggregator, and for HIPT the ViT${}_{4096}-\text{WSI}$
component were considered.

###  3.3 Interpretation of Model Predictions

Following weakly-supervised training, we generated visual attention maps,
which facilitated optimal collaboration with the pathologist to interrogate
the models’ predictions.

DSMIL attention maps are visualised using the per-patch attention weights of
MIL aggregator as a heatmap overlay on the WSI Li et al. (2021a). A weight
value close to 1.0 indicates that the patch contributes heavily to the final
prediction compared with a patch that has a score close to 0.0. We generated
these attention maps for all WSIs in our dataset.

HIPT attention maps are a natural part of HIPT due to its transformer backbone
Chen et al. (2022). The attention maps in the original paper demonstrated that
they can highlight unique cancer-relevant tissue morphologies. Therefore, we
also sought to understand what IBD-relevant tissue morphologies HIPT could
learn.

HALO Link Integration was used to make the attention map review process
interactive and straightforward for the pathologist. HALO is an image analysis
platform for quantitative tissue analysis in digital pathology. We used HALO
Link for sharing the WSIs and our predictions with the pathologists. The
pathologist was not aware of any slide level endoscopic labels prior to
assessment. For each slide, we asked the pathologist to describe the critical
histopathological features for each slide including identification of
inflammatory regions cellular composition, and morphological/structural
features of the epithelium and stroma. We used HALO
Link111https://indicalab.com/halolink/ to compare the attention maps and the
pathologist comments side by side to interpret features learnt by the models.

A cell prediction model for detection and classification of six cell types
trained on the publicly available annotated CoNIC H&E dataset Graham et al.
(2021) was used for understanding 6 cell populations. The CoNIC dataset was
part of a grand medical challenge with the aim of predicting six types of
cells on H&E slides, including epithelial cells, neutrophils, lymphocytes,
plasma cells, connective cells, and eosinophils. We aggregated the cell
predictions into statistics per patch and then visualised these as the
heatmaps (Appendix E). The heatmaps were compared with the attention maps
produced by DSMIL and HIPT and several representative examples were reviewed
by the pathologist. Finally, several human interpretable features (HIFs) were
calculated from the readouts and correlated with the endoscopic scores for
both CD and UC.

##  4 Experimental Setup

We included 1394 H&E-stained biopsies from 418 CD and 218 UC patients enrolled
in a multi-centred longitudinal Study of a Prospective Adult Research Cohort
with IBD (SPARC IBD). The samples reside in the IBD Plexus database, provided
with informed consent by the Crohn’s and Colitis Foundation Raffals et al.
(2022). The total tissue area in 1394 biopsies is $\approx$7500µm2 whereas the
TCGA-lung resection dataset used to train DSMIL in the original work Li et al.
(2021a) has a total area of $\approx$65000µm2 across 1054 slides, making SPARC
IBD a relatively challenging dataset for classification.

We used 3 clinically-relevant labels from the SPARC IBD dataset that were
acquired from endoscopy to define weakly-supervised classification tasks.
These included: i) disease diagnosis (num. WSIs: CD=903 vs UC=491), ii)
macroscopic tissue appearance (num. WSIs: normal=922 vs lesional=472), and
iii) CD and UC-specific endoscopic severity scores (num. WSIs: low (CD=714,
UC=134) vs high (CD=281, UC=131) score) - see Appendix A for more details on
the exact splits. We used 4 NVIDIA A100 GPUs for training.

The QC model was trained on representative subset of 19 WSIs. We reserved a
set of 10 WSIs for testing on which the Dice score of the QC model in
segmenting background, good tissue and artefacts was 0.730. The QC model was
then applied to all SPARC IBD WSIs. We excluded 11 slides with $>$50% rejected
tissue and then DSMIL and HIPT were trained on patches extracted at 40x
(Appendix A).

In training, we initially used DSMIL and HIPT models, which were pre-trained
on the TCGA dataset and fine-tuned on our dataset, in order to assess the
transferability of knowledge from TCGA cancer resections to IBD biopsies. In
fine-tuning, the weights of the self-supervised components of both models
(SimCLR in DSMIL and ViT${}_{16}-256$ & ViT${}_{256}-4096$ in HIPT) were
frozen and only the weakly-supervised classification components (MIL
aggregator in DSMIL and ViT${}_{4096}-\text{WSI}$ in HIPT) were trained to
predict the weakly-supervised classification tasks in SPARC IBD. We labelled
the models that were fine-tuned on SPARC IBD as DSMIL-F and HIPT-F.

We also explored training DSMIL and HIPT end-to-end (E2E) on SPARC IBD. The
self-supervised and weakly-supervised components of both models were trained
on SPARC IBD from scratch to compare with the performance of fine-tuned models
mentioned previously. We labeled the models that were trained end-to-end on
SPARC IBD as DSMIL-E2E and HIPT-E2E (see Appendix A). In self-supervised and
weakly-supervised training, we closely followed the training settings in the
original papers Li et al. (2021a); Chen et al. (2022).

In weakly-supervised training, we performed 5-fold cross validation, with an
80:20 split in all three prediction categories stratified on patients and
ensuring that distributions of disease diagnosis, macroscopic appearance,
biopsy location and the target label were consistent across all training and
testing splits (more details in Appendix A). The same cross-validation splits
were used for both DSMIL and HIPT models in all experiments. We used area
under the receiver-operator characteristic curve (AUROC) to measure predictive
performance. For further evaluation, the visual attention maps for all slides
were generated and compared with the cell prediction heatmaps from the CoNIC
model. To obtain a quantitative evaluation within reasonable pathologist
effort, we sorted the WSIs by tissue size and selected 10 WSIs from CD and UC
patients, each containing 5 normal and 5 lesional.

##  5 Results and Discussion

Figure 2: H&E image, Attention maps, cell level predictions and pathologist
comments for UC patient with lesional macroscopic appearance and high
endoscopic score.

tab:auc_performance  Model Disease
Diagnosis Macroscopic
Appearance Endoscopic Score (CD) Endoscopic Score (UC) DSMIL-F 0.656±0.007
0.522±0.008 0.528±0.016 0.592±0.007 DSMIL-E2E 0.692±0.010 0.750±0.006
0.740±0.009 0.634±0.017 HIPT-F 0.825±0.017 0.780±0.012 0.766±0.026 0.788±0.034
HIPT-E2E 0.865±0.019 0.814±0.008 0.786±0.017 0.802±0.014

Table 1: Mean AUROC ± standard error (5-fold cross validation) of trained
models across different prediction tasks.

tab:extra_experiments  Model All data Just CD Just UC Just Ileum Just Colon
DSMIL-E2E 0.750±0.006 0.689±0.010 0.808±0.013 0.586±0.010 0.775±0.019 HIPT-E2E
0.814±0.008 0.804±0.010 0.837±0.015 0.739±0.032 0.823±0.036

Table 2: Performance of DSMIL-E2E and HIPT-E2E trained on clinically-relevant
subsets of SPARC IBD for the task of predicting macroscopic appearance.

Model Performance. In Table LABEL:tab:auc_performance we compare the
performance of end-to-end and fine-tuned DSMIL and HIPT models for predicting
3 weakly-supervised tasks in SPARC IBD. HIPT-E2E significantly outperforms
DSMIL-E2E across all tasks (two-tailed t-test: diagnosis - p$<$0.0001,
macroscopic appearance - p$<$0.0005, endoscopic scores (CD) - p$<$0.05 and
(UC) - p$<$0.0001). It is likely that both the spatial patterns among cells
and the context of the tissue microenvironment are well captured by HIPT’s
transformer backbone, leading to improved performance. In addition, we found
that E2E training on SPARC IBD was the optimal strategy for both DSMIL and
HIPT training across all tasks. We suggest that the embeddings learned by E2E
models through self-supervised pre-training on SPARC IBD were more useful for
downstream prediction tasks than pre-trained embeddings from TCGA. This
difference can be attributed to the differences between IBD and cancer
morphologies as well as the difference between TCGA resections vs IBD
biopsies.

Pathologist Evaluation. Through the pathologist evaluation of 20 WSIs, there
were 8 TP, 7 TN, 4 FP and 1 FN when using the endoscopic labels whereas there
were 9 TP, 7 TN, 3 FP and 1 FN when using the pathologist’s evaluation as
ground truth. Therefore, while not perfectly related, endoscopic and
histological labels are strongly correlated, and the model was able to learn
histologically relevant features via supervised training on endoscopic labels.
Additionally, we obtained pathologist annotations of inflammatory regions for
1 WSI (Figure 2) and computed the Dice coefficient between the annotation and
the thresholded attention map, obtaining 0.625, suggesting that inflammatory
infiltrate was well localised in this WSI.

Cell Prediction Model. From the cell model’s predictions, we found that the
most indicative HIF in UC was the ratio of neutrophils to all cells in the
tissue (p=0.0007), while in CD it was the ratio of eosinophils to all cells in
the tissue (p=0.0002). The additional HIFs are listed in Appendix E. These
findings are reflected in Alhmoud et al. (2020). Qualitatively, there are
similarities between the cell prediction maps and the pathologist’s
description of inflammation (for example, Figure 5). However, the quantitative
performance at the cell level must be explored further. As this is a labor-
intensive and time-consuming task for pathologists, we suggest combining
weakly-supervised learning with supervised models is an alternative, more
efficient, approach.

Applications. The trained models can be used to automatically rank biopsies by
disease severity allowing pathologists to save time by prioritising more
severe biopsies. This is important since the majority of biopsies received by
pathologists are within normal limits. For biopsies outside of normal limits,
the attention maps can be used for guiding the pathologist during microscopic
evaluation of potentially malignant areas. We intend to quantitatively assess
the speed up of pathologist workflow in future works. For clinical
applications, we explored training DSMIL and HIPT using only CD or only UC
biopsies as well as exclusively ileum or colon biopsies for predicting
macroscopic appearance since disease features can present differently within
these subsets (Table LABEL:tab:extra_experiments). We found that model
performance can be improved by training on just UC and just colon subsets,
while performance decreases when training on just ileum, suggesting that this
is a more difficult task due to fewer ileum biopsies. We consider these types
of modifications critical in applying these models to clinical trials and will
be further explored in future works.

##  6 Conclusion and Future Work

We demonstrated that weakly supervised learning can produce accurate models
trained on H&E-stained biopsies with endoscopic labels only. We observed
potential in the publicly available CoNIC model and the attention maps.
However, further quantitative and qualitative assessment of their
generalisability is needed. Through reviewing the attention maps with
pathologists, we can potentially better understand the relationship between
histological and endoscopic labels. We plan to validate the trained SPARC IBD
models with the corresponding attention maps on external IBD datasets from
clincal trials.

We would like to thank AstraZeneca for sponsoring and supporting this research
and Crohn’s & Colitis Foundation for providing the SPARC IBD dataset.

## References

  * Abbas et al. (2019) Mostafa Abbas, John Matta, Thanh Le, Halima Bensmail, Tayo Obafemi-Ajayi, Vasant Honavar, and Yasser El-Manzalawy.  Biomarker discovery in inflammatory bowel diseases using network-based feature selection.  _PloS one_ , 14(11):e0225382, 2019.
  * Alhmoud et al. (2020) Tarik Alhmoud, Anas Gremida, Diego Colom Steele, Imaneh Fallahi, Wael Tuqan, Nina Nandy, Mahmoud Ismail, Barakat Aburajab Altamimi, Meng-Jun Xiong, Audra Kerwin, and David Martin.  Outcomes of inflammatory bowel disease in patients with eosinophil-predominant colonic inflammation.  _BMJ Open Gastroenterology_ , 16(7):1–7, 2020\.
  * Augustine (2022) Tanya N. Augustine.  Weakly-supervised deep learning models in computational pathology.  _eBioMedicine_ , 81:104117, 2022.
  * Baumgart and Carding (2007) Daniel C Baumgart and Simon R Carding.  Inflammatory bowel disease: cause and immunobiology.  _The Lancet_ , 369(9573):1627–1640, 2007.
  * Chen et al. (2022) Richard J Chen, Chengkuan Chen, Yicong Li, Tiffany Y Chen, Andrew D Trister, Rahul G Krishnan, and Faisal Mahmood.  Scaling vision transformers to gigapixel images via hierarchical self-supervised learning.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 16144–16155, 2022.
  * Chen et al. (2020) Ting Chen, Simon Kornblith, Mohammad Norouzi, and Geoffrey Hinton.  A simple framework for contrastive learning of visual representations.  In _International conference on machine learning_ , pages 1597–1607. PMLR, 2020.
  * DeRoche et al. (2014) Tom C DeRoche, Shu-Yuan Xiao, and Xiuli Liu.  Histological evaluation in ulcerative colitis.  _Gastroenterology report_ , 2(3):178–192, 2014\.
  * Graham et al. (2021) Simon Graham, Mostafa Jahanifar, Quoc Dang Vu, Giorgos Hadjigeorghiou, Thomas Leech, David Snead, Shan E Ahmed Raza, Fayyaz Minhas, and Nasir Rajpoot.  Conic: Colon nuclei identification and counting challenge 2022.  _arXiv preprint arXiv:2111.14485_ , 2021.
  * Huang et al. (2017) Gao Huang, Zhuang Liu, Laurens Van Der Maaten, and Kilian Q Weinberger.  Densely connected convolutional networks.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 4700–4708, 2017.
  * Irani et al. (2018) Nazneen R Irani, Lai Mun Wang, Gary S Collins, Satish Keshav, and Simon PL Travis.  Correlation between endoscopic and histological activity in ulcerative colitis using validated indices.  _Journal of Crohn’s and Colitis_ , 12(10):1151–1157, 2018.
  * Lemmens et al. (2013) Bart Lemmens, Ingrid Arijs, Gert Van Assche, Xavier Sagaert, Karel Geboes, Marc Ferrante, Paul Rutgeerts, Séverine Vermeire, and Gert De Hertogh.  Correlation between the endoscopic and histologic score in assessing the activity of ulcerative colitis.  _Inflammatory bowel diseases_ , 19(6):1194–1201, 2013.
  * Li et al. (2021a) Bin Li, Yin Li, and Kevin W Eliceiri.  Dual-stream multiple instance learning network for whole slide image classification with self-supervised contrastive learning.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , pages 14318–14328, 2021a.
  * Li et al. (2021b) Jiayun Li, Wenyuan Li, Anthony Sisk, Huihui Ye, W Dean Wallace, William Speier, and Corey W Arnold.  A multi-resolution model for histopathology image classification and localization with multiple instance learning.  _Computers in biology and medicine_ , 131:104253, 2021b.
  * Lu et al. (2021) Ming Y Lu, Drew FK Williamson, Tiffany Y Chen, Richard J Chen, Matteo Barbieri, and Faisal Mahmood.  Data-efficient and weakly supervised computational pathology on whole-slide images.  _Nature biomedical engineering_ , 5(6):555–570, 2021.
  * Marcolini et al. (2022) Alessia Marcolini, Ernesto Arbitrio, and Nicole Bussola.  Histolab, 2022.  URL https://github.com/histolab/histolab.
  * Najdawi et al. (2022) Fedaa Najdawi, Kathleen Sucipto, Pratik Mistry, Stephanie Hennek, Christina Jayson, Mary Lin, Darren Fahy, Shawn Kinsey, Ilan Wapinski, Andrew H Beck, et al.  Artificial intelligence enables quantitative assessment of ulcerative colitis histology.  _medRxiv_ , 2022.
  * Popa et al. (2020) Iolanda Valentina Popa, Alexandru Burlacu, Catalina Mihai, and Cristina Cijevschi Prelipcean.  A machine learning model accurately predicts ulcerative colitis activity at one year in patients treated with anti-tumour necrosis factor $\alpha$ agents.  _Medicina_ , 56(11):628, 2020.
  * Raffals et al. (2022) Laura E Raffals, Sumona Saha, Meenakshi Bewtra, Cecile Norris, Angela Dobes, Caren Heller, Sirimon O’Charoen, Tara Fehlmann, Sara Sweeney, Alandra Weaver, et al.  The development and initial findings of a study of a prospective adult research cohort with inflammatory bowel disease (sparc ibd).  _Inflammatory bowel diseases_ , 28(2):192–199, 2022.
  * Reddy et al. (2019) Bhargava K Reddy, Dursun Delen, and Rupesh K Agrawal.  Predicting and explaining inflammation in crohn’s disease patients using predictive analytics methods and electronic medical record data.  _Health informatics journal_ , 25(4):1201–1218, 2019.
  * Shanahan (2001) Fergus Shanahan.  Inflammatory bowel disease: immunodiagnostics, immunotherapeutics, and ecotherapeutics.  _Gastroenterology_ , 120(3):622–635, 2001.
  * Shao et al. (2021) Zhuchen Shao, Hao Bian, Yang Chen, Yifeng Wang, Jian Zhang, Xiangyang Ji, et al.  Transmil: Transformer based correlated multiple instance learning for whole slide image classification.  _Advances in Neural Information Processing Systems_ , 34:2136–2147, 2021.
  * Tomita et al. (2019) Naofumi Tomita, Behnaz Abdollahi, Jason Wei, Bing Ren, Arief Suriawinata, and Saeed Hassanpour.  Attention-based deep neural networks for detection of cancerous and precancerous esophagus tissue on histopathological slides.  _JAMA network open_ , 2(11):e1914645–e1914645, 2019.

##  Appendix A Supplementary Methods

###  A.1 QC and Patch Extraction

From applying our QC model to all 1394 WSIs, we found that 6, 11, and 40
slides had, respectively, $>$75, $>$50, and $>$25% rejected areas (Figure 3).
Following QC, patches were extracted at 40x for training DSMIL-E2E (224 by 244
resolution) and HIPT-E2E (4096 by 4096 resolution), following the methods of
Li et al. (2021a) and Chen et al. (2022). 2,054,901 patches were extracted for
DSMIL and 13,952 patches were extracted for HIPT.

Figure 3: QC summary statistics on SPARC IBD and QC result on an additional
slide

###  A.2 Cross-Validation

5-fold cross validation is performed on in all weakly-supervised
classification tasks in SPARC IBD. Splits are performed at the patient level,
maintaining the distributions of biopsy location, disease diagnosis and the
target label - see Figure 4 for the exact splits for all tasks. This was done
to ensure that models are tested on a representative distribution of biopsies.
The same cross validation splits are used for both DSMIL and HIPT experiments
to allow them to be compared. For predicting macroscopic appearance, the
“erosions/ulcers” and “inflammation” labels were combined into a “lesional”
class. For endoscopic score, the median of SES score for CD and modified Mayo
endoscopic score for UC were 0, so we treated the negative class as a score of
0 and positive class as score $>$0.

Figure 4: Proportion of patients (n=638) in original, train and test sets
across 5 cross validation folds and across biopsy location, disease diagnosis
and macroscopic appearance for all tasks. Proportions of patients in each
category are kept consistent across all train and test splits.

###  A.3 Model Training

The hyperparameters used in self-supervised pretraining and weakly-supervised
training for all DSMIL and HIPT models are shown in Table
LABEL:tab:dsmil_hipt_parameters. To follow closely the methods of Li et al.
(2021a) and Chen et al. (2022), we use the default parameters for both models
across the board, but for HIPT we found that self-supervised pretraining
converged after 30 epochs and in weakly-supervised learning we use only 1
transformer encoder layer. The same cross validation splits are used for both
DSMIL and HIPT in all experiments.

tab:dsmil_hipt_parameters  Model Self-supervised pretraining Weakly-supervised
training DSMIL-F Frozen from TCGA epochs: 200; embedding size: 512; learning
rate: 0.0002; weight decay: 0.005 DSMIL-E2E epochs: 100; learning rate: 0.001;
weight decay: 0.000001; batch size: 1024; arch: ResNet18 epochs: 200;
embedding size: 512; learning rate: 0.0002; weight decay: 0.005 HIPT-F Frozen
from TCGA epochs: 20; layers: 1; heads: 3; dropout: 0.25, learning rate:
0.0003 HIPT-E2E epochs: 30; weight decay: 0.04; learning rate: 0.0005, warmup
epochs: 10 epochs: 20; layers: 1; heads: 3; dropout: 0.25, learning rate:
0.0003

Table 3: Hyperparameters used to train DSMIL and HIPT models.

##  Appendix B Additional Experiments

Since DSMIL can take as input patches extracted at different magnifications,
we assessed the performance of DSMIL in predicting disease diagnosis with
different magnification datasets. These results are summarised in Table
LABEL:tab:dsmil_mag_experiments. In order to be consistent with HIPT we used
40x patches to train DSMIL-E2E for all main experiments (Table
LABEL:tab:auc_performance). However, we find that for predicting macroscopic
appearance, slightly improved performance can be achieved with lower
magnifications, although there is no statistically significant difference
between these models trained at different magnifications and HIPT-E2E still
significantly outperforms DSMIL-E2E at all magnifications.

tab:dsmil_mag_experiments  Dataset AUROC ± 1 SE 5x 0.645±0.008 10x 0.752±0.006
20x 0.758±0.006 40x 0.750±0.006

Table 4: Comparison of DSMIL-E2E trained on different magnifications in
predicting macroscopic appearance.

##  Appendix C Additional Attention Maps and Pathologist Feedback

Figure 5: H&E image, lymphocytes density by the cell level model, attention
maps from two models, at WSI level and cell level resolutions, and pathologist
comments.

##  Appendix D Misclassified Cases

Figure 6: Two slides classified as “normal” from macroscopic appearance but
predicted as “lesional” by DSMIl-E2E. Both slides were confirmed by the
pathologist to contain inflammation and hence should be considered “lesional”.

##  Appendix E Supervised Cell Model

Figure 7: Six cell types predictions overlaid over two CoNIC patches.. CoNIC
examples have ground truth, annotations by pathologists available in the
dataset. Figure 8: Four immune cell density heatmaps overlaid over an
IBDplexus WSI. Figure 9: An example how human interpretable features (HIFs)
are calculated from six class cells predictions.


# From Stress Chains to Acoustic Emission

Ke Gao<EMAIL_ADDRESS>Geophysics, Los Alamos National Laboratory, Los
Alamos, NM, USA  Robert Guyer  Geophysics, Los Alamos National Laboratory, Los
Alamos, NM, USA  Physics, University of Nevada-Reno, NV, USA  Esteban Rougier
Geophysics, Los Alamos National Laboratory, Los Alamos, NM, USA  Christopher
X. Ren  Geophysics, Los Alamos National Laboratory, Los Alamos, NM, USA  Paul
A. Johnson  Geophysics, Los Alamos National Laboratory, Los Alamos, NM, USA

###### Abstract

A numerical scheme using the combined finite-discrete element methods (FDEM)
is employed to study a model of an earthquake system comprising a granular
layer embedded in a formation. When the formation is driven so as to shear the
granular layer, a system of stress chains emerges. The stress chains endow the
layer with resistance to shear and on failure launch broadcasts into the
formation. These broadcasts, received as acoustic emission, provide a remote
monitor of the state of the granular layer, of the earthquake system.

###### pacs:

123

Low frequency earthquakes (LFEs) Shelly  _et al._ (2006), non-volcanic tremor
Nadeau and Dolenc (2005); Beroza and Ide (2011); Rogers and Dragert (2003) and
acoustic emission Scholz (2002); Rouet-Leduc  _et al._ (2018) are examples of
weak seismic signals that may serve as harbingers of a major seismic event, an
earthquake Obara and Kato (2016). In a laboratory setting Kaproth and Marone
(2013), in which stick-slip events are simulated, acoustic emissions are
detected away from the stick-slip events Johnson  _et al._ (2013). Recent
machine learning and related studies Rouet-Leduc  _et al._ (2017); Brzinski
and Daniels (2018) of the acoustic emissions are able to use them to predict
location in the seismic cycle, i.e., the evolution of the acoustic emission as
a stick-slip scenario unfolds allows one to follow the scenario and anticipate
the subsequent earthquake. In all cases, field and laboratory, the signals
involved are sourced or detected in a volume remote from the volume that
spawns the earthquake. These findings bring to the fore the question of the
causal relationship between signals detected on passive, remote monitors and
the dynamics of the elastic structures that launch important seismic events.
In this paper we introduce a numerical model that lets us follow this
causality, i.e., examine/connect the dynamics in a granular system (fault
gouge), to signals detected on passive remote monitors. We find that stress
chains Peters  _et al._ (2005); Daniels and Hayman (2008) are the principals
in the dynamics of the granular system and that their dynamics are the source
of acoustic emissions. The numerical model necessarily combines discrete
element methods (DEM), to describe grain-grain interactions, and finite
element methods (FEM), to describe elasticity within grains and wave
propagation away from the granular system Morgan and Boettcher (1999);
Dorostkar  _et al._ (2017); Munjiza (2004).

The numerical model, two dimensional, comprises of a gouge layer of
approximately $50\times 1$ cm2 sandwiched between two $50\times 25$ cm2 pieces
of formation Geller  _et al._ (2015) (Fig. 1). The gouge layer is filled with
particles (disks) of two diameters, $1.2$ and $1.6$ mm, and each particle is
represented by 24 approximately equal size constant strain triangular FEM
elements. The particles are in contact with one another and with the formation
via a penalty function based contact interaction algorithm Munjiza  _et al._
(2011). The bottom edge of the lower formation piece is essentially rigid and
fixed in $x$, and upon which a normal force, $N$, is applied. The top edge of
the upper formation piece, also essentially rigid and fixed in $y$, is driven
uniformly at constant velocity $V_{0}$ in the $x$-direction. The material
density and elastic properties of gouge, formation, and top (bottom) edge of
the formation piece are recorded in Table I of the Supplementary Material.

Typically the state of the system is set by a choice of $V_{0}$ and $N$: in
the example here $V_{0}=0.50$ mm/sec and $N=28$ kPa. The basic output is the
shear and normal forces between the formation and gouge Gao  _et al._ (2018).
This is shown in Fig. 2 as a coefficient of friction (black), $\mu=$ (shear
force)$/$(normal force) vs time. This stress-time pattern has several
qualities of note. There is on average a linear stress-time relationship that
is punctuated by many small and a few large stress drops. Compressive
structures form in the gouge that push back against the effort of the upper
formation piece to drag the lower formation piece along Peters  _et al._
(2005); Daniels and Hayman (2008). The three green lines, all with the same
slope, illustrate the underlying spring-like character of the elasticity of
the gouge. Large stress drops occur irregularly in time. The stress values at
which large stress drops occur vary markedly, i.e., there is no single stress
at which the system fails.

To begin to look at the dynamics of the system we look at the behavior of the
formation immediately adjacent to the gouge. At the top of the gouge (T) (at
the bottom of the upper formation piece) and at the bottom of the gouge (B)
(at the top of the lower formation piece) there are 80 uniformly spaced sensor
pairs that detect the motion of the formation, i.e., the velocity of the
formation point with which the sensor is identified. In Fig. 2 we plot the
average $x$-velocity of 80 T-sensors (sensors just above the gouge and near
its center) as a function of time, $\overline{V}_{X}^{T}(t)$ (blue)

$\overline{V}_{X}^{T}(t)=\frac{1}{80}\sum_{n=1}^{80}V_{X}^{T}(t,n),$ (1)

where $V_{X}^{T}(t,n)$ is the $x$-component of the velocity of T sensor $n$ at
time $t$. The velocity $\overline{V}_{X}^{T}(t)$ has a background value of
approximately $0.25$ mm/sec ($V_{0}/2$) from which there are spikes in
velocity that coincide with sharp, large ($1$ mm/msec) stress changes. In Fig.
2 we also plot the average $x$-velocity of the 80 associated B-sensors
(sensors just below the gouge and near its center) as a function of time,
$\overline{V}_{X}^{B}(t)$ (red). The velocity $\overline{V}_{X}^{B}(t)$ also
has a background velocity of approximately $0.25$ mm/sec. From this background
there are spikes in $\overline{V}_{X}^{B}(t)$ that coincide with sharp, large
stress changes. The spikes in $\overline{V}_{X}^{B}(t)$ are opposite in sign
from the spikes in $\overline{V}_{X}^{T}(t)$. That is, attending sharp, large
stress changes are sharp, large velocity dipoles delivered to the formation.
We focus on these velocity dipoles. To characterize them we form the velocity
dipole field

$D(t,n)=\frac{V_{X}^{T}(t,n)-V_{X}^{B}(t,n)}{2}$ (2)

that has a value at each moment of time $t$ at each sensor $n$, see the
Supplementary Material. The spectrum of values of $D(t,n)$, shown in Fig.
3(a), is broadly distributed, $10^{-8}<D(t,n)<0.2$. We use the magnitude of
the velocity dipole field to form a detailed picture of the space-time
structure of events in the gouge. This is illustrated in Figs. 3(b) and (c)
for the large stress drop ($\Delta\mu/\mu\approx 0.2$) near $11000$ msec in
Fig. 2. Figure 3(b) is a zoom of the stress drop; Fig. 3(c) shows the space-
time points at which the velocity dipole field is large. The sharp drop in
Fig. 2 is ragged in close up, Fig. 3(b), and involves a complex set of events
throughout the gouge, Fig. 3(c).

We can also use the large values of $D(t,n)$ advantageously to locate and
examine small, simple events. In Fig. 4(a) there are 3 stress drops of modest
size ($\Delta\mu/\mu\approx 0.002$, about $1\%$ of that in Fig. 3). We look at
the third stress drop in Fig. 4(a). This stress drop occurs rapidly in time,
lasting approximately $3$ msec, and locally in space, involving large $D(t,n)$
at about 20 sensors, Fig. 4(b). To see more details we examine
$V_{X}^{T}(t,n)$ and $V_{X}^{B}(t,n)$ separately in Fig 4(c). We see that
while $V_{X}^{T}(t,n)$ and $V_{X}^{B}(t,n)$ are simultaneous in time they are
structured in space. The B sensors contributing to the velocity dipole field
are further along in $x$ than the associated T sensors. The structural feature
in the gouge, that upon failure delivers the forces that produce the velocity
dipole field in the formation, acts like a strut. It is tethered to the
formation so that it can deliver a shear stress that crosses the gouge and
gives the gouge (an unconsolidated material) a shear modulus. We characterize
the spatial structure of the velocity dipole field associated with small
events like that in Fig. 4(b) with a measure of the separation of the points
of tether, $\Delta n=-n_{X}^{B}-n_{X}^{T}$,

$n_{X}^{*}=\frac{\sum_{n}~{}n~{}d_{X}^{*}(n)}{\sum_{n}~{}d_{X}^{*}(n)},~{}~{}~{}d_{X}^{*}(n)=\int
dt~{}V_{X}^{*}(t,n),$ (3)

where the integral on $t$ and the sum on $n$ are over a domain that surrounds
the small event, $d_{X}^{*}(n)$ is the displacement of sensor $n$, and
$n_{X}^{B}$ ($n_{X}^{T}$) is typically negative (positive). For the spectrum
of $r_{F}=0.24\Delta n/H$ we find the result shown in Fig. 5, where $H$ is the
nominal gouge layer thickness and $0.24$ mm is the spacing between sensors.
While the average value of $r_{F}$ is $\approx 0.5$, and there are both $+$
and $-$ values of $r_{F}$.

We turn from examination of the motion of the formation, as evidenced in the
behavior of the velocity dipole field, to the examination of forces in the
gouge. To do this, for each element in a gouge particle at each time $t$, we
find the maximum principal stress. We take this stress to characterize the
force the element is carrying. When the elements carrying large stress are
illuminated an oriented fabric of forces crossing the gouge is revealed, Fig.
6. Isolated stress chains are an important component of this fabric. These
chains are oriented approximately as suggested by the velocity dipole field,
e.g., Fig. 4(c). To make comparison of the force fabric with the velocity
dipole field we characterize isolated stress chains that cross the gouge with
$r_{G}=\Delta X/H$, where $\Delta X$ is the projection of the chain figure
onto the $x$-axis and $H$ is the width of the gouge. The spectrum of $r_{G}$
for isolated stress chains is shown in Fig. 5.

The velocity dipole field has been used to examine events that are local in
space and time. The two components of the velocity dipole field reveal detail
about the structure in the gouge that produced the event. Often forces are
delivered to the T sensors that are upstream from the forces delivered to the
B sensors. This orientation mirrors the orientation of the fabric of forces in
the gouge, Fig. 6. The most extreme components of this fabric are the isolated
linkages that reach from T to B. From Fig. 5 and Fig. 6 there are more complex
structures in the force fabric. On failure these structures deliver forces to
T and B that are not easily identified with a simple geometry. Comparison of
the spectrum of $r_{G}$ and the spectrum of $r_{F}$ confirms that the geometry
of the stress chains is closely related to the geometry of the velocity dipole
field.

We complete the argument, stress chain $\rightarrow$ acoustic emission, by
examining the motion of points in the formation far from the gouge, i.e., far
from the velocity dipole field. This is done in Fig. S4 of the Supplementary
Material, where it is seen that points in the formation remote from the gouge
move much the same way as the velocity dipole field. The causality is stress
chain $\rightarrow$ velocity dipole field $\rightarrow$ far field signal,
i.e., acoustic emission. The complexity of the numerical model, gouge plus
formation, means that the simulations are quite long. Consequently the
formation does not have a true far field points. The explicit near field
demonstration in Fig. S4 is complemented by use of the representation theorems
in Aki and Richards Aki and Richards (1980) that would have the velocity
dipole field as source for broadcast into the formation.

By the evidence of Fig. 6(b) one might estimate the stress chains to be
separated by approximately the average of their projection onto the x-axis; a
number of order $H$. This suggests that on average about $50$ stress chains
live in the gouge. A single stress chain on failure releases stress that is a
few percent of the stress drop associated with a large slip. This is
consistent with a large slip involving almost all of the extant stress chains,
comparing Fig. 3 and Fig. 4. In addition to assigning stress chains
responsibility for time-sharp events we assign the elasticity evidenced in the
average rise in stress seen in Fig. 2 to their compression.

To conclude, we have employed a FDEM treatment of a numerical model of a
sheared gouge layer. The formation adjacent to the gouge and the interior of
the particles that comprise the gouge have linear elasticity that is resolved
with an FEM treatment. Interactions of gouge particles with one another and
with the formation, involving repulsive-only forces, are resolved with a DEM
treatment. When at fixed normal force the system is sheared at constant drive
velocity the basic stress-time behavior shows intermittent large stress drops
separated by intervals of approximately uniform elastic compression. The
principal actors in the gouge are stress chains. Two sets of sensors at the
gouge-formation interfaces allow determination of the velocity dipole field
which describes the disturbance to the formation caused by the activity in the
gouge. The motion of the formation far from the gouge, taken to be the
acoustic emission, is monitored. We establish the linkage between stress
chains, the velocity dipole field, and the acoustic emission. That is, the
relationship between the activity in the gouge and the far field. This is an
example of a causal scenario that one might hope to encounter and take
advantage of in geophysical systems. Non volcanic tremor, LFEs, etc. are far
fields that may have a role in such a scenario. The great complexity of
geophysical elastic systems makes the demonstration of such a scenario very
difficult. The rapid increase in the quality and quantity of reliable
geophysical monitoring and the infusion of new, sophisticated analysis methods
offer promise that such scenarios are foreseeable. It is worth noting that
acoustic emissions are coming dominantly from the force chain breakage but not
necessarily entirely. Also, the system studied is representative of faults
containing gouge; however, in real faults, acoustic emissions may also arise
from block to block asperity slip and potentially other mechanisms such as
material breakage as the slip front propagates.

###### Acknowledgements.

We acknowledge funding from Institutional Support (LDRD) at Los Alamos
National Laboratory, as well as funding from the US DOE Office of Science,
Geosciences. Technical support and computational resources from the Los Alamos
National Laboratory Institutional Computing Program are highly appreciated.

## References

  * Shelly  _et al._ (2006) D. R. Shelly, G. C. Beroza, S. Ide, and S. Nakamula, Nature 442, 188 (2006).
  * Nadeau and Dolenc (2005) R. M. Nadeau and D. Dolenc, Science 307, 389 (2005).
  * Beroza and Ide (2011) G. C. Beroza and S. Ide, Annu. Rev. Earth Planet. Sci. 39, 271 (2011).
  * Rogers and Dragert (2003) G. Rogers and H. Dragert, Science 300, 1942 (2003).
  * Scholz (2002) C. H. Scholz,  _The mechanics of earthquakes and faulting_ (Cambridge university press, 2002).
  * Rouet-Leduc  _et al._ (2018) B. Rouet-Leduc, C. Hulbert, and P. A. Johnson, Nature Geoscience 12, 75 (2018).
  * Obara and Kato (2016) K. Obara and A. Kato, Science 353, 253 (2016).
  * Kaproth and Marone (2013) B. M. Kaproth and C. Marone, Science 341, 1229 (2013).
  * Johnson  _et al._ (2013) P. A. Johnson, B. Ferdowsi, B. M. Kaproth, M. Scuderi, M. Griffa, J. Carmeliet, R. A. Guyer, P.-Y. Le Bas, D. T. Trugman, and C. Marone, Geophys. Res. Lett. 40, 5627 (2013).
  * Rouet-Leduc  _et al._ (2017) B. Rouet-Leduc, C. Hulbert, N. Lubbers, K. Barros, C. J. Humphreys, and P. A. Johnson, Geophys. Res. Lett. 44, 9276 (2017).
  * Brzinski and Daniels (2018) T. A. Brzinski and K. E. Daniels, Phys. Rev. Lett. 120, 218003 (2018).
  * Peters  _et al._ (2005) J. F. Peters, M. Muthuswamy, J. Wibowo, and A. Tordesillas, Phys. Rev. E 72, 041307 (2005).
  * Daniels and Hayman (2008) K. E. Daniels and N. W. Hayman, J. Geophys. Res. Solid Earth 113 (2008), 10.1029/2008JB005781.
  * Morgan and Boettcher (1999) J. K. Morgan and M. S. Boettcher, J. Geophys. Res. Solid Earth 104, 2703 (1999).
  * Dorostkar  _et al._ (2017) O. Dorostkar, R. A. Guyer, P. A. Johnson, C. Marone, and J. Carmeliet, Geophys. Res. Lett. 44, 6101 (2017).
  * Munjiza (2004) A. A. Munjiza,  _The combined finite-discrete element method_ (John Wiley & Sons, 2004).
  * Geller  _et al._ (2015) D. A. Geller, R. E. Ecke, K. A. Dahmen, and S. Backhaus, Phys. Rev. E 92, 060201 (2015).
  * Munjiza  _et al._ (2011) A. A. Munjiza, E. E. Knight, and E. Rougier,  _Computational mechanics of discontinua_ (John Wiley & Sons, 2011).
  * Gao  _et al._ (2018) K. Gao, B. J. Euser, E. Rougier, R. A. Guyer, Z. Lei, E. E. Knight, J. Carmeliet, and P. A. Johnson, J. Geophys. Res. Solid Earth 123 (2018), 10.1029/2018JB015668.
  * Aki and Richards (1980) K. Aki and P. Richards,  _Quantitative Seismology_ (WH Freeman San Francisco, 1980).

Figure 1: Model system. The formation and the interior of the particles is
described with a FEM scheme. The contacts of gouge particles with one another
and with the formation are described with a DEM scheme. The formation has a
rough surface profile (lower right). Sensors in the formation surface,
adjacent to the gouge, are spaced by $0.24$ mm. The colors of particles
represent the particle number. The elastic properties of gouge and formation
are listed in the Supplementary Material. Figure 2: Primary outputs. The
coefficient of friction (black) as a function of time, left hand scale. The
average velocity of the top (T, blue) and bottom (B, red) sensors as a
function of time, Eq. (1), right hand scale. The green lines, all of the same
slope, show the average elasticity of the gouge. The stress drop near $11000$
msec (marked by red rectangle) is examined in Fig. 3. Figure 3: Spectrum of
velocity dipole strengths, $D$. (a) The spectrum of velocity dipole strengths,
for all space-time points, Eq. (2), is broadly distributed from $10^{-8}$ to
approximately $0.2$. The space-time points $(n,t)$ with large $D$ (in the 3
intervals red, green, blue on the lower right in (a)) are used to identify
events. The stress near $11000$ msec in Fig. 2 is shown in detail in (b). The
velocity dipole strength at $(n,t)$ for this time interval is shown in (c)
where the time axis is the same as in (b) and $n$ is on the vertical axis. The
velocity dipole strength is shown with the color coding from (a). Figure 4:
Details of a small, simple event. (a) Friction coefficient vs time for 3
adjacent small events that have stress drop a few percent of a major event
such as the one in Fig. 3(b). (b) The velocity dipole strengths associated
with the third event in (a), color coded from Fig. 3(a). (c) Velocity dipole
composition: the velocity of the top sensor set (blue) and, displaced to the
right for clarity, and the velocity of the bottom sensor set (red). Figure 5:
Probability distribution functions. Probability distribution function of
$r_{F}$ for the separation of the two components of the velocity dipole field
(small events), blue. Probability distribution function of $r_{G}$ for
isolated stress chains, red. Figure 6: Stress chains. (a) The spectrum of
maximum principal stress of all elements in the gouge. (b) Elements in the FEM
description of the particles, with maximum principal stress greater than $0.4$
MPa, are shown in black according to their centroid positions. Two chains of
elements at the right, 1 and 2, are judged to be isolated as is the chain 3 on
the left. These chains are among the chains used to form the pdf of $r_{G}$ in
Fig. 5.