# De novo reconstruction of satellite repeat units from sequence data

Yujie Zhang1 Justin Chu2,3 Haoyu Cheng2,3 Heng Li2,3 1Harvard School of Public
Health, 677 Huntington Avenue, Boston, MA 02115, USA, 2Department of Data
Science, Dana-Farber Cancer Institute, 450 Brookline Ave, Boston, MA 02215,
USA, 3Department of Biomedical Informatics, Harvard Medical School, 10
Shattuck St, Boston, MA 02115, USA

(2022; 2022)

###### Abstract

Satellite DNA are long tandemly repeating sequences in a genome and may be
organized as high-order repeats (HORs). They are enriched in centromeres and
are challenging to assemble. Existing algorithms for identifying satellite
repeats either require the complete assembly of satellites or only work for
simple repeat structures without HORs. Here we describe Satellite Repeat
Finder (SRF), a new algorithm for reconstructing satellite repeat units and
HORs from accurate reads or assemblies without prior knowledge on repeat
structures. Applying SRF to real sequence data, we showed that SRF could
reconstruct known satellites in human and well-studied model organisms. We
also found satellite repeats are pervasive in various other species,
accounting for up to 12% of their genome contents but are often
underrepresented in assemblies. With the rapid progress on genome sequencing,
SRF will help the annotation of new genomes and the study of satellite DNA
evolution even if such repeats are not fully assembled.

##  1 Introduction

Satellite DNA (SatDNA) are long tandemly repeating sequences that look like
“BBBBBB$\cdots$”, where each symbol “B” represents a repeat unit, also known
as a _monomer_. A monomer “B” could range from a few basepairs (bp) to
thousands of bp in length and an entire SatDNA could span megabases in large
genomes. Several percent of the human genome, or a couple of hundred megabases
in total, is composed of SatDNA (Altemose et al.,, 2022). Monomers in a SatDNA
array are similar in sequence but often not identical due to random mutations.

In some species, SatDNA may be organized as high-order repeats (HORs; Miga,
2019). For example, the centromere of human chromosome 2 has a pattern like
“ABCDABCDABCD$\cdots$”. Letters “A”–“D” correspond to four diverged alpha
repeat monomers of $\sim$171bp each, respectively, and the “ABCD” unit is
repeated many times in the centromere with all copies being similar to each
other. Researchers who study centromere repeats usually say “ABCD” is a 4-mer
HOR unit. Because we will often mention short nucleotide sequences in this
article, we will call “ABCD” in this example as 4-monomer HOR unit to avoid
confusion.

SatDNA is often not assembled in long contigs due to its repetitiveness. We
would have to reconstruct SatDNA from raw sequence reads in this case (Lower
et al.,, 2018). Wei et al (Wei et al.,, 2014) developed k-Seek to study SatDNA
consisting of 2–10bp repeat units. Melters et al (Melters et al.,, 2013)
applied Tandem Repeat Finder (TRF; Benson, 1999) to Sanger reads and
fragmented short-read contigs to find the most common monomer in each species.
TAREAN (Novák et al.,, 2017) does all-vs-all comparison between short reads,
clusters the reads and then identifies circular structures from the cluster
graphs. These methods can reconstruct unknown monomers but they are unable to
reveal HOR structures. On the contrary, Alpha-CENTAURI (Sevim et al.,, 2016)
reconstructs HORs from long reads but it requires known monomer sequences.

With improved sequencing technologies, it is now possible to assemble through
human centromeres (Nurk et al.,, 2022). More recent methods, including
NTRprism (Altemose et al.,, 2022), HORmon (Kunyavskaya et al.,, 2022) and
HiCAT (Gao et al.,, 2022), can identify detailed chromosome-specific HOR
patterns from complete SatDNA sequences in human. These methods demand high-
quality assembly and rich prior knowledge on SatDNA in the studied species.
However, the finished human genome, CHM13, was derived from a near homozygous
molar cell line that is easier to assemble. For a normal diploid human
individual, we could only assemble through a fraction of SatDNA even with the
best possible data and algorithm (Rautiainen et al.,, 2023). The complete
assembly of other species is even rarer. This has limited the application of
such assembly-based SatDNA reconstruction algorithms.

In this article, we will describe a new algorithm, Satellite Repeat Finder
(SRF), for assembling SatDNA repeat units. SRF overcomes the limitation of
previous methods. It is applicable to both accurate reads and high-quality
assembly and is able to automatically reconstruct HORs with no prior knowledge
on monomer sequences.

##  2 Results

###  2.1 The SRF algorithm

In a SatDNA array “BBBBBB$\cdots$”, suppose every monomer “B” is identical to
each other. Under a long enough $k$, the $k$-mer de Bruijn graph of the SatDNA
array will be a single cycle. When there are basepair differences between
monomers, the de Bruijn graph will not be a simple cycle. If there are many
different monomers, the de Bruijn graph can become very complex and cannot be
resolved with classical graph cleaning algorithms (Zerbino and Birney,, 2008).

Our intuition is that if there are many copies of the monomer, we may still be
able to find a cycle composed of highly abundant $k$-mers in the de Bruijn
graph. We can start with the most abundant $k$-mer and at each bifurcation in
the graph, we greedily choose the $k$-mer of the highest occurrence. We repeat
this process until we go back the starting $k$-mer, which will reconstruct a
repeat unit, or come to a deadend, which will be discarded. Algorithm 1
provides more details. Here, $t\to s$ indicates $k$-mer $t$ and $s$ are
adjacent in the de Bruijn graph. For simplicity, this algorithm traverses a
unidirected de Bruijn graph. In SRF, we implemented a bidirected de Bruijn
graph such that we will not find a repeat unit on both strands.

SRF works with Illumina short reads, PacBio HiFi long reads and high-quality
assembly contigs and can identify HORs (Table 2.1). When assembling satellite
repeats from PacBio HiFi reads in this article, we counted 151-mers with KMC
(Kokot et al.,, 2017) and collected 151-mers occurring $\geq$10 times over the
average read coverage. K-mer counting may take a few tens of minutes for a
high-coverage human dataset and is the performance bottleneck. SRF only takes
seconds to reconstruct all repeat units after k-mer counting.

Features of user-facing tools for SatDNA reconstruction  Tool Reads Contigs De
novo HORs Alpha-CENTAURI Yes No No Yes HiCAT No Yes No Yes HORmon No Yes No
Yes NTRprism No Yes Yes Yes SRF (this work) Yes Yes Yes Yes TAREAN Yes No Yes
No “Reads”: whether the tool works with unassembled reads. “Contigs”: whether
the tool works with high-quality contigs. “De novo”: whether the tool works
without known monomer sequence as input. “HORs”: whether the tool can identify
high-order repeats (HORs). References to these tools can be found in the
Introduction section.

###  2.2 Estimating satellite abundance

The SRF algorithm does not provide a good estimate of repeat abundance. We
mapped all input sequences against reconstructed repeat units to measure the
total length of each repeat unit. For human CHM13 data, we observed many
diverged hits between HORs and scattered monomers in pericentromeric regions.
We thus filtered hits of identity below 90% to get more accurate HOR length
estimates. The effect of the identity is determined by the repeat structure in
a species. For example, switching off the filter would increase the total
abundance estimate by 40% for human but only by 4% for _A.thaliana_. We still
applied this filter to all datasets even though this may lead to
underestimates for some species.

SRF may reconstruct repeat units similar in sequence. The similar repeat units
may be mapped the same genomic locus. To remove redundancy, we only select the
hit of the highest identity among hits overlapping on an input sequence. With
this procedure, we map each base on an input sequence to at most one repeat
unit.

Occasionally a small number of long terminal repeats (LTRs) may occur tandemly
in a few region. SRF may identify such LTRs even though they do not form long
tandem arrays. When estimating abundance, we additionally filter out repeat
with $<$2 tandem copies in the middle of a sequence or with $<$1.5 tandem
copies when the repeat-to-read alignment reaches the end of a read. This
filter is reliable when we apply SRF to assemblies but may miss long repeat
units when applied to reads. We again opted for conservative estimates.

###  2.3 A brief introduction to human satellites

In the human genome, the most abundant satellite family is alpha satellites
with most of them present in long alpha HORs ($\alpha$HORs). The consensus of
the minimal alpha repeat unit is 171bp in length. The active centromeric
regions that centromeric proteins bind to are primarily composed of
$\alpha$HORs (Altemose et al.,, 2022). Conversely, though, not all
$\alpha$HORs are present in the active regions. These inactive $\alpha$HORs
tend to be shorter than active ones. Alpha repeat monomers are also present in
pericentrometic regions without clear HOR structures. In addition to alpha
repeat, the human genome is also enriched with three types of human satellites
(HSat1–3), contributing to a few percent of human genome (Altemose,, 2022).
Almost all these satellites are located around centromeres or on the long arm
of the Y chromosome.

In the human reference genome GRCh38 (Schneider et al.,, 2017), $\alpha$HORs
were computationally generated from a Hidden Markov Model (Miga et al.,,
2014); HSats are underrepresented. At present, only the T2T-CHM13 assembly
(Nurk et al.,, 2022) provides a complete representation of all satellite
arrays.

Input: $k$-mer set $\mathcal{S}$ and count function
$\phi:\mathcal{S}\to\mathbb{N}$

Output: SatDNA repeat units

while  _$\mathcal{S}\not=\emptyset$_ do

$t_{0}\leftarrow\operatorname*{argmax}_{s\in\mathcal{S}}\phi(s)$

$\triangleright$ Most abundant $k$-mer in $\mathcal{S}$

$\mathcal{H}\leftarrow\\{t_{0}\\}$

while  _$\mathcal{H}\not=\emptyset$_ do

$t\leftarrow\operatorname*{argmax}_{s\in\mathcal{H}}\phi(s)$

$\triangleright$ Most abundant $k$-mer in $\mathcal{H}$

$\mathcal{H}\leftarrow\mathcal{H}\setminus\\{t\\}$

$\triangleright$ Remove $t$ from $\mathcal{H}$

if  _$t=t_{0}$ and $t_{0}\not\in\mathcal{S}$_ then $\triangleright$ Back to
the starting $k$-mer

Trace back from $t$ and report a repeat unit

break

for  _$s\in\mathcal{S}$ and $t\to s$_ do $\triangleright$ Traverse $t$’s
neighbors

$\mathcal{S}\leftarrow\mathcal{S}\setminus\\{s\\}$

$\mathcal{H}\leftarrow\mathcal{H}\cup\\{s\\}$

$\triangleright$ Move $s$ from $\mathcal{S}$ to $\mathcal{H}$

${\rm parent}(s)=t$

$\triangleright$ Keep parent for traceback

$\mathcal{S}\leftarrow\mathcal{S}\setminus\\{t_{0}\\}$

$\triangleright$ Remove $t_{0}$ from $\mathcal{S}$

Algorithm 1 Assemble SatDNA repeat units

Human chromosome-specific high-order alpha repeats ($\alpha$HORs)  chr HORmon
a HiCAT b SRF ($k$=171) SRF/171 SRF/171 SRF/101 SRF ($k$=171) centromere
centromere chromosome c assembly d HiFi reads d Illumina d HPRC assembly e 1 6
2 6 (4.2); 11 (0.5) 6 (2.0) 6 (3.5) 2 (2.2) 6 [89] 2 4 4 4 (2.3) 4 (2.3) 4
(2.2) 4 (2.2) 4 [94] 3 17 17 17 (1.4) 17 (1.4) 17 (1.4) 17 (1.4) 17 [94] 4 19
19 19 (3.5) 19 (3.5) 19 (2.9) 19 (3.4) 19 [94] 5 6 12 8 (2.5) 8 (1.8) 4 (1.9)
missing 8 [43]; 4 [37] 6 18 18 18 (2.0) 18 (2.0) 18 (2.0) 18 (2.0) 18 [93] 7 6
6 6 (3.3) 6 (3.2) 6 (3.2) 6 (3.2) 6 [92]; 12 [2] 8 11 15 7 (1.1) 7 (1.1) 7
(1.0) 11 (1.0) 7 [61]; 8 [33] 9 7 11 4 (1.8) 4 (1.4) 11 (2.0) 4 (1.7) 4 [77];
11 [17] 10 8 6 8 (2.1) 8 (2.1) 8 (1.7) 8 (2.1) 6 [66]; 8 [28] 11 5 5 5 (3.4) 5
(3.3) 5 (3.4) 5 (3.4) 5 [94] 12 8 8 8 (2.6) 8 (2.6) 8 (2.6) 8 (2.6) 8 [94] 13
11 7 4 (0.4) 4 (0.4) 7 (1.5) 7 (1.5) 4 [55]; 11 [23]; 7 [16] 14 8 8 8 (2.6)
missing missing missing missing 15 11 15 11 (0.8); 20 (0.5) 11 (0.8) 11 (0.8)
11 (0.8) 11 [94] 16 10 10 10 (2.0) 10 (1.9) 10 (1.9) missing 10 [94] 17 16 14
16 (3.3) 16 (3.3) 16 (3.5) 16 (3.5) 16 [56]; 13 [38] 18 12 12 8 (3.6) 8 (3.8)
12 (4.9) missing 12 [66]; 8 [19] 19 2 2 4 (0.4); 2 (0.4) missing 13 (0.5)
missing 13 [29]; 32 [4] 20 16 16 16 (2.1) 16 (2.1) 16 (2.1) 8 (0.5) 16 [94] 21
11 11 11 (0.3) missing missing missing missing 22 8 8 8 (2.9); 20 (0.5) 8
(2.8) 8 (2.6) 8 (2.9) 8 [94] X 12 12 12 (3.1) 12 (3.1) 12 (3.1) 12 (3.1) 12
[76] Y 34 No Y 34 (0.3) 34 (0.3) No Y No Y 34 [18] a $\alpha$HOR lengths in
the monomer unit in the CHM13 v2.0 genome, retrieved from Kunyavskaya et al.,
(2022). b length of “top 1” $\alpha$HOR from each chromosome retrieved from
Gao et al., (2022). Both HORmon and HiCAT were applicable to extracted
centromeric sequences only.c SRF applied to each CHM13 chromosome separately.
In a format “$m$ ($L$)”, $m$ denotes the length of an HOR in the monomer unit
and $L$ is its span on the CHM13 assembly in megabases. d SRF applied to CHM13
assembly, PacBio High-Fidelity (HiFi) reads and Illumina short reads,
respectively. $k$=101 used for Illumina reads. CHM13 reads do not contain
chrY. e SRF applied to 94 phased haploid assemblies produced by the Human
Pangenome Reference Consortium (HPRC). In a format “$m$ [$n$]”, $m$ is the
monomer length and $n$ is the number of samples with the HOR according to
manual inspection.

###  2.4 Annotating satellite repeats in human T2T-CHM13

We first ran SRF on each T2T-CHM13 chromosome separately and compared the
results to existing annotations by HORmon (Kunyavskaya et al.,, 2022) and
HiCAT (Gao et al.,, 2022). HORmon reports the same $\alpha$HOR lengths as
Altemose et al., (2022). In Table 2.3, column “SRF ($k$=171) chromosome” shows
the lengths of HORs identified by running SRF on individual chromosomes. SRF
reported the same $\alpha$HOR lengths as HORmon except for chromosome 5, 8, 9,
13 and 18. The chromosome 8 of T2T-CHM13 has been well studied by Logsdon et
al., (2021). Although the 11-monomer is the most abundant, it is interleaved
with 4-, 7- and 8-monomers that are derived from the 11-monomer. The 7-monomer
is the second most abundant array and forms the longest $\alpha$HOR array in
the middle of the centromere. The greedy SRF algorithm chooses the 7-monomer
over the 11-monomer possibly because the 7-monomer has a more conservative
consensus. The SRF-HORmon differences in other chromosomes may have a similar
cause.

Unlike HORmon and HiCAT which require users to provide the monomer sequence
and prepare centromeric sequences, SRF was directly applied to whole
chromosome sequences with no prior knowledge. In addition to active
$\alpha$HORs, SRF identified shorter $\alpha$HOR arrays outside the active
regions. It also found many long non-alpha satellite arrays including a repeat
unit of 1,814 bp on chr15, of 6,112 bp on chr16, of 3,569 bp on chrY and of
2,420 bp on chrY as well. These span over one megabase and have been reported
previously (Altemose,, 2022).

SRF further found a satellite array on the long arm of chromosome 1 between
coordinate 227,746,662 and 228,024,151. The repeat unit is 2,240 bp in length,
composed of an AluY repeat, a 5S-RNA and di-nucleotide repeats. This is the
only non-centromeric array in T2T-CHM13 longer than 100kb.

The SRF inference on the whole T2T-CHM13 genome (column “SRF/171 assembly” in
Table 2.3) is close to the inference on individual chromsomes. SRF missed the
$\alpha$HOR array on chr14 and chr21 because chr22 and chr13, respectively,
have very similar arrays which are merged during the whole-genome inference.

SRF works on sequence reads which HORmon, HiCAT and NTRprism are not
applicable to. On PacBio High-Fidelity (HiFi) reads, SRF reconstructed
$\alpha$HORs similar to the whole-genome reconstruction (column “SRF/171 HiFi
reads”). It can also identify the majority of $\alpha$HORs from Illumina short
reads (column “SRF/101 Illumina” in Table 2.3), though the use of shorter
101-mer reduces the sensitivity to some arrays.

From Table 2.3 we can see that some $\alpha$HOR arrays, such as those on chr3
and chr11, can be consistently reconstructed by various tools on different
types of input data. However, some other arrays, such as those on chr8 and
chr19, are intrinsically harder to reconstruct. These are probably because
monomers in a HOR may be connected in different ways, as is shown by
Kunyavskaya et al., (2022).

###  2.5 Satellite repeats in multiple human assemblies

We applied SRF independently to each phased human assembly produced by the
Human Pangenome Reference Consortium (HPRC). We identified $\alpha$HORs with
dna-brnn (Li,, 2019), aligned them to the T2T-CHM13 genome and HORmon
consensus and manually assigned the $\alpha$HORs to chromosomes based on the
similarity to existing annotations. The last column of Table 2.3 shows the
$\alpha$HOR lengths and their frequencies. We consider two $\alpha$HORs are
different if they have different lengths and the shorter $\alpha$HOR cannot be
aligned into the longer one at $<$2% sequence divergence.

There are 47 diploid samples and 94 haploid assemblies. We could find
$\alpha$HORs assigned to individual chromosomes in most cases. We sometimes
see $\alpha$HORs of different lengths assigned to the same chromosome but
their sequence divergences are small. This again could be caused by the
different ways HOR monomers are connected in individual samples (Logsdon et
al.,, 2021; Kunyavskaya et al.,, 2022).

We also applied SRF to the pool of all haploid assemblies. This procedure may
miss infrequent satellite arrays but it helps to simplify the study of shared
arrays. In addition to active $\alpha$HORs, we identified supposedly inactive
$\alpha$HORs that are at $>$10% divergence from existing annotated
$\alpha$HORs. Notably, there is a 20-monomer $\alpha$HOR that is mapped to the
chromosome 15 of T2T-CHM13 and spans several hundred kilobases in most
samples. There are other examples like this. SRF also found long HSat arrays
and non-HOR satellites, including those found in CHM13, which are easy to
identify as they do not have internal structures.

###  2.6 Satellite repeats in _Arabidopsis thaliana_

We obtained three HiFi datasets (Table 2.6) and downsampled them to about
40-fold coverage each. Col-0N and Ey15-2R were sequenced from a pool of
multiple samples and Col-0R from a single sample.

_A. thaliana_ PacBio HiFi datasets  Sample Strain SRA Accession Source Col-0N
Col-0 ERR6210723 Naish et al., (2021) Col-0R Col-0 ERR8666127 Rabanal et al.,
(2022) Ey15-2R Eyach15-2 ERR8666125 Rabanal et al., (2022)

The _A. thaliana_ centromeres are composed of CEN180 satellites which are
known to have high-order organizations (Naish et al.,, 2021). When we pooled
the two Col-0 datasets and applied SRF, we reconstructed one 1-monomer, two
2-monomers, one 3-monomer and one 6-monomer HORs. However, this result is
unstable. If we ran SRF on each Col-0 dataset separately or on the full
coverage, we would get different HORs. We thus ignored the HOR structures and
focused on the abundance of CEN180 only.

Figure 1: Normalized counts of 179-mers in three _A. thaliana_ read datasets.
Raw 179-mer counts in reads are normalized by coverage. A 179-mer is selected
in the plot if it matches the CEN180 satellite and if its normalized count is
at least 50 in one of the datasets. (a) Counts between two different samples
from the same strain. (b) Counts between two different strains.

SRF estimated that 5.2% of read bases in Col-0N and 7.5% in Col-0R are
composed of CEN180 satellites. Col-0R has 44% more CEN180 satellites than
Col-0N. To check whether this large difference is caused by an artifact, we
inferred the relative CEN180 abundance with a different approach as follows.
We collected 179-mers matching CEN180 and occurring $\geq$400 times in the
pooled Col-0 dataset. We counted the total numbers of these 179-mers in Col-0N
and Col-0R separately and normalized the counts by read bases. We found Col-0R
has 41% more CEN180 179-mers than Col-0N, broadly in line with the earlier
analysis. Centromere contents may differ greatly even between two samples from
the same strain.

SRF estimated that 11.5% of Ey15-2R is composed of CEN180 satellites, higher
than both Col-0R and Col-0N. Furthermore, while Col-0R and Col-0N share
similar high-occurrence 179-mers that match CEN180 (Fig. 1a), Col-0 and Ey15-2
share few common 179-mers (Fig. 1b). The centromere sequences between strains
are distinct both in content and in length.

In addition to the CEN180 satellite, SRF also reconstructed a 10,067 bp rDNA
unit from the two Col-0 datasets. It has 3.1% abundance in Col-0R and 1.7% in
Col-0N. If we assume the _A. thaliana_ genome is 132Mb in length according to
the Col-0N assembly, Col-0R has $\sim$400 copies of this rDNA unit while
Col-0N has $\sim$220 copies. The Col-0N assembly (Naish et al.,, 2021) only
has seven copies, located towards the telomeric ends of chr2 or chr4 short
arms. SRF did not reconstruct an rDNA unit from Ey15-2R. We mapped the Col-0
rDNA unit to Ey15-2R reads and estimated that Ey15-2R has $\sim$200 copies.

###  2.7 Satellite repeats in other model organisms

We applied SRF to the HiFi reads of three model organisms (Hon et al.,, 2020):
the reference C57BL/6J strain of _Mus musculus_ (mouse; AC:SRR11606870), the
F1 generation of the reference ISO1 strain and the A4 strain of _Drosophila
melanogaster_ (AC:SRR10238607), and the B73 strain of _Zea mays_ (maize;
AC:SRR11606869).

In mouse, SRF identified two satellite units. The second most abundant repeat
is the 234 bp major satellite around centromeres (Arora et al.,, 2021; Thakur
et al.,, 2021). The first is 1,199 bp in length, composed of 10 copies of the
120 bp minor satellite unit. This confirms the high-order organization of
minor satellites observed by Pertile et al., (2009). The full-length hits of
this repeat in the mouse reference genome mostly come from the sex chromosomes
and are all below 75% in identity. Nonetheless, this repeat is abundant in
reads with the majority of alignments at 95% identity or higher. To further
investigate this repeat, we assembled the HiFi reads with hifiasm (Cheng et
al.,, 2021). We can find long tandem arrays of this repeat on multiple
contigs, all shorter than 1.1Mb. Hifiasm keeps the repeat content but is
unable to assemble this satellite.

In _Drosophila_ , the most abundant satellite SRF identified is a 358 bp
repeat unit hitting 0.90% of read bases. It belongs to the 1.688 family (Khost
et al.,, 2017). The abundance of the 358 bp repeat is lower in the BDGP6
reference genome, at 0.26% only. SRF assembled the 240 bp Intergenic spacer
(IGS; Shatskikh et al., 2020) into two sequences, at 240 bp and 239 bp,
respectively. The edit distance between the two IGS sequences is 5. They hit
to 0.43% of read bases in total but are depleted in the reference at $<$0.01%
only. SRF also found other known satellite repeats such as (AAGAC)n, (AACAC)n,
(AATAG)n, (GGTCCCGTACT)n and (AATAACATAG)n (Shatskikh et al.,, 2020; Thakur et
al.,, 2021). There are more copies of these repeat units but because they are
short, they contribute less to the genome in comparison to the 1.688 and IGS
satellites.

SRF reconstructed a 5,045 bp repeat unit at 0.34% abundance in reads and 0.08%
in the reference genome. It harbors histone genes and is located in a small
region on chromosome 2L. To investigate further, we assembled the HiFi reads
using the hifiasm trio-binning mode with ISO1 and A4 short reads from
SRR6702604, SRR457665, SRR457666 and SRR457707. When aligning the ISO1
haplotype assembly to the reference genome, we see a clean 242kb insertion
entirely composed of the 5,045 bp histone repeat. The insertion has 48 tandem
copies at $>$99% identity between the copies. The BDGP6 reference genome might
have misassembled this region.

In maize, SRF reconstructed a 741 bp repeat unit at 0.25% abundance. It
matches the SAT1_ZM record in RepBase. This SRF unit includes four copies of a
180 bp knob-associated repeat (Ananiev et al., 1998b, ). In the NAM-5.0
reference genome or the hifiasm assembly, this repeat tends to be present in
short contigs and towards ends of long contigs. It is not assembled well. SRF
also identified many potential repeat units at $<$0.09% abundance in reads.
Nonetheless, none of them form long tandem arrays. Meanwhile, under the
151-mer setting, SRF failed to identify the 156 bp CentC repeat (Ananiev et
al., 1998a, ). SRF could find this repeat if we counted 101-mers. Only 0.045%
of read bases were mapped to CentC. Low-abundance SatDNA is harder to assemble
correctly.

###  2.8 Comparison to TAREAN

TAREAN (Novák et al.,, 2017) can identify novel satellite repeats from
sequence reads. Its developers recommend to use reads at up to 0.5-fold
coverage to avoid redundancy between reads sequenced from the same loci. We
ran TAREAN on simulated short reads at 0.2-fold from the Dropsophila HiFi
dataset described above, without introducing additional sequencing errors.
TAREAN found six high-confidence satellite repeats, including the 1.688 family
and the histone cluster, (GGTCCCGTACT)n and (AATAACATAG)n. The other two
TAREAN repeats also hit to SRF contigs. SRF assembled eight more SatDNA repeat
units at $>$0.05% abundance. Manually inspecting the alignment of SRF contigs
to raw HiFi reads, we observed tandem pattern for all of them, suggesting they
were real SatDNA.

To evaluate whether TAREAN can reconstruct HORs, we ran TAREAN on 0.2-fold
CHM13 reads randomly sampled from SRR2088062. TAREAN took 5 hours and found
four high-confidence satellite repeats, including a 2-monomer alpha repeat at
1.0% abundance, a HSat2 repeat at 0.9%, a SAR satellite and a beta satellite.
TAREAN did not identify other HORs.

###  2.9 Satellite repeats in other species

HiFi datasets for non-model organisms  Species Common name Source _A.
bisporus_ cultivated mushroom PRJEB52214 _A. ruthenus_ sterlet PRJEB19273 _C.
elaphus_ red deer Pemberton et al., (2021) _C. lupus_ grey wolf Sinding et
al., (2021) _D. lineata_ orange-striped anemone Wood et al., (2022) _F.
$\times$ ananassa_ royal royce strawberry Hon et al., (2020) _H. helix_ ivy
PRJEB47300 _L. sulphureus_ chicken mushroom Wright et al., (2022) _M.
domestica_ apple Könyves et al., (2022) _M. meles_ Eurasian badger Newman et
al., (2022) _M. sylvestris_ crab apple Ruhsam et al., (2022) _M. glacialis_
spiny starfish Lawniczak et al., (2021) _O. orca_ killer whale Foote et al.,
(2022) _O. sativa_ rice SRR10238608 _P. pellucida_ blue-rayed limpet Lawniczak
et al., (2022) _R. muscosa_ yellow-legged frog Hon et al., (2020) _V.
atalanta_ red admiral butterfly Lohse et al., (2021)

Figure 2: Abundance of satellite DNA in 21 species.

We randomly selected 14 species from the Darwin Tree of Life project and
collected two species from Hon et al., (2020) (Table 2.9). We assembled SatDNA
in these and several other species described in earlier sections. SRF may
reconstruct mitochondria or chloroplast from sequence reads. We manually
removed them based on NCBI BLAST against the nt database. We then estimated
the abundance of SatDNA in each of these species (Fig. 2).

Red deer (_C. elaphus_) has the highest abundance at 11.9%. A single 796 bp
repeat unit accounts for 10.3% of satellite DNA. Killer whale (_O. orca_) in
the same order is also enriched with satellite DNA. Yellow-legged frog (_R.
muscosa_) is next to killer whale. SRF reconstructed many distant variants of
a 131 bp repeat unit. On the other extreme, apples (_M. domestica_ and _M.
sylvestris_) barely have satellite repeats partly because they have
transposon-rich centromeres (Zhang et al.,, 2019). Plants and fungi are
generally depleted of satellites.

It is worth noting that our abundance estimate may be an underestimate due to
the additional filters we used. For example, chicken mushroom (_L.
sulphureus_) had a repeat unit of 9,659 bp at 0.9% abundance. As we discarded
alignments shorter than 1.5 times 9,659 bp, we filtered out many HiFi reads
shorter than this threshold even if entire reads were aligned to the repeat.
The abundance estimate would be doubled without this filter. Such long repeat
units are infrequent in the species we studied.

To investigate what satellites are organized as HOR, we ran TRF (Benson,,
1999) on SRF-assembled repeat motif. A repeat motif is considered to have a
high-order structure if TRF identifies a tandem repeat repeating at least
three times and covering 90% of the motif. With this criterion, 98.5% of human
satellites are HORs with a variety of number of monomers. 8.6% of satellites
in Eurasian badger (_M. meles_) are HORs of a 138 bp monomer, contributing to
0.14% of the genome. The other species in our survey either do not have
multiple HORs composed of similar monomers or only have HORs at $<$0.1%
abundance. Consistent with our observation, Melters et al., (2013) rarely
identified HORs consisting of $\geq$3 monomers. The authors attributed this to
the limited Sanger read length. Based on longer reads and a different
algorithm, our result suggests that most species do not exhibit rich HOR
structures.

##  3 Discussions

SRF is a de novo assembler for reconstructing SatDNA repeat units and can
identify most known HORs and SatDNA in well-studied species without prior
knowledge on monomer sequences or repeat structures. It is the only de novo
algorithm for reconstructing HORs from sequence reads as well as high-quality
assemblies. SRF only depends on a third-party k-mer counter. It is easy to run
and fast to execute.

SRF uses a greedy algorithm to assemble SatDNA repeat units. When two repeat
units share long similar sequences, the one of lower abundance and higher
diversity may be missed. We plan to improve the current algorithm by reporting
multiple overlapping cycles. This may be able to find a more complete
collection of HORs in the human genome.

Meanwhile, although SRF can reconstruct known HORs in human, it may report
incidental HORs in species, such as mouse and _A. thaliana_ , that only have
weak high-order patterns. We need to run TRF (Benson,, 1999) on SRF contigs to
obtain minimal repeat units. SRF may also assemble the same class of repeat
into multiple similar but not identical copies. We can align assembled repeat
units to identify such redundancy. Manual curation is recommended for a deeper
insight into the SatDNA structure of a new species.

Estimating the abundance of SatDNA is challenging. Sometimes ancient SatDNA
repeats may be too diverged from the assembled repeat consensus to be aligned
confidently. In human, whether to count scattered monomers in pericentromeric
regions as long SatDNA arrays would affect the estimate as well. In addition,
occasionally SatDNA units can be $>$5 kb in length. We may not observe clear
tandem patterns in $\sim$10 kb HiFi reads, which would lead to underestimate.
We do not have an automated algorithm to provide accurate abundance estimate
in corner cases.

SatDNA is pervasive in many species. It is however often underrepresented in
current reference genomes such as the human GRCh38 genome and the _Drosophila_
BDGP6 genome. Even with improved sequencing technologies and assembly
algorithms, the assembly of SatDNA is often fragmented. With thousands of
species sequenced recently (Challis et al.,, 2020; Rhie et al.,, 2021) and
more to come in future, SRF may become an important tool to identity and
annotate SatDNA in these species. It may also supplement RepeatModeler (Flynn
et al.,, 2020) to provide a more comprehensive repeat library for masking
SatDNA in assembled genomes.

##  4 Methods

###  4.1 Running SRF for human assemblies

We counted 171-mers occurring 20 times or more with KMC, using command line
kmc -fm -k171 -ci20 -cs100000 and extracted the 171-mer counts with kmc_dump.
SRF is directly applied to the output of kmc_dump output in the default
setting.

###  4.2 Running SRF on sequence reads

We estimated the approximate read depth by dividing the total number of read
bases by the number of bases in the reference genome or the corresponding read
assembly. We counted 151-mers with kmc -fq -k151 -ciXX -cs1000000, where XX is
10 times the average read depth of each sample.

###  4.3 Estimating the abundance of SatDNA

We aligned reconstructed repeat units to HiFi reads or contigs with minimap2
(Li,, 2018), using command line minimap2 -c -N1000000 -f1000 -r100,100
<(srfutils.js enlong srf.fa), where srfutils.js is a companion script along
with the SRF tool. Option -N1000000 asks minimap2 to report up to a million
hits per query sequence; -f1000 considers high-occurrence seeds; -r100,100
enables a small bandwidth of 100 bp during alignment.

After the alignment, we used srfutils.js paf2bed to filter poor alignments and
to merge adjacent alignments, and then used srfutils.js bed2abun to calculate
the abundance of each repeat unit.

###  4.4 Running TAREAN

For human CHM13, we used real reads short reads reads. We ran TAREAN with
singularity exec --bind ${PWD}:/data/ shub://repeatexplorer/repex_tarean
seqclust -p -c 32 -r 50000000. For _Drosorphila_ , we simulated 125 bp paired-
end reads from HiFi reads with dwgsim -N 146000 -1 125 -2 125 -y0 -e0 -E0 -r0
-F0 -R0. This command line did not add additional sequencing errors; the
simulated reads only carried real sequencing errors on the original HiFi
reads.

##  5 Data access

The SRF implementation and associated analysis scripts are provided at
https://github.com/lh3/srf. A modified TRF with an alternative command-line
interface is available at https://github.com/lh3/TRF-mod. Assembled repeat
units and their abundance estimates can be found at
https://zenodo.org/record/7814465.

##  6 Competing interest statement

H.L. is a consualtant for Integrated DNA Technologies, Inc.

##  7 Acknowledgements

This work is supported by US National Human Genome Research Institute (NHGRI)
grant R01HG010040 and U01HG010961 to H.L.

## References

  * Altemose, (2022) Altemose, N. (2022).  A classical revival: Human satellite dnas enter the genomics era.  Semin Cell Dev Biol, 128:2–14.
  * Altemose et al., (2022) Altemose, N., Logsdon, G. A., Bzikadze, A. V., Sidhwani, P., Langley, S. A., Caldas, G. V., Hoyt, S. J., Uralsky, L., Ryabov, F. D., Shew, C. J., et al. (2022).  Complete genomic and epigenetic maps of human centromeres.  Science, 376:eabl4178.
  * (3) Ananiev, E. V., Phillips, R. L., and Rines, H. W. (1998a).  Chromosome-specific molecular organization of maize (zea mays l.) centromeric regions.  Proc Natl Acad Sci U S A, 95:13073–8.
  * (4) Ananiev, E. V., Phillips, R. L., and Rines, H. W. (1998b).  Complex structure of knob dna on maize chromosome 9. retrotransposon invasion into heterochromatin.  Genetics, 149(4):2025–37.
  * Arora et al., (2021) Arora, U. P., Charlebois, C., Lawal, R. A., and Dumont, B. L. (2021).  Population and subspecies diversity at mouse centromere satellites.  BMC Genomics, 22:279.
  * Benson, (1999) Benson, G. (1999).  Tandem repeats finder: a program to analyze DNA sequences.  Nucleic Acids Res, 27:573–80.
  * Challis et al., (2020) Challis, R., Richards, E., Rajan, J., Cochrane, G., and Blaxter, M. (2020).  BlobToolKit - interactive quality assessment of genome assemblies.  G3 (Bethesda), 10(4):1361–1374.
  * Cheng et al., (2021) Cheng, H., Concepcion, G. T., Feng, X., Zhang, H., and Li, H. (2021).  Haplotype-resolved de novo assembly using phased assembly graphs with hifiasm.  Nat Methods, 18:170–175.
  * Flynn et al., (2020) Flynn, J. M., Hubley, R., Goubert, C., Rosen, J., Clark, A. G., Feschotte, C., and Smit, A. F. (2020).  Repeatmodeler2 for automated genomic discovery of transposable element families.  Proc Natl Acad Sci U S A, 117(17):9451–9457.
  * Foote et al., (2022) Foote, A., Bunskoek, P., Wellcome Sanger Institute Tree of Life programme, Wellcome Sanger Institute Scientific Operations: DNA Pipelines collective, Tree of Life Core Informatics collective, and Darwin Tree of Life Consortium (2022).  The genome sequence of the killer whale, Orcinus orca (Linnaeus, 1758).  Wellcome Open Research, 7:250.
  * Gao et al., (2022) Gao, S., Yang, X., Zhao, X., Wang, B., and Ye, K. (2022).  HiCAT: A tool for automatic annotation of centromere structure.  bioRxiv.
  * Hon et al., (2020) Hon, T., Mars, K., Young, G., Tsai, Y.-C., Karalius, J. W., Landolin, J. M., Maurer, N., Kudrna, D., Hardigan, M. A., Steiner, C. C., Knapp, S. J., Ware, D., Shapiro, B., Peluso, P., and Rank, D. R. (2020).  Highly accurate long-read HiFi sequencing data for five complex genomes.  Sci Data, 7:399.
  * Khost et al., (2017) Khost, D. E., Eickbush, D. G., and Larracuente, A. M. (2017).  Single-molecule sequencing resolves the detailed structure of complex satellite dna loci in drosophila melanogaster.  Genome Res, 27:709–721.
  * Kokot et al., (2017) Kokot, M., Dlugosz, M., and Deorowicz, S. (2017).  KMC 3: counting and manipulating k-mer statistics.  Bioinformatics, 33(17):2759–2761.
  * Kunyavskaya et al., (2022) Kunyavskaya, O., Dvorkina, T., Bzikadze, A. V., Alexandrov, I. A., and Pevzner, P. A. (2022).  Automated annotation of human centromeres with HORmon.  Genome Res, 32:1137–1151.
  * Könyves et al., (2022) Könyves, K., Mian, S., Johns, J., Royal Botanic Garden Edinburgh Genome Acquisition Lab, Royal Botanic Gardens Kew Genome Acquisition Lab, Darwin Tree of Life Barcoding collective, Wellcome Sanger Institute Tree of Life programme, Wellcome Sanger Institute Scientific Operations: DNA Pipelines collective, Tree of Life Core Informatics collective, Ruhsam, M., Leitch, I. J., and Darwin Tree of Life Consortium (2022).  The genome sequence of the apple, Malus domestica (Suckow) Borkh., 1803.  Wellcome Open Research, 7:297.
  * Lawniczak et al., (2021) Lawniczak, M. K., Darwin Tree of Life Barcoding collective, Wellcome Sanger Institute Tree of Life programme, Wellcome Sanger Institute Scientific Operations: DNA Pipelines collective, Tree of Life Core Informatics collective, and Darwin Tree of Life Consortium (2021).  The genome sequence of the spiny starfish, Marthasterias glacialis (Linnaeus, 1758).  Wellcome Open Research, 6:295.
  * Lawniczak et al., (2022) Lawniczak, M. K., Darwin Tree of Life Barcoding collective, Wellcome Sanger Institute Tree of Life programme, Wellcome Sanger Institute Scientific Operations: DNA Pipelines collective, Tree of Life Core Informatics collective, and Darwin Tree of Life Consortium (2022).  The genome sequence of the blue-rayed limpet, Patella pellucida Linnaeus, 1758.  Wellcome Open Research, 7:126.
  * Li, (2018) Li, H. (2018).  Minimap2: pairwise alignment for nucleotide sequences.  Bioinformatics, 34(18):3094–3100.
  * Li, (2019) Li, H. (2019).  Identifying centromeric satellites with dna-brnn.  Bioinformatics, 35:4408–4410.
  * Logsdon et al., (2021) Logsdon, G. A., Vollger, M. R., Hsieh, P., Mao, Y., Liskovykh, M. A., Koren, S., Nurk, S., Mercuri, L., Dishuck, P. C., Rhie, A., de Lima, L. G., Dvorkina, T., Porubsky, D., Harvey, W. T., Mikheenko, A., Bzikadze, A. V., Kremitzki, M., Graves-Lindsay, T. A., Jain, C., Hoekzema, K., Murali, S. C., Munson, K. M., Baker, C., Sorensen, M., Lewis, A. M., Surti, U., Gerton, J. L., Larionov, V., Ventura, M., Miga, K. H., Phillippy, A. M., and Eichler, E. E. (2021).  The structure, function and evolution of a complete human chromosome 8\.  Nature, 593:101–107.
  * Lohse et al., (2021) Lohse, K., García-Berro, A., Talavera, G., Darwin Tree of Life Barcoding collective, Wellcome Sanger Institute Tree of Life programme, Wellcome Sanger Institute Scientific Operations: DNA Pipelines collective, Tree of Life Core Informatics collective, and Darwin Tree of Life Consortium (2021).  The genome sequence of the red admiral, Vanessa atalanta (Linnaeus, 1758).  Wellcome Open Research, 6:356.
  * Lower et al., (2018) Lower, S. S., McGurk, M. P., Clark, A. G., and Barbash, D. A. (2018).  Satellite dna evolution: old ideas, new approaches.  Curr Opin Genet Dev, 49:70–78.
  * Melters et al., (2013) Melters, D. P., Bradnam, K. R., Young, H. A., Telis, N., May, M. R., Ruby, J. G., Sebra, R., Peluso, P., Eid, J., Rank, D., et al. (2013).  Comparative analysis of tandem repeats from hundreds of species reveals unique insights into centromere evolution.  Genome Biol, 14:R10.
  * Miga, (2019) Miga, K. H. (2019).  Centromeric satellite dnas: Hidden sequence variation in the human population.  Genes (Basel), 10.
  * Miga et al., (2014) Miga, K. H., Newton, Y., Jain, M., Altemose, N., Willard, H. F., and Kent, W. J. (2014).  Centromere reference models for human chromosomes x and y satellite arrays.  Genome Res, 24:697–707.
  * Naish et al., (2021) Naish, M., Alonge, M., Wlodzimierz, P., Tock, A. J., Abramson, B. W., Schmücker, A., Mandáková, T., Jamge, B., Lambing, C., Kuo, P., Yelina, N., Hartwick, N., Colt, K., Smith, L. M., Ton, J., Kakutani, T., Martienssen, R. A., Schneeberger, K., Lysak, M. A., Berger, F., Bousios, A., Michael, T. P., Schatz, M. C., and Henderson, I. R. (2021).  The genetic and epigenetic landscape of the arabidopsis centromeres.  Science, 374:eabi7489.
  * Newman et al., (2022) Newman, C., Tsai, M.-s., Buesching, C. D., Holland, P. W. H., Macdonald, D. W., Darwin Tree of Life Consortium, University of Oxford and Wytham Woods Genome Acquisition Lab, Wellcome Sanger Institute Tree of Life programme, Wellcome Sanger Institute Scientific Operations: DNA Pipelines collective, and Tree of Life Core Informatics collective (2022).  The genome sequence of the European badger, Meles meles (Linnaeus, 1758).  Wellcome Open Research, 7:239.
  * Novák et al., (2017) Novák, P., Ávila Robledillo, L., Koblížková, A., Vrbová, I., Neumann, P., and Macas, J. (2017).  TAREAN: a computational tool for identification and characterization of satellite dna from unassembled short reads.  Nucleic Acids Res, 45:e111.
  * Nurk et al., (2022) Nurk, S., Koren, S., Rhie, A., Rautiainen, M., Bzikadze, A. V., Mikheenko, A., Vollger, M. R., Altemose, N., Uralsky, L., Gershman, A., et al. (2022).  The complete sequence of a human genome.  Science, 376:44–53.
  * Pemberton et al., (2021) Pemberton, J., Johnston, S. E., Fletcher, T. J., Darwin Tree of Life Barcoding collective, Wellcome Sanger Institute Tree of Life programme, Wellcome Sanger Institute Scientific Operations: DNA Pipelines collective, Tree of Life Core Informatics collective, and Darwin Tree of Life Consortium (2021).  The genome sequence of the red deer, Cervus elaphus Linnaeus 1758\.  Wellcome Open Research, 6:336.
  * Pertile et al., (2009) Pertile, M. D., Graham, A. N., Choo, K. H. A., and Kalitsis, P. (2009).  Rapid evolution of mouse y centromere repeat dna belies recent sequence stability.  Genome Res, 19:2202–13.
  * Rabanal et al., (2022) Rabanal, F. A., Gräff, M., Lanz, C., Fritschi, K., Llaca, V., Lang, M., Carbonell-Bejerano, P., Henderson, I., and Weigel, D. (2022).  Pushing the limits of hifi assemblies reveals centromere diversity between two arabidopsis thaliana genomes.  Nucleic Acids Res, 50:12309–12327.
  * Rautiainen et al., (2023) Rautiainen, M., Nurk, S., Walenz, B. P., Logsdon, G. A., Porubsky, D., Rhie, A., Eichler, E. E., Phillippy, A. M., and Koren, S. (2023).  Telomere-to-telomere assembly of diploid chromosomes with verkko.  Nat Biotechnol.
  * Rhie et al., (2021) Rhie, A., McCarthy, S. A., Fedrigo, O., Damas, J., Formenti, G., Koren, S., Uliano-Silva, M., Chow, W., Fungtammasan, A., Kim, J., Lee, C., Ko, B. J., Chaisson, M., Gedman, G. L., Cantin, L. J., Thibaud-Nissen, F., Haggerty, L., Bista, I., Smith, M., Haase, B., Mountcastle, J., Winkler, S., Paez, S., Howard, J., Vernes, S. C., Lama, T. M., Grutzner, F., Warren, W. C., Balakrishnan, C. N., Burt, D., George, J. M., Biegler, M. T., Iorns, D., Digby, A., Eason, D., Robertson, B., Edwards, T., Wilkinson, M., Turner, G., Meyer, A., Kautt, A. F., Franchini, P., Detrich, 3rd, H. W., Svardal, H., Wagner, M., Naylor, G. J. P., Pippel, M., Malinsky, M., Mooney, M., Simbirsky, M., Hannigan, B. T., Pesout, T., Houck, M., Misuraca, A., Kingan, S. B., Hall, R., Kronenberg, Z., Sović, I., Dunn, C., Ning, Z., Hastie, A., Lee, J., Selvaraj, S., Green, R. E., Putnam, N. H., Gut, I., Ghurye, J., Garrison, E., Sims, Y., Collins, J., Pelan, S., Torrance, J., Tracey, A., Wood, J., Dagnew, R. E., Guan, D., London, S. E., Clayton, D. F., Mello, C. V., Friedrich, S. R., Lovell, P. V., Osipova, E., Al-Ajli, F. O., Secomandi, S., Kim, H., Theofanopoulou, C., Hiller, M., Zhou, Y., Harris, R. S., Makova, K. D., Medvedev, P., Hoffman, J., Masterson, P., Clark, K., Martin, F., Howe, K., Flicek, P., Walenz, B. P., Kwak, W., Clawson, H., Diekhans, M., Nassar, L., Paten, B., Kraus, R. H. S., Crawford, A. J., Gilbert, M. T. P., Zhang, G., Venkatesh, B., Murphy, R. W., Koepfli, K.-P., Shapiro, B., Johnson, W. E., Di Palma, F., Marques-Bonet, T., Teeling, E. C., Warnow, T., Graves, J. M., Ryder, O. A., Haussler, D., O’Brien, S. J., Korlach, J., Lewin, H. A., Howe, K., Myers, E. W., Durbin, R., Phillippy, A. M., and Jarvis, E. D. (2021).  Towards complete and error-free genome assemblies of all vertebrate species.  Nature, 592(7856):737–746.
  * Ruhsam et al., (2022) Ruhsam, M., Bell, D., Hart, M., Hollingsworth, P., Royal Botanic Garden Edinburgh Genome Acquisition Lab, Darwin Tree of Life Barcoding collective, Wellcome Sanger Institute Tree of Life programme, Wellcome Sanger Institute Scientific Operations: DNA Pipelines collective, Tree of Life Core Informatics collective, and Darwin Tree of Life Consortium (2022).  The genome sequence of the European crab apple, Malus sylvestris (L.) Mill., 1768.  Wellcome Open Research, 7:296.
  * Schneider et al., (2017) Schneider, V. A., Graves-Lindsay, T., Howe, K., Bouk, N., Chen, H.-C., Kitts, P. A., Murphy, T. D., Pruitt, K. D., Thibaud-Nissen, F., Albracht, D., Fulton, R. S., Kremitzki, M., Magrini, V., Markovic, C., McGrath, S., Steinberg, K. M., Auger, K., Chow, W., Collins, J., Harden, G., Hubbard, T., Pelan, S., Simpson, J. T., Threadgold, G., Torrance, J., Wood, J. M., Clarke, L., Koren, S., Boitano, M., Peluso, P., Li, H., Chin, C.-S., Phillippy, A. M., Durbin, R., Wilson, R. K., Flicek, P., Eichler, E. E., and Church, D. M. (2017).  Evaluation of GRCh38 and de novo haploid genome assemblies demonstrates the enduring quality of the reference assembly.  Genome Res, 27:849–864.
  * Sevim et al., (2016) Sevim, V., Bashir, A., Chin, C.-S., and Miga, K. H. (2016).  Alpha-CENTAURI: assessing novel centromeric repeat sequence variation with long read sequencing.  Bioinformatics, 32:1921–1924.
  * Shatskikh et al., (2020) Shatskikh, A. S., Kotov, A. A., Adashev, V. E., Bazylev, S. S., and Olenina, L. V. (2020).  Functional significance of satellite dnas: Insights from drosophila.  Front Cell Dev Biol, 8:312.
  * Sinding et al., (2021) Sinding, M.-H. S., Gopalakrishnan, S., Raundrup, K., Dalén, L., Threlfall, J., Darwin Tree of Life Barcoding collective, Wellcome Sanger Institute Tree of Life programme, Wellcome Sanger Institute Scientific Operations: DNA Pipelines collective, Tree of Life Core Informatics collective, Darwin Tree of Life Consortium, and Gilbert, T. (2021).  The genome sequence of the grey wolf, Canis lupus Linnaeus 1758.  Wellcome Open Research, 6:310.
  * Thakur et al., (2021) Thakur, J., Packiaraj, J., and Henikoff, S. (2021).  Sequence, chromatin and evolution of satellite dna.  Int J Mol Sci, 22.
  * Wei et al., (2014) Wei, K. H.-C., Grenier, J. K., Barbash, D. A., and Clark, A. G. (2014).  Correlated variation and population differentiation in satellite DNA abundance among lines of drosophila melanogaster.  Proc Natl Acad Sci U S A, 111:18793–8.
  * Wood et al., (2022) Wood, C., Bishop, J., Harley, J., Mrowicki, R., Marine Biological Association Genome Acquisition Lab, Darwin Tree of Life Barcoding collective, Wellcome Sanger Institute Tree of Life programme, Wellcome Sanger Institute Scientific Operations: DNA Pipelines collective, Tree of Life Core Informatics collective, and Darwin Tree of Life Consortium (2022).  The genome sequence of the orange-striped anemone, Diadumene lineata (Verrill, 1869).  Wellcome Open Research, 7:93.
  * Wright et al., (2022) Wright, R., Woof, K., Douglas, B., Gaya, E., Royal Botanic Gardens Kew Genome Acquisition Lab, Darwin Tree of Life Barcoding collective, Wellcome Sanger Institute Tree of Life programme, Wellcome Sanger Institute Scientific Operations: DNA Pipelines collective, Tree of Life Core Informatics collective, and Darwin Tree of Life Consortium (2022).  The genome sequence of the chicken of the woods fungus, Laetiporus sulphureus (Bull.) Murrill, 1920.  Wellcome Open Research, 7:83.
  * Zerbino and Birney, (2008) Zerbino, D. R. and Birney, E. (2008).  Velvet: algorithms for de novo short read assembly using de bruijn graphs.  Genome Res, 18:821–9.
  * Zhang et al., (2019) Zhang, L., Hu, J., Han, X., Li, J., Gao, Y., Richards, C. M., Zhang, C., Tian, Y., Liu, G., Gul, H., Wang, D., Tian, Y., Yang, C., Meng, M., Yuan, G., Kang, G., Wu, Y., Wang, K., Zhang, H., Wang, D., and Cong, P. (2019).  A high-quality apple genome assembly reveals the association of a retrotransposon and red fruit colour.  Nat Commun, 10(1):1494.


# Bounding the difference between model predictive control and neural networks

Ross Drummond<EMAIL_ADDRESS>
Stephen R. Duncan<EMAIL_ADDRESS>
Department of Engineering Science  University of Oxford  17 Parks Road  OX1
3PJ  Oxford  United Kingdom
Matthew C. Turner<EMAIL_ADDRESS>
Department of Electronics and Computer Science  University of Southampton
Southampton  UK  SO17 1BJ  Patricia Pauli<EMAIL_ADDRESS>
Frank Allgöwer<EMAIL_ADDRESS>
Institute for Systems Theory and Automatic Control  University of Stuttgart
70569 Stuttgart  Germany

###### Abstract

There is a growing debate on whether the future of feedback control systems
will be dominated by data-driven or model-driven approaches. Each of these two
approaches has their own complimentary set of advantages and disadvantages,
however, only limited attempts have, so far, been developed to bridge the gap
between them. To address this issue, this paper introduces a method to bound
the worst-case error between feedback control policies based upon model
predictive control (MPC) and neural networks (NNs). This result is leveraged
into an approach to automatically synthesize MPC policies minimising the
worst-case error with respect to a NN. Numerical examples highlight the
application of the bounds, with the goal of the paper being to encourage a
more quantitative understanding of the relationship between data-driven and
model-driven control.

###### keywords:

Neural network robustness, model predictive control, semi-definite
programming.

##  1 Introduction

Two distinct approaches for the feedback control of dynamical systems are
emerging, categorised into _model-driven_ and _data-driven_ methods. Model-
driven approaches rely upon a physical model of the system to make predictions
about its future behaviour and then compute a control input to shape this
future response ([Green and Limebeer (2012)]). Data-driven methods, in
particular model-free reinforcement learning, forgo the need of a predictive
model and instead rely upon identifying patterns within the trajectory data
itself to determine control policies ([Sutton and Barto (2018)]). The benefits
and limitations of these two approaches are also quite distinct. Model-driven
approaches are often equipped with performance certificates for properties
such as closed-loop stability, optimality and robustness ([Green and Limebeer
(2012)]), with these certificates making them more suited for safety-critical
systems. But model-driven methods generally struggle in applications where
existing models remain underdeveloped and when applied to large-scale systems
as computing the control polices can become intractable ([Šiljak and Zečević
(2005)]).

Data-driven methods, by contrast, do not require an accurate model of the
system ([Sutton and Barto (2018)]) and can find control policies which would
be unobtainable using model-based approaches (as they avoid the biases
inherent within models), a point illustrated by the famous “move 37” from
AlphaGo ([Silver et al. (2017)]) which confounded experts’ analysis. A main
drawback of most data-driven methods is their general lack of rigorous
guarantees, e.g. robustness (Szegedy et al., 2013), as small changes in their
inputs can lead to large changes in their outputs, and explaining their
behaviour remains an open problem. Robustness issues are a primary reason why
we have yet to fully see the heralded successes of data-driven methods
transfer from controlled, laboratory settings to real-life practical
applications.

The distinction between model and data driven control methods has driven
interest in bridging the gap between them. In particular, there is a need to
better characterise the relationship between each approach’s benchmark methods
which are, arguably, model predictive control (MPC) for model-based control
and deep neural networks (NNs) for the data-driven approach. Relating these
two approaches will not only bring greater understanding but will also
encourage the development of new control methods that combine each approach’s
complementary advantages.

#### Contribution:

The two main results of this paper are:

  1. 1.

A method to bound the difference between a given MPC and a NN controller
(Theorem 4.1).

  2. 2.

A method to synthesise MPCs that robustly approximate neural networks (Theorem
5.1).

The overriding theme of these results is to quantify the similarity between
MPC and neural network control actions, with the bounds holding robustly, in
the sense that they hold for fairly generic neural network architectures and
MPC cost functions considering pre-defined hyper-rectangles as input
constraints. The bounds are obtained by solving a semi-definite program (SDP),
a class of convex optimisation programs which have been extensively used for
analysing the robustness of both NNs and MPC, as in Fazlyab et al. (2020); Li
et al. (2006).

#### Literature:

Several existing results have also been concerned with relating the disparate
control schemes of MPC and NNs. Of particular note are the results on explicit
MPC where the user takes advantage of the fact that the optimal solution to
the linear MPC problem can be expressed as a piece-wise linear function
(Bemporad et al. (2002)). By identifying which region of this piece-wise
linear function the current state of the system is in, known as solving the
“point location problem”, then the MPC action can be obtained without online
computation of the finite horizon optimal control problem. The main drawback
of explicit MPC is that solving the point location problem can be challenging
(since the complexity of the piece-wise affine function, measured by the
number of piece-wise regions, scales exponentially with the state dimension
and the number of constraints, Alessio and Bemporad (2009)), which has limited
its use in practice. To avoid this computational bottleneck, it has been
proposed to approximate the piece-wise linear optimal MPC policy with a neural
network, Parisini and Zoppoli (1995); Zhang et al. (2019), as the
approximating NN will also be a piece-wise linear function if a ReLU function
is used for the activation function green!50!blue(Darup, 2020). Evaluating a
NN is generally much simpler than solving the point location problem but can
introduce an error which has proven difficult to bound. The bounds formulated
in this paper are able to quantify the worst-case error of this approximation.

The presented results can be understood within the more general context of
applying control theory to verify the robustness of neural networks. The
theoretical framework for these methods are now classic, emerging from studies
such as Chu and Glover (1999), Barabanov and Prokhorov (2002), Angeli (2009),
and have been extended in recent years in studies such as Fazlyab et al.
(2020). The presented results contrast with these existing robustness problems
in that they are concerned with the robustness of a neural network with
respect to MPC, not to itself. Some of the authors used similar ideas to
quantify the similarity of one neural network with respect to another (Li et
al. (2021)), which allowed the approximation errors of pruned, quantised and
reduced-order neural networks to be bounded (Li et al. (2021); Drummond et al.
(2021)). The recently released study Fabiani and Goulart (2021) which was also
concerned with bounding the error between MPC and a NN must also be
highlighted. The method developed in that paper involved solving a mixed
integer optimisation problem to bound the error between MPC and NNs with ReLU
activation functions in particular. In contrast, the results of this paper
solve a SDP to verify the robustness of fairly generic NN architectures with
respect to MPC. Finally, it is highlighted how this paper is concerned with
quantifying the similarity between MPC and a NN but it is known, that, in
certain circumstances, these two policies can be shown to be equivalent, Darup
(2020).

Examples of NN robustness failures (Szegedy et al. (2013)) have driven
mounting interest in using semi-definite programming for NN analysis and
synthesis Fazlyab et al. (2020). These results include the estimation of
Lipschitz constants for NN mappings (Fazlyab et al., 2019; Pauli et al.,
2021b) and the stability analysis of closed-loop systems using NN controllers
Yin et al. (2021a). The stability analysis problem in particular can be linked
to the problem of absolute stability (Desoer and Vidyasagar (2009); Barabanov
and Prokhorov (2002); Chu and Glover (1999)), as a feed-forward neural network
is a static, memoryless nonlinearity Pauli et al. (2021a). Within this
context, Yin et al. (2021b) proposed approximating a MPC using NN controllers
while satisfying LMI constraints to ensure closed-loop stability. Enforcing
stability on a large-scale structure like a NN can be challenging, which
motivates a workaround that finds and enforces bounds relating the two
controller types. Doing so may lead to more elegant stability constraints to
be formulated, as meeting these bounds have already been shown to imply
stability guarantees Hertneck et al. (2018).

### Notation

The set of real numbers is denoted $\mathbb{R}$ and the natural numbers are
$\mathbb{N}$. For dimension $n$, the set of real vectors is $\mathbb{R}^{n}$,
the set of vectors with non-negative elements is denoted $\mathbb{R}^{n}_{+}$
and the zero vector is ${\bm{0}_{n}}$. Element-wise positivity of a vector is
denoted $>$ while positive definiteness of a matrix is denoted $\succ$, with
similar definitions for negativity and semi-definiteness of matrices. The set
of real matrices of dimension $n\times m$ is $\mathbb{R}^{n\times m}$ and the
zero matrix is $\bm{0}_{n\times m}$. For dimension $n$, the space of positive
(semi-)definite matrices is ($\mathcal{S}^{n}_{\succeq 0}$)
$\mathcal{S}^{n}_{\succ 0}$, non-negative diagonal matrices is
$\mathbb{D}^{n}_{\succeq 0}$ and the identity matrix is $I_{n}$. We use the
$\star$ notation for symmetric matrices, as in $\begin{bmatrix}A&B\\\
B^{\top}&C\end{bmatrix}=\begin{bmatrix}A^{\top}&B^{\top}\\\
B&C^{\top}\end{bmatrix}=\begin{bmatrix}A&B\\\ \star&C\end{bmatrix}.$ Several
of the proofs of the lemmas are immediate and so are not shown.

##  2 Problem Set-up

This section introduces the general problem considered in the paper of
computing bounds and approximations between MPCs and NNs.

###  2.1 Model predictive control

Consider input-constrained model predictive control with linear dynamics and a
quadratic cost.

######  Definition 2.1.

Consider a system with linear dynamics

$\displaystyle x[k+1]=Ax[k]+Bu[k],$ (1)

with state $x[k]\in\mathbb{R}^{n_{x}}$ and input $u[k]\in\mathbb{R}$. For some
finite horizon $N\in\mathbb{N}$, the input constrained MPC control law with
linear dynamics is defined as the solution to

$\displaystyle z_{\text{MPC}}(x[k])=$
$\displaystyle\,\text{arg}\,\min_{u}\begin{bmatrix}x[k+1]\\\ \vdots\\\
x[k+N]\end{bmatrix}^{\top}Q\begin{bmatrix}x[k+1]\\\ \vdots\\\
x[k+N]\end{bmatrix}+\begin{bmatrix}u[k]\\\ \vdots\\\
u[k+N-1]\end{bmatrix}^{\top}R\begin{bmatrix}u[k]\\\ \vdots\\\
u[k+N-1]\end{bmatrix},$ (2a) $\displaystyle\text{subject to }Lu\leq b,~{}b\geq
0,$ (2b)

for some $Q\in\mathbb{S}^{Nn_{x}}_{\succeq 0}$, $R\in\mathbb{S}^{N}_{\succeq
0}$, $L\in\mathbb{R}^{n_{in}\times N}$ and $b\in\mathbb{R}^{n_{in}}_{+}$.

The linear dynamics of (1) mean that the optimisation problem (2) can be
converted into a quadratic program (QP).

######  Definition 2.2.

The MPC control law defined by Definition 2.1 can be expressed as the solution
to a constrained quadratic program

$\displaystyle z_{\text{MPC}}(x[k])=$
$\displaystyle~{}\text{arg}\,\min_{u}u^{\top}Hu+u^{\top}hx[k]$ (3a)
$\displaystyle\text{subject to }Lu\leq b,~{}b\geq 0,$ (3b)

for some $H\in\mathcal{S}^{N}_{\succ 0}$, $h\in\mathbb{R}^{N\times n_{x}}$
parameterised by the matrices ($A,\,B,\,Q,\,R$) of Definition 2.1 and
$b\in\mathbb{R}^{n_{in}}_{+}$. The first instant of the QP solution from (3)
is taken as the MPC control action, as in

$\displaystyle
u_{\text{MPC}}(x[k])=\begin{bmatrix}1,&{\bm{0}_{N-1}}^{\top}\end{bmatrix}z_{\text{MPC}}(x[k])=C_{\text{MPC}}z_{\text{MPC}}(x[k]).$
(4)

###  2.2 Neural networks

The class of neural networks considered in this paper are those which can be
structured in an implicit form (El Ghaoui et al. (2021)).

######  Definition 2.3.

The considered class of neural networks
$u_{\text{NN}}(x[k]):\mathbb{R}^{n_{x}}\to\mathbb{R}$ are those defined by the
implicit form

$\displaystyle{z}_{\text{NN}}$
$\displaystyle=\phi(W{z}_{\text{NN}}+W_{0}x[k]+\beta),\quad\phi(\cdot):\mathbb{R}^{M}\mapsto\mathbb{R}^{M},$
(5a) $\displaystyle u_{\text{NN}}(x[k])$
$\displaystyle=C_{\text{NN}}{z}_{\text{NN}}+D_{\text{NN}},$ (5b) where
$C_{\text{NN}}=\begin{bmatrix}\bm{0}_{n_{\ell}}^{\top},&\ldots\,,&\bm{0}_{n_{\ell}}^{\top},&W^{o}\end{bmatrix},$
$D_{\text{NN}}\in\mathbb{R}$ and $\displaystyle{z}_{\text{NN}}$
$\displaystyle=\begin{bmatrix}z^{1}_{\text{NN}}\\\ z^{2}_{\text{NN}}\\\
\vdots\\\
z^{\ell}_{\text{NN}}\end{bmatrix},\,W=\begin{bmatrix}W^{1,1}&W^{1,2}&\ldots&W^{1,\ell}\\\
W^{2,1}&W^{2,2}&\ddots&\vdots\\\ \vdots&\ddots&\ddots&W^{\ell-1,\ell}\\\
W^{\ell,1}&\dots&W^{\ell,\ell-1}&W^{\ell,\ell}\end{bmatrix},\,W_{0}=\begin{bmatrix}W^{1,0}\\\
W^{2,0}\\\ \vdots\\\
W^{\ell,0}\end{bmatrix},\,\beta=\begin{bmatrix}\beta^{0}\\\ \beta^{1}\\\
\vdots\\\ \beta^{\ell}\end{bmatrix}.$ (5c)

In the above, $\phi(\cdot)$ are the nonlinear activation functions which could
be any function satisfying one or more of the function properties of
Definition 7.1 (Appendix 1), acting element-wise upon their arguments. The NN
will be assumed to have $\ell$ hidden layers and, for the sake of notational
simplicity, it will be assumed that each layer will be of dimension
$n_{\ell}$, as in $z^{j}_{\text{NN}}\in\mathbb{R}^{n_{\ell}}$,
$W^{i,j}\in\mathbb{R}^{n_{\ell}\times n_{\ell}}$,
$W^{i,0}\in\mathbb{R}^{n_{\ell}\times n_{x}}$,
$\beta^{j}\in\mathbb{R}^{n_{\ell}}$ and
${W^{o}}^{\top}\in\mathbb{R}^{n_{\ell}}$ for $i=1,\,\dots,\,\ell$,
$j=1,\,\dots,\,\ell$. The total length of the vector
$z_{\text{NN}}\in\mathbb{R}^{M}$ is defined as $M=n_{\ell}\times\ell$.

Adopting the implicit NN structure of (5) in Definition 2.3, also known as an
equilibrium network Revay et al. (2020); Bai et al. (2019), helps simplify the
notation whilst also ensuring the results are applicable to a wide class of NN
architectures, including deep recurrent and feed-forward NNs.

###  2.3 Bounds between MPC and NN

Bounds quantifying the similarity between the NN and MPC control policies
described in Definitions 2.2 and 2.3 are sought. These bounds are formalised
as solutions to the following problem.

######  Problem 2.4.

For any state $x[k]\in\mathcal{X}$, find $\gamma_{x}\geq 0$ and $\gamma\geq 0$
such that the error between the neural network and MPC control actions, namely
$u_{\text{NN}}(x[k])-u_{\text{MPC}}(x[k])$, is bounded by

$\displaystyle\|u_{\text{NN}}(x[k])-u_{\text{MPC}}(x[k])\|_{2}^{2}\leq\gamma_{x}\|x[k]\|_{2}^{2}+\gamma,\quad\forall
x[k]\in\mathcal{X}.$ (6)

######  Remark 2.5.

If both controllers enforce $x[k]=0$ as an equilibrium point of the dynamical
system, then it is possible to obtain $\gamma=0$ since the two control actions
should be equivalent at that point.

##  3 Quadratic Constraints

For the bounds of Problem 2.4 to hold robustly, as in for all inputs $x[k]$ in
some set $\mathcal{X}$, characterisations of both the nonlinear activation
functions of the neural network and the solution structure of the MPC control
law have to be included within the problem formulation. Here, this information
is incorporated using quadratic constraints.

###  3.1 Quadratic constraints for the neural network

The first step of this process is to determine the quadratic constraint of the
nonlinear activation functions $\phi(.)$. Many different activation functions
have been implemented within neural networks, including the tanh and ReLU
functions, with each having their own characteristics. In general, each of
these candidate activation functions satisfies certain properties which can be
incorporated within a robust optimisation framework to compute the error
bounds.

Examples of these function properties are given in Definition 7.1 (Appendix
1). Most commonly adopted activation functions, including the tanh and ReLU,
satisfy at least one of these properties. For example, the tanh function is
slope-restricted, bounded and sector bounded, the ReLU function is slope-
restricted while both it and its complement are positive and satisfy the
complementarity conditions.

By satisfying one or more of these function properties, the activation
functions also satisfy quadratic constraints, as detailed in Lemma 7.2
(Appendix 2). These quadratic constraints are defined with the multipliers
${\bf T}^{i}$
($i\in\left\\{\text{s},\text{sl},+,c+,B,0,\times,\otimes\right\\}$), and in
the following, the notation ${\bf T}\in\mathbb{T}$ will be used to collect all
of the relevant ${\bf T}^{i}$’s
($i\in\left\\{\text{s},\text{sl},+,c+,B,0,\times,\otimes\right\\}$) from Lemma
7.2 (Appendix 2), with the set $\mathbb{T}$ characterising the positivity
conditions of the ${\bf T}^{i}$’s.

If the neural network’s activation function satisfies more than one of these
quadratic constraints, as in it satisfies one or more of the properties in
Lemma 7.2 (Appendix 2), then each of these quadratic constraints can be
combined into a single quadratic constraint defined by a matrix $\lambda$, as
considered previously in results such as Drummond et al. (2021); Fazlyab et
al. (2020). In this way, information about the nonlinear activation functions
of the neural network can be incorporated into the robust optimisation problem
for the error bounds.

######  Lemma 3.1.

If the activation function $\phi(\cdot)$ satisfies one or more of the
quadratic constraints of Lemma 7.2 (Appendix 2), then there exists matrices

$\displaystyle\hat{\lambda}({\bf
T})=\begin{bmatrix}\hat{\lambda}_{11}&\hat{\lambda}_{12}&\hat{\lambda}_{13}\\\
\star&\hat{\lambda}_{22}&\hat{\lambda}_{23}\\\
\star&\star&\hat{\lambda}_{33}\end{bmatrix},\quad{\lambda}({\bf
T})=\begin{bmatrix}{\lambda}_{11}&{\lambda}_{12}&{\lambda}_{13}\\\
\star&{\lambda}_{22}&{\lambda}_{23}\\\
\star&\star&{\lambda}_{33}\end{bmatrix},$ (7)

defined by the ${\bf T}^{i}$’s
$(i\in\left\\{\text{s},\text{sl},+,c+,B,0,\times,\otimes\right\\})$ of Lemma
7.2 (Appendix 2), such that with the vectors

$\displaystyle\hat{\mu}(x[k])$
$\displaystyle=\begin{bmatrix}(Wz_{\text{NN}}+W_{0}x[k]+\beta)^{\top},&{z_{\text{NN}}}^{\top},&1\end{bmatrix}^{\top},$
(8a) $\displaystyle{\mu}(x[k])$
$\displaystyle=\begin{bmatrix}{x[k]}^{\top},&{z_{\text{NN}}}^{\top},&1\end{bmatrix}^{\top},$
(8b)

the following inequality holds

$\displaystyle s_{\text{NN}}(x[k])=\hat{\mu}(x[k])^{T}\hat{\lambda}({\bf
T})\hat{\mu}(x[k])=\mu(x[k])^{T}\lambda({\bf T})\mu(x[k])\geq 0.$ (9)

######  Proof 3.2.

Immediate from applying the S-procedure (Jönsson (2001)) on Lemma 7.2
(Appendix 2).

###  3.2 Quadratic constraints for model predictive control

Following (Li et al., 2006, Result 5), a quadratic constraint can also be
characterised for the solution of an MPC control law defined by the QP of
Definition 2.2. In Li et al. (2006), this quadratic constraint was introduced
to characterise the closed-loop stability of linear systems controlled by MPC,
using Zames-Falb multipliers to provide the stability certificates (Zames and
Falb (1968); Turner and Drummond (2021)). Here, we exploit this quadratic
constraint to include information about the solution structure of the MPC
control law into the error bound computation.

######  Theorem 3.3.

[Li et al. (2006), Result 5] The MPC control law characterised by the QP of
Definition 2.2 satisfies, for any $\tau_{\text{QP}}\geq 0$,

$\displaystyle
s_{\text{QP}}(x[k])=-\tau_{\text{QP}}\left({z_{\text{MPC}}}^{\top}Hz_{\text{MPC}}+{z_{\text{MPC}}}^{\top}hx[k]\right)\geq
0.$ (10)

On top of this quadratic constraint for the QP cost function (10), another
quadratic constraint for the inequality bounds (3b) can also be devised since,
for any positive vector $q\in\mathbb{R}^{n_{in}}_{+}$,

$\displaystyle
s_{in}(x[k])=q^{\top}\left(b-Lz_{MPC}\right)+\left(b-Lz_{MPC}\right)^{\top}q\geq
0.$ (11)

The function properties of Definition 7.1 (Appendix 1) can be used to build
satisfactory $q$. For example, if $\phi(.)$ and its complement are positive,
as for the ReLU, then, for some $\tau_{in}\in\mathbb{R}_{+}^{n_{in}}$,
$(\tau_{in}^{+},\,\tau_{in}^{c+})\in\mathbb{R}_{+}^{n_{in}\times M}$, one
could define

$\displaystyle
q={\tau_{in}}+(\tau_{in}^{+}+\tau_{in}^{c+}(I_{M}-W))z_{\text{NN}}-\tau_{in}^{c+}W_{0}x[k]-\tau_{in}^{c+}\beta$
(12)

which links the NN and MPC vectors, $z_{\text{NN}}$ and $z_{\text{MPC}}$
within the quadratic constraint. The quadratic constraint for the inequality
constraints $s_{in}(x[k])$ can then be combined with that of the cost function
from Theorem 3.3 to give the following characterisation of the MPC solution.

######  Lemma 3.4.

For any $x[k]\in\mathbb{R}^{n_{x}}$, the MPC solution of the QP given in
Definition 2.2 satisfies

$\displaystyle s_{\text{MPC}}(x[k])=~{}$ $\displaystyle
s_{\text{QP}}(x[k])+s_{\text{in}}(x[k])\geq 0,$ (13)

with $s_{\text{QP}}(x[k])$ defined in (10) and $s_{\text{in}}(x[k])$ in (11).

###  3.3 Quadratic constraints for the input

Information about the input constraint set $x[k]\in\mathcal{X}$ must also be
included within the robust optimisation problem for the bounds. By restricting
the input space to a hyper-rectangle, the following lemma defines a quadratic
constraint for this set containment.

######  Definition 3.5.

For some upper $\overline{x}$ and lower $\underline{x}$ bound, define
$\mathcal{X}:=\\{x[k]:\,x[k]\geq\underline{x},\,x[k]\leq\overline{x}\\}$.

######  Lemma 3.6 (Fazlyab et al. (2020)).

If $x[k]\in\mathcal{X}$, then, for any $\tau_{x}\in\mathbb{D}_{+}^{n_{x}}$,
the following holds

$\displaystyle
s_{\mathcal{X}}(x[k])=\begin{bmatrix}x[k]^{\top}&1\end{bmatrix}\begin{bmatrix}-2\tau_{x}&\tau_{x}(\overline{x}+\underline{x})\\\
\star&-(\overline{x}^{\top}\tau_{x}\underline{x}+\underline{x}^{\top}\tau_{x}\overline{x})\end{bmatrix}\begin{bmatrix}x[k]\\\
1\end{bmatrix}\geq 0.$ (14)

##  4 Error bound between MPC and NN

The quadratic constraints of Lemmas 3.1, 3.4 and 3.6 allow information about
the nonlinear activation functions, the solution structure of the MPC
controller and the set containment $x[k]\in\mathcal{X}$ to be included within
the error bound optimisation. In this way, the error bounds of Problem 2.4 can
be computed by solving the following SDP.

######  Theorem 4.1.

Consider a given NN structure of Definition 2.3 and MPC control law of
Definition 2.2. For given weights $\omega_{x}$ and $\omega$, if there exists a
solution to

$\displaystyle\min_{{\bf
T},\,\tau_{QP},\,\tau_{in},\,\gamma_{x},\,\gamma,\,\tau_{x}}$
$\displaystyle\quad\omega_{x}\gamma_{x}+\omega\gamma$ (15a) subject to:
$\displaystyle\quad\Lambda_{\text{NN}}({\bf
T})+\Lambda_{\text{MPC}}(\tau_{QP},\tau_{in})+X(\tau_{x})+\Omega(\gamma_{x},\,\gamma)\preceq
0,$ (15b) $\displaystyle\quad{\bf
T}\in\mathbb{T},\tau_{x}\in\mathbb{D}^{n_{x}}_{+},\,\tau_{QP}\geq
0,\tau_{in}\in\mathbb{R}_{+}^{n_{in}},\,\gamma_{x}\geq 0,\gamma\geq 0,$ (15c)

where the matrices $\Lambda_{\text{NN}}({\bf T})$,
$\Lambda_{\text{MPC}}(\tau_{QP},\tau_{in})$, $X(\tau_{x})$ and
$\Omega(\gamma_{x},\,\gamma)$ are defined in (23) of Appendix 3, then

$\displaystyle\|u_{\text{NN}}(x[k])-u_{\text{MPC}}(x[k])\|_{2}^{2}\leq\gamma_{x}\|x[k]\|_{2}^{2}+\gamma,\quad\forall
x[k]\in\mathcal{X}.$ (16)

######  Proof 4.2.

Multiplying (15b) on the left by

$\displaystyle\zeta(x[k])=\begin{bmatrix}x[k]^{\top},&{z_{\text{NN}}}^{\top},&{z_{\text{MPC}}}^{\top},&1\end{bmatrix}^{\top}$
(17)

and on the right by it’s transpose implies that, when $x[k]\in\mathcal{X}$,
then

$\displaystyle
s_{NN}(x[k])+s_{\mathcal{X}}(x[k])+s_{MPC}(x[k])+\|u_{\text{NN}}(x[k])-u_{\text{MPC}}(x[k])\|_{2}^{2}-\gamma_{x}\|x[k]\|_{2}^{2}-\gamma\leq
0.$ (18)

Lemmas 3.1, 3.4 and 3.6 state that $s_{NN}(x[k])\geq 0$,
$s_{\mathcal{X}}(x[k])\geq 0$, $s_{MPC}(x[k])\geq 0$ when
$x[k]\in\mathcal{X}$. Hence, the bound (16) must hold.

##  5 Automatic synthesis of an MPC approximating a NN

Theorem 4.1 allows the error between _given_ MPC and NN control policies to be
bounded. This idea can be extended to _synthesize_ MPC control laws that
minimise the worst-case error to the NN.

######  Theorem 5.1.

Consider a given NN structure of (5). For the given weights $\omega_{x}$ and
$\omega$, if there exists a solution to

$\displaystyle\min_{{\bf T},\,H,\,h,\,L,\,b,\,\gamma_{x},\,\gamma,\,\tau_{x}}$
$\displaystyle\quad\omega_{x}\gamma_{x}+\omega\gamma$ (19a) subject to:
$\displaystyle\quad\Lambda_{\text{NN}}({\bf
T})+\check{\Lambda}_{\text{MPC}}(H,h,L,b)+X(\tau_{x})+\Omega(\gamma_{x},\,\gamma)\preceq
0,$ (19b) $\displaystyle\quad{\bf
T}\in\mathbb{T},\tau_{x}\in\mathbb{D}_{+}^{n_{x}},\,\gamma_{x}\geq
0,\gamma\geq 0,\,H\in\mathcal{S}^{N}_{\succ 0},\,h\in\mathbb{R}^{N\times
n_{x}},\,b\in\mathbb{D}_{+}^{n_{in}},L\in\mathbb{R}^{{n_{in}}\times N}$ (19c)

with the matrices $\Lambda_{\text{NN}}({\bf T})$,
$\check{\Lambda}_{\text{MPC}}(H,h,L,b)$, $X(\tau_{x})$ and
$\Omega(\gamma_{x},\,\gamma)$ defined in Appendix 3, then

$\displaystyle\|u_{\text{NN}}(x[k])-u_{\text{MPC}}(x[k])\|_{2}^{2}\leq\gamma_{x}\|x[k]\|_{2}^{2}+\gamma,\quad\forall
x[k]\in\mathcal{X}.$ (20)

######  Proof 5.2.

From Lemma 3.4, the computed $(H,h,L,b)$ define a quadratic program from
Definition 2.2 for which $s_{MPC}(x[k])\geq 0$ using $\tau_{\text{QP}}=1$ and
each element of $\tau_{\text{in}}$ being one. The proof is then immediate from
the argument of Theorem 4.1.

######  Remark 5.3.

The main advantage of using Theorem 5.1 to generate an approximating MPC is
that it includes the worst-case approximation error directly within the cost
function being minimised in (19).

Figure 1: Bounds $\gamma$ and $\gamma_{x}$ and empirical lower bound on
$\gamma_{x}$ for respectively 10 NNs with $\\{3,5,10,20,100\\}$ neurons in the
hidden layer.

##  6 Numerical examples

In the following, we apply Theorem 4.1 to bound the difference between an MPC
and a NN controller. All NNs were trained using Pytorch and all SDPs were
solved in Matlab using the parser YALMIP (Lofberg (2004)) and solver MOSEK
(Andersen and Andersen (2000)).

We consider the following linear time-invariant system

$x[k+1]=\begin{bmatrix}4/3&-2/3\\\ 1&0\end{bmatrix}x[k]+\begin{bmatrix}0\\\
1\end{bmatrix}u[k],$ (21)

and solve (3) choosing the horizon length $N=10$ and the input constraints
$-0.1\leq u_{\text{MPC}}\leq~{}0.1$ such that
$L=\begin{bmatrix}1&-1\end{bmatrix}^{\top}$,
$b=\begin{bmatrix}0.1&0.1\end{bmatrix}^{\top}$. The MPC quadratic cost
function of (2) was parameterised by

$\widetilde{P}=\begin{bmatrix}7.1667&-4.2222\\\
-4.2222&4.6852\end{bmatrix},~{}\widetilde{Q}=\begin{bmatrix}1&-2/3\\\
-2/3&3/2\end{bmatrix},~{}\widetilde{R}=1,$

yielding
$Q=\operatorname{blkdiag}(\widetilde{Q},\,\widetilde{Q},\,\dots,\widetilde{Q},\,\widetilde{P})$
and $R=\operatorname{blkdiag}(\tilde{R},\,\dots,\,\tilde{R})$. Subsequently,
feed-forward neural networks with one hidden layer and ReLU activation were
trained using back-propagation to approximate this MPC controller. The results
are detailed in Figure 1 where upper and lower bounds for the obtained
approximation error are shown. The upper bound, defined in terms of $\gamma$
and $\gamma_{x}$, was obtained by solving Theorem 4.1 while the empirical
lower bound was obtained by sampling. For different numbers of neurons in the
hidden layer, ten randomly initialised NNs were trained and their
approximation errors are plotted in the figure. For all cases, it was found
that $\gamma=0$, agreeing with the prediction of Remark 2.5. The computed
upper bound was always above the empirical lower bound, with the bounds being
relatively tight for the small NNs. The larger gap between the upper and lower
bounds for the larger NN examples is due to the increasing number of neurons
causing an increased abstraction of the neural network mapping by the
quadratic constraints, and hence an increase in the conservatism of the
problem.

## Conclusions

A method to bound the worst-case error between feedback control policies based
upon neural networks (NNs) and model-predictive control (MPC) was introduced.
Using these bounds, a method to synthesize MPC policies minimising the worst-
case error with respect to a NN was developed. The overriding goal of this
work was to improve our understanding of, and even quantify, the relationship
between MPC and NNs, helping to bridge the gap between model-driven and data-
driven control. Future work will explore using the bounds to generate closed-
loop stability certificates of both the MPC and NN controllers, reducing the
conservatism of the bounds and applying them to the verification of control
systems used in practice.

##  7 Appendix

###  7.1 Appendix 1: Activation function properties

######  Definition 7.1 (Drummond et al. (2021)).

The activation function $\phi(s):\mathcal{S}\subset\mathbb{R}\to\mathbb{R}$
satisfying $\phi(0)=0$ is said to be _sector bounded_ if
$\frac{\phi(s)}{s}~{}\in~{}[0,\delta]~{}\forall s\in\mathcal{S},~{}\delta>0,$
and _slope restricted_ if
$\frac{d\phi(s)}{ds}~{}\in~{}[\underline{\beta},\beta],~{}\forall
s\in\mathcal{S},~{}\beta>0.$ If $\underline{\beta}=0$ then the nonlinearity is
_monotonic_ and if $\phi(s)$ is slope restricted then it is also sector
bounded. The activation function $\phi(s)$ is _bounded_ if
$\phi(s)~{}\in~{}[\underline{c},\overline{c}],~{}\forall s\in\mathcal{S},$ it
is _positive_ if $\phi(s)\geq 0,~{}\forall s\in\mathcal{S},$ its _complement
is positive_ if $\phi(s)-s\geq 0,~{}\forall s\in\mathcal{S},$ and it satisfies
the _complementarity condition_ if $(\phi(s)-s)\phi(s)=0,~{}\forall
s\in\mathcal{S}.$

###  7.2 Appendix 2: Quadratic constraints for the activation functions

######  Lemma 7.2 (Drummond et al. (2021)).

Consider the vectors $y,y_{1}\in\mathbb{R}^{n_{y}}$, and
$\upsilon\in\mathbb{R}^{n_{\upsilon}}$ that are mapped component-wise through
the activation functions $\phi(\cdot):\mathbb{R}^{n_{y}}\to\mathbb{R}^{n_{y}}$
and $\tilde{\phi}(\cdot):\mathbb{R}^{n_{v}}\to\mathbb{R}^{n_{v}}$. If
$\phi(y)$ is _sector-bounded_ , then

$\displaystyle(\delta y-\phi(y))^{T}{\bf T}^{\text{s}}\phi(y)$
$\displaystyle\geq 0,\quad\forall y\in\mathbb{R}^{n_{y}},~{}{\bf
T}^{\text{s}}\in\mathbb{D}^{n_{y}}_{+};$ (22a) _slope-restricted_ then
$(\beta(y-y_{1})-(\phi(y)-\phi(y_{1}))^{T}{\bf
T}^{\text{sl}}(\phi(y)-\phi(y_{1})-\underline{\beta}(y-y_{1}))\geq
0;\quad\forall\\{y,\,y_{1}\\}\in\mathbb{R}^{n_{y}},\,{\bf
T}^{\text{sl}}\in\mathbb{D}^{n_{y}}_{+};$ (22b) _bounded_ then
$\displaystyle(\overline{c}-\phi(y))^{T}{\bf T}^{B}(\phi(y)-\underline{c})$
$\displaystyle\geq 0,\quad\forall y\in\mathbb{R}^{n_{y}},\,{\bf
T}^{B}\in\mathbb{D}^{n_{y}}_{+};$ (22c) _positive_ then $\displaystyle({\bf
T}^{+})^{T}\phi(y)\geq 0,\quad\forall y\in\mathbb{R}^{n_{y}},\,{\bf
T}^{+}\in\mathbb{R}_{+}^{n_{y}};$ (22d) such that its _complement is positive_
then $\displaystyle({\bf T}^{c+})^{T}(\phi(y)-y)\geq 0,\quad\forall
y\in\mathbb{R}^{n_{y}},\,{\bf T}^{c+}\in\mathbb{R}^{n_{y}}_{+}.$ (22e) If
$\phi(y)$ satisfies the _complementary_ condition then
$\displaystyle(\phi(y)-y)^{T}{\bf T}^{0}\phi(y)=0,\quad\forall
y\in\mathbb{R}^{n_{y}},\,{\bf T}^{0}\in\mathbb{D}^{n_{y}}.$ (22f)
Additionally, if both $\phi(y)$ and $\tilde{\phi}(\upsilon)$ and their
complements are positive then so are the _cross terms_
$\displaystyle\small\tilde{\phi}(\upsilon)^{T}{\bf T}^{\times}(\phi(y)-y)\geq
0,\;\forall\upsilon\in\mathbb{R}^{n_{v}},y\in\mathbb{R}^{n_{y}},\,{\bf
T}^{\times}\in\mathbb{R}_{+}^{n_{\upsilon}\times n_{y}},$ (22g)
$\displaystyle\tilde{\phi}(\upsilon)^{T}{\bf T}^{\otimes}\phi(y)\geq
0,\;\forall\upsilon\in\mathbb{R}^{n_{v}},y\in\mathbb{R}^{n_{y}},\,{\bf
T}^{\otimes}\in\mathbb{R}_{+}^{n_{\upsilon}\times n_{y}}.$ (22h)

###  7.3 Appendix 3: Matrix definitions

$\displaystyle\Omega=\begin{bmatrix}-\gamma_{x}I_{n_{x}}&\bm{0}_{n_{x}\times
M}&\bm{0}_{n_{x}\times N}&\bm{0}_{n_{x}}\\\
\star&{C_{\text{NN}}}^{\top}C_{\text{NN}}&-{C_{\text{NN}}}^{\top}C_{\text{MPC}}&\bm{0}_{M}\\\
\star&\star&{C_{\text{MPC}}}^{\top}C_{\text{MPC}}&-{C_{\text{MPC}}}^{\top}{D_{\text{NN}}}\\\
\star&\star&\star&-\gamma\end{bmatrix},~{}\Lambda_{NN}=\begin{bmatrix}{\lambda}_{11}&{\lambda}_{12}&\bm{0}_{n_{x}\times
N}&{\lambda}_{13}\\\ \star&{\lambda}_{22}&\bm{0}_{M\times N}&{\lambda}_{23}\\\
\star&\star&\bm{0}_{N\times N}&\bm{0}_{N}\\\
\star&\star&\star&{\lambda}_{33}\end{bmatrix},$ (23a)
$\displaystyle\Lambda_{MPC}=\begin{bmatrix}\bm{0}_{n_{x}\times
n_{x}}&\bm{0}_{n_{x}\times M}&-(\tau_{QP}h^{\top})/2&\bm{0}_{n_{x}}\\\
\star&\bm{0}_{M\times M}&\bm{0}_{M\times N}&\bm{0}_{M}\\\
\star&\star&-\tau_{QP}H&-L^{\top}\tau_{in}\\\
\star&\star&\star&b^{\top}\tau_{in}+\tau_{in}^{\top}b\end{bmatrix},$ (23b)
$\displaystyle X=\begin{bmatrix}-2\tau_{x}&\bm{0}_{n_{x}\times
M}&\bm{0}_{n_{x}\times N}&\tau_{x}(\overline{x}+\underline{x})\\\
\star&\bm{0}_{M\times M}&\bm{0}_{M\times N}&\bm{0}_{M}\\\
\star&\star&\bm{0}_{N\times N}&\bm{0}_{N}\\\
\star&\star&\star&-(\overline{x}^{\top}\tau_{x}\underline{x}+\underline{x}^{\top}\tau_{x}\overline{x})\end{bmatrix},~{}\check{\Lambda}_{\text{MPC}}=\begin{bmatrix}\bm{0}_{n_{x}\times
n_{x}}&\bm{0}_{n_{x}\times M}&-h^{\top}/2&\bm{0}_{n_{x}}\\\
\star&\bm{0}_{M\times M}&\bm{0}_{M\times N}&\bm{0}_{M}\\\
\star&\star&-H&-L^{\top}\\\ \star&\star&\star&b+b^{\top}\end{bmatrix}.$ (23c)

Ross Drummond would like to thank the Royal academy of engineering for funding
this research through a UKIC Fellowship. This work was partly funded by
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under
Germany’s Excellence Strategy - EXC 2075 - 390740016. We acknowledge the
support by the Stuttgart Center for Simulation Science (SimTech). The authors
thank the International Max Planck Research School for Intelligent Systems
(IMPRS-IS) for supporting Patricia Pauli.

## References

  * Alessio and Bemporad (2009) Alessandro Alessio and Alberto Bemporad.  A survey on explicit model predictive control.  In _Nonlinear model predictive control_ , pages 345–369. Springer, 2009.
  * Andersen and Andersen (2000) Erling D Andersen and Knud D Andersen.  The MOSEK interior point optimizer for linear programming: an implementation of the homogeneous algorithm.  In _High Performance Optimization_ , pages 197–232. Springer, 2000\.
  * Angeli (2009) David Angeli.  Convergence in networks with counterclockwise neural dynamics.  _IEEE Transactions on Neural Networks_ , 20(5):794–804, 2009.
  * Bai et al. (2019) Shaojie Bai, J Zico Kolter, and Vladlen Koltun.  Deep equilibrium models.  _arXiv preprint arXiv:1909.01377_ , 2019.
  * Barabanov and Prokhorov (2002) Nikita E Barabanov and Danil V Prokhorov.  Stability analysis of discrete-time recurrent neural networks.  _IEEE Transactions on Neural Networks_ , 13(2):292–303, 2002.
  * Bemporad et al. (2002) Alberto Bemporad, Manfred Morari, Vivek Dua, and Efstratios N Pistikopoulos.  The explicit linear quadratic regulator for constrained systems.  _Automatica_ , 38(1):3–20, 2002.
  * Chu and Glover (1999) Yun-Chung Chu and Keith Glover.  Bounds of the induced norm and model reduction errors for systems with repeated scalar nonlinearities.  _IEEE Transactions on Automatic Control_ , 44(3):471–483, 1999.
  * Darup (2020) Moritz Schulze Darup.  Exact representation of piecewise affine functions via neural networks.  In _Procs. of the European Control Conference_ , pages 1073–1078. IEEE, 2020.
  * Desoer and Vidyasagar (2009) Charles A Desoer and Mathukumalli Vidyasagar.  _Feedback systems: Input-output properties_.  SIAM, 2009.
  * Drummond et al. (2021) Ross Drummond, Mathew C Turner, and Stephen R Duncan.  Reduced-order neural network synthesis with robustness guarantees.  _arXiv preprint arXiv:2102.09284_ , 2021.
  * El Ghaoui et al. (2021) Laurent El Ghaoui, Fangda Gu, Bertrand Travacca, Armin Askari, and Alicia Tsai.  Implicit deep learning.  _SIAM Journal on Mathematics of Data Science_ , 3(3):930–958, 2021.
  * Fabiani and Goulart (2021) Filippo Fabiani and Paul J Goulart.  Reliably-stabilizing piecewise-affine neural network controllers.  _arXiv preprint arXiv:2111.07183_ , 2021.
  * Fazlyab et al. (2019) Mahyar Fazlyab, Alexander Robey, Hamed Hassani, Manfred Morari, and George J Pappas.  Efficient and accurate estimation of Lipschitz constants for deep neural networks.  _arXiv preprint arXiv:1906.04893_ , 2019.
  * Fazlyab et al. (2020) Mahyar Fazlyab, Manfred Morari, and George J Pappas.  Safety verification and robustness analysis of neural networks via quadratic constraints and semidefinite programming.  _IEEE Transactions on Automatic Control_ , 2020.
  * Green and Limebeer (2012) Michael Green and David JN Limebeer.  _Linear robust control_.  Courier Corporation, 2012.
  * Hertneck et al. (2018) Michael Hertneck, Johannes Köhler, Sebastian Trimpe, and Frank Allgöwer.  Learning an approximate model predictive controller with guarantees.  _IEEE Control Systems Letters_ , 2(3):543–548, 2018.
  * Jönsson (2001) Ulf Jönsson.  Lecture notes on integral quadratic constraints.  2001\.
  * Li et al. (2006) Guang Li, William P Heath, and Barry Lennox.  The stability analysis of systems with nonlinear feedback expressed by a quadratic program.  In _Procs. of the Conference on Decision and Control_ , pages 4247–4252. IEEE, 2006.
  * Li et al. (2021) Jiaqi Li, Ross Drummond, and Stephen R Duncan.  Robust error bounds for quantised and pruned neural networks.  In _Learning for Dynamics and Control_ , pages 361–372. PMLR, 2021\.
  * Lofberg (2004) Johan Lofberg.  YALMIP: A toolbox for modeling and optimization in MATLAB.  In _International Conference on Robotics and Automation_ , pages 284–289. IEEE, 2004.
  * Parisini and Zoppoli (1995) Thomas Parisini and Riccardo Zoppoli.  A receding-horizon regulator for nonlinear systems and a neural approximation.  _Automatica_ , 31(10):1443–1451, 1995.
  * Pauli et al. (2021a) Patricia Pauli, Dennis Gramlich, Julian Berberich, and Frank Allgöwer.  Linear systems with neural network nonlinearities: Improved stability analysis via acausal Zames-Falb multipliers.  _arXiv preprint arXiv:2103.17106_ , 2021a.
  * Pauli et al. (2021b) Patricia Pauli, Anne Koch, Julian Berberich, Paul Kohler, and Frank Allgower.  Training robust neural networks using Lipschitz bounds.  _IEEE Control Systems Letters_ , 2021b.
  * Revay et al. (2020) Max Revay, Ruigang Wang, and Ian R Manchester.  Lipschitz bounded equilibrium networks.  _arXiv preprint arXiv:2010.01732_ , 2020.
  * Šiljak and Zečević (2005) Dragoslav D Šiljak and AI Zečević.  Control of large-scale systems: Beyond decentralized feedback.  _Annual Reviews in Control_ , 29(2):169–179, 2005\.
  * Silver et al. (2017) David Silver, Julian Schrittwieser, Karen Simonyan, Ioannis Antonoglou, Aja Huang, Arthur Guez, Thomas Hubert, Lucas Baker, Matthew Lai, Adrian Bolton, et al.  Mastering the game of Go without human knowledge.  _Nature_ , 550(7676):354–359, 2017.
  * Sutton and Barto (2018) Richard S Sutton and Andrew G Barto.  _Reinforcement learning: An introduction_.  MIT press, 2018.
  * Szegedy et al. (2013) Christian Szegedy, Wojciech Zaremba, Ilya Sutskever, Joan Bruna, Dumitru Erhan, Ian Goodfellow, and Rob Fergus.  Intriguing properties of neural networks.  _arXiv preprint arXiv:1312.6199_ , 2013.
  * Turner and Drummond (2021) Matthew C Turner and Ross Drummond.  Discrete-time systems with slope restricted nonlinearities: Zames-Falb multiplier analysis using external positivity.  _International Journal of Robust and Nonlinear Control_ , 31(6):2255–2273, 2021.
  * Yin et al. (2021a) He Yin, Peter Seiler, and Murat Arcak.  Stability analysis using quadratic constraints for systems with neural network controllers.  _IEEE Transactions on Automatic Control_ , 2021a.
  * Yin et al. (2021b) He Yin, Peter Seiler, Ming Jin, and Murat Arcak.  Imitation learning with stability and safety guarantees.  _IEEE Control Systems Letters_ , 2021b.
  * Zames and Falb (1968) George Zames and PL Falb.  Stability conditions for systems with monotone and slope-restricted nonlinearities.  _SIAM Journal on Control_ , 6(1):89–108, 1968\.
  * Zhang et al. (2019) Xiaojing Zhang, Monimoy Bujarbaruah, and Francesco Borrelli.  Safe and near-optimal policy learning for model predictive control using primal-dual neural networks.  In _Procs. of the American Control Conference_ , pages 354–359. IEEE, 2019.