\section*{Results}

In section \ref{public}, we will analyze  the diversity of \CDR sequences as well as the fraction of public sequences in the TCR repertoire.
In section \ref{network_structure}, we will investigate network properties of the TCR repertoire, namely the average neighbor degree distribution and  the robustness of the system against
removal of T cells. In section \ref{virus}, we will finally investigate  the network properties of antigen-cognate \CDR sequences and determine the a.a. composition of the \CDR sequences.
\subsection{Public sequences}
\label{public}
During V(D)J recombination each unique  \CDR sequence has a certain probability of being generated.

Thus, each individual has a probability $p_{ind}(\boldsymbol{\sigma})$ of possessing a TCR with a unique \CDR sequence  $\boldsymbol{\sigma}$, which can be high or low.
Hence, we do not just want to know whether a sequence occurs in
only one individual (private) or in every individual (public).
Instead, the possibility arises that the sequence exists in multiple, but not
all individuals. That prompts the question whether these sequences,
which occur in multiple but not all individuals,
are private or public. This topic was investigated in
previous publications~\cite{elhanati2018predicting, quigley2010convergent, robins2009comprehensive, freeman2009profiling, benichou2012rep, six2013past, robins2010overlap, venturi2011mechanism, elhanati2014quantifying, zvyagin2014distinctive,  pogorelyy2017persisting, madi2014t, soto2020high}, still
 there are different definitions of public, ranging from
 defining any sequence that occurs twice or more~\cite{miho2019large}
 to restricting the public category to sequences which occur in
 all individuals~\cite{madi2017t}.
In order to obtain a full view, Fig. \ref{fig:fig_2} (a) shows
the frequency distribution of sequences, as explained below.
\begin{figure}
    \includegraphics[scale = 1.4]{ims/public_scheme.pdf}
    \caption{{\bf Conceptual illustration of the sharing number analysis for $n=4$.} Each human
    $i$ possesses a unique TCR repertoire with a set $T_i$ of N \CDR  sequences, labeled by $\sigma_{ij}$ ($j^{th}$ sequence of the $i^{th}$ human). The set $A = \{T_1, ..., T_4\}$ of sets contains all TCR repertoires.}
    \label{fig:public_scheme}
\end{figure}

\subsubsection{Frequency decreases with sharing number}
May the set of $n$ TCR repertoires, each from a unique individual,  be represented as $A = \{ T_i~|~0<i\leq n \}$ (see Fig. \ref{fig:public_scheme}), where $T_i = \{ {\boldsymbol \sigma_{ij} ~|~ 0<j\leq N}\}$ is a TCR repertoire comprising $N$ sequences  and we assume that all TCR repertoires $T_i$ have the same size $N$, where {\it size} of a repertoire is defined as the number of \CDR sequences within the repertoire.

If for a given  sequence $\boldsymbol \sigma_{ij}$  another sequence $\boldsymbol \sigma_{i'j'}$ exists with $\boldsymbol \sigma_{ij} = \boldsymbol \sigma_{i'j'}$ and $i \neq i'$, that means
that this sequence occurs in two different individuals.
Hence, according to the
definition in Ref.~\cite{miho2019large}, it is public. Let us define the {\it sharing number}
$n_{{\boldsymbol \sigma}}$ as the number
of individuals in which a sequence equal to ${\boldsymbol \sigma}$ can be found~\cite{elhanati2018predicting}.
If $n_{{\boldsymbol \sigma}} > 1$, the sequence is public, otherwise private. Moreover,
let us define $M_k$ as the number of different sequences $\sigma_{ij} \in A$ \footnote{strictly speaking $\sigma_{ij}\in \bigcup A := \{{\boldsymbol  \sigma}|\exists T \in A: {\boldsymbol \sigma} \in T\}$}
with
$n_{{\boldsymbol \sigma_{ij}}} \geq k$, where $k$ may be any non-negative integer.
In panel (a) of Fig. \ref{fig:fig_2} the normalized frequency $f_N(k) := \frac{M_k}{M_1}$ is plotted in dependence of $k$ on a logarithmic scale for $n=4$
TCR repertoires, each comprising $N=10^6$ \CDR sequences, generated by SONIA (green for human, blue for mouse) and the same number of
empirical human \CDR sequences.
First of all, we observe that for humans $f_N$ declines sharply with
increasing $k$.
A decline can be expected since it is less probable to find a sequence in many samples than in just a few.
The empirical data exhibits a
slower decline, though the qualitative  behavior is similar to that of generated
sequences, and is in qualitative agreement with the results of a previous  study~\cite{elhanati2018predicting} with more human subjects but fewer \CDR sequences per repertoire.

\subsubsection{Sharing number depends on repertoire size}

In Fig. \ref{fig:fig_2} (b) we show the distribution $f_N(k)$ from Fig. \ref{fig:fig_2} (a) again, this time for different $N \in \{1/4\times 10^{5}, 1/2\times 10^{5}, 1\times 10^{5}, 2\times 10^{5}, 4\times 10^{5}\}$.
The distribution shifts upwards for higher $N$. The inset of Fig. \ref{fig:fig_2} (b) plots $f_{N}(k=2)$ in dependence of $N$ with logarithmic $N$-axis.
For comparison, we show a straight line.

For empirical and generated sequences we  observe a logarithmic increase of $f_N(k=2)$ with $N$ and therefore can confirm the predictions of Ref. ~\cite{elhanati2018predicting} up to  $N=400,000$ sequences. The increase is expected to slow down and saturate for higher $N$.

In agreement with previous studies~\cite{elhanati2018predicting}, it can be concluded that the fraction of public frequencies in a TCR repertoire also depends on the size $N$ of the repertoire itself.
\subsubsection{Selection effects}
 V(D)J recombination spawns every unique \CDR sequence ${\boldsymbol \sigma}$ with a designated probability $p$.
 Generating a set
of sequences by V(D)J recombination and subsequent thymus selection is thus a weighted combination $C_{{\boldsymbol p}}(N_{max}, N)$ with repetitions (for more thorough definition see sections \ref{Combination} and \ref{Combination of amino acids} under Methods),
where $N_{max}$ denotes the number of theoretically possible different sequences (sequence types), which is estimated to
be up to more than $10^{15}$~\cite{miles2011bias, laydon2015estimating},
though other estimates suggest that the real number of possibilities lies only around  $10^{6}$~\cite{robins2009comprehensive}.
$N$ is again the number of generated sequences for a TCR repertoire and ${\boldsymbol p} = \{ p_i~|~ 0<i\leq N_{max}\}$ refers to the set of probabilities $p_i$ for generating a certain sequence type ${\boldsymbol \sigma_i}$.
 This can be compared to an unweighted combination $C_{\mathbb{1}}(N_{max}, N)$,
 where all probabilities $p_i$ are equal.
After a combination each sequence type $\boldsymbol{\sigma}$ occurs with a frequency $f(\boldsymbol{\sigma})$. Let's denote the number of different sequence  types $\boldsymbol{\sigma}$ that have the same frequency $f(\boldsymbol{\sigma})$ with $n(f)$.
Fig.~\ref{fig:fig_2}~(c)
plots $n(f)$ over $f$ for a TCR repertoire of size $N=10^{6}$ (blue and green dots for SONIA generated sequences of mouse and human respectively, green diamond markers for empirical human sequences).
The generation of this
repertoire can be viewed as a weighted combination $C_{\boldsymbol p}(N_{max}, N=10^{6})$, where $N_{max}$
can range from $10^{6}$ to more than $10^{15}$.
For comparison, $N=10^{6}$
numbers were randomly chosen from the range of the first $N_{max}=10^{6}$ natural numbers, i.e. from the set $\{i~|~ 0<i\leq 10^6\}$,
 in an unweighted combination  $C_{\mathbb{1}}(10^{6}, 10^{6})$.
 The resulting function $n_{\mathbb{1}}(f)$ (red curve)
 declines significantly faster than that of the weighted
 combination $C_{\boldsymbol p}(N_{max}, N)$. Hence,
 $n_{\mathbb{1}}$ has more \enquote{types} with a low
 frequency but less  with a high one, compared to the
 distribution $n(f)$ of real sequence  types. This confirms that
 the distribution of the TCR repertoire does not just result
 from random selection but instead from enhancement and
 suppression during V(D)J  recombination and thymus selection.
Some \CDR sequences are chosen with a probability above the random average, thus they are enhanced.
Others stay below the average, they are suppressed. We mark the crossing point in red, where \CDR sequences are neither suppressed nor enhanced.
The empirical data (diamond shaped green marks) show a similar decline, compared to the generated \CDR sequences.
Of course, this distinction remains only schematic, since
the exact value of $N_{max}$ in the TCR repertoire is not known exactly and estimates deviate by several orders of magnitude. But it clearly shows that the histogram mainly reflects the selection effects from V(D)J recombination and Thymus selection.

\subsection{Network structure of TCR repertoires}
\label{network_structure}
The ample size and complexity of the TCR repertoire impedes the understanding of its structure and  biological implications.  The decades-old field of network
science provides a set of tools to overcome this hurdle. Hence, in the following section, we focus on basic
network properties, namely the degree distribution, and the LCC of a TCR network. We employed our TCR network construction framework TCR-NET to create large networks, which form the base for this analysis. Nevertheless, also the
calculation of network properties for a given network can require significant computational effort,
and varies greatly between different network properties. Therefore, we will
use different network sizes in dependence of the network properties that we want to determine, varying from several thousand \CDR sequences in computationally expensive properties such as the betweenness distribution or the assortativity to up to $8\times 10^5$ \CDR sequences for the robustness analysis. Doing that enables us to always use the largest networks possible for each analysis.

\subsubsection{Assortativity and mean degree}
In Fig. \ref{fig:fig_2} (d) we analyzed the
assortativity and the mean degree of the TCR network.
We repeated the analysis of Fig. \ref{fig:fig_2}  (a) for $n=4$ TCR repertoires, but this time each of size $N = 8\times 10^{3}$
and selected all \CDR sequences (in total $N_{total}=32\times 10^3$)
to construct the \enquote{total network} with $l_{max}=2$.
Then we selected only the sequences,
labeled as \enquote{public} and constructed the
network with them, labelling it the \enquote{public network}
as illustrated in Fig. \ref{fig:fig_2} (d), choosing $l_{max} = 2$ and having $N_{sim} = 364$ and $N_{exp} = 646$ empirical and generated  public sequences respectively. For generated sequences as well as experimentally obtained sequences, the public network has a similar or even lower assortativity but a higher
mean degree than the total network, reflecting
the denser structure. This indicates, that public
sequences are more similar to each other, while
private sequences cover a more diverse range of sequences for immune protection.
\subsubsection{Betweenness centrality and neighbor degree distribution}
\label{neighbor}
Panel (e) of Fig. \ref{fig:fig_2} shows the normalized histogram of the betweenness centralities of nodes for a network constructed from $N=8\times 10^3$ sequences and with $l_{max} = 1$.
There are a relatively small number of nodes with a
high centrality compared to the majority of nodes with much lower centrality.
In panel (f) of Fig. \ref{fig:fig_2} the average neighbor degree $\langle d_{nn}(d) \rangle$
(see section \ref{complex_networks}) is plotted
in dependence of the degree $d$ of a node itself, with \NSeqPublicAnalysis
and $l_{max} = 1$.
We observe that a higher degree is associated with a higher average neighbor degree.
The latter increases proportionally with the square root
of the degree of the node itself (see  solid curves in panel (f) of Fig. \ref{fig:fig_2}) for generated,
as well as empirical sequences.

From the positive slope and the positive assortativity (for the latter see Fig. \ref{fig:fig_2} (d)) we conclude that the TCR network is an  assortative network.
Thus, there are dense areas, where nodes are highly interconnected. On
the other hand, the nodes with low degree are connected to neighbors
which also have a low degree.
This suggests that the immune system prepares against a certain group  of  epitopes by creating a dense
net of interconnected TCRs while against more \enquote{unusual} antigens there is
only a relatively weak protection.
  To visualize that, the inset
 shows a TCR network of $N = 8\times 10^{3}$ \CDR sequences with $l_{max} = 1$, where the
 radial position of each node is inversely proportional
 to its degree, while the polar angles
 are chosen randomly. Each node is colored depending on its average neighbor degree, yellow indicating a high value, while dark blue a low value. It can be observed, that the nodes in the center, i.e. having a high degree, are also connected to higher degree nodes, i.e. they have a higher average neighbor degree.

A goal for the future would be to investigate, how the immune system
responds to the aforementioned \enquote{unusual antigens}. Perhaps they are really so
scarce that the necessity to prepare an immune  response does not arise.
It also remains unclear if TCRs in these empty areas
are perhaps more versatile and able to
respond to a more diverse range of antigens, differing among each other by a higher Levenshtein distance (higher degeneracy).
The observed assortativity is unusual, since many biological networks tend to be more disassortative~\cite{newman2002assortative}.
Nevertheless, positive degree assortativity
is also found in other networks outside of the area of immunology, e.g.
the \enquote{e-printarchive coauthorship network} from arxiv.org~ \cite{lee2006statistical}. They all form the aforementioned group of assortative networks.

\subsubsection{Robustness of TCR repertoires}
Spontaneous mutations in the TCR repertoire of humans and animals can alter their immune systems'
response to antigens.
Moreover, the average lifetime of T cells is estimated to be not more than
a few weeks~\cite{westera2013closing}.  Therefore, the structure of the TCR repertoire changes constantly. Moreover,  diseases such as HIV can lead to a  reduction of T cells~\cite{philip2021cd8, hazenberg2000t}. This raises the question, how long characteristic structural properties stay preserved.
We simulate a change in the repertoire with the removal of a TCR, thus the removal
of a node, which can be caused, e.g., by apoptosis, lysis, or anergy
of T cells.
To study the robustness of the TCR network, we then analyze how  certain TCR network properties change under node removal, which  is also commonly performed for other network
types in graph theory~\cite{callaway2000network, schneider2013algorithm} and can be understood as a site percolation problem (see section \ref{transition}).

In Fig. \ref{fig:fig_2} (g) we show the normalized relative size (for definition see section \ref{complex_networks} and \ref{computing_LCC} under Methods) of the largest connected component (LCC) as a function of $1-p$, where $p := \frac{N_{r}}{N}$ denotes the fraction of removed nodes, $N_{r}$ the total number of removed nodes, and $N$ the total number of nodes in the network before any node removal. We chose the parameters  \NSeqRobustness and $l_{max} = 1$.
The relative size of the LCC continuously  declines during node removal. But even after removing half of all nodes ($p=0.5$), the relative size has decreased by less than twenty percent, indicating
a strong robustness of the TCR network and low susceptibility to damages. The main cluster stays preserved. This  is in agreement with the results found
in Ref.~\cite{miho2019large} for antibodies
and necessary in order to maintain the
immunological footprint of an individual.
The empirical
data show a similar robustness of the immune system, with only slight deviations.
On the other hand, a loss of more than eighty percent of all nodes strongly damages the network. In contrast to other network types, such as Erd{\H o}s--R{\'e}nyi graphs~\cite{tishby2018revealing},
 the TCR network does not exhibit a phasetransition at some point $p_C<1$ where the LCC abruptly vanishes.
Instead, a significant portion
always remains, which continuously decreases toward zero when $1-p$ becomes close to
zero, similar to the Barabási--Albert model~\cite{PhysRevLett.85.4626} (see section \ref{complex_networks} under Methods).
This corresponds to the limes
\begin{equation}
    \lim\limits_{p \rightarrow 1} \textnormal{LCC} = 0
\end{equation}
The spontaneous jumps in the relative size of the LCC can be explained with a large connected component, only being connected to the rest of the LCC by one remaining edge. If this last edge is removed due to the removal of a node, then the size of the LCC can shrink by a large portion, since a large number of nodes is suddenly cut of.
In Fig. \ref{fig:fig_2} (h) we plot the normalized degree distributions (see section \ref{complex_networks} and \ref{normalization})
corresponding to the node-removal probabilities
$p\in \{ 0.01, 0.49, 0.89, 0.95, 0.98\}$
for the same network parameters as in Fig. \ref{fig:fig_2} (g). A high $p$ is represented by a bright green line, a low $p$ by a dark green tone for generated sequences and analoguously in blue for empirical sequences.  The normalized frequency of each degree value is plotted in dependence of the respective degree as a histogram.
As expected, for empirical as well as generated sequences, the frequency decreases with higher $p$, i.e. higher number of removed nodes, but the
negative slope also increases. The networks degree distribution shows similarity with an algebraic decline, which
is also found in several  other
graph types (see also section \ref{complex_networks}).

\subsubsection{LCC transition}
A TCR recognises an antigen, e.g. from a virus or a cancer cell,
 by detecting a characteristic site of the antigen, called {\itshape epitope},
  which matches the {\itshape paratope} of the TCR.
Usually, the paratope is not specific to only one epitope of an antigen (for a detailed introduction to paratopes and epitopes see Ref.~\cite{e2003fundamental} chapter 4). During the early 2000s
 it became concensus that there is a significant degeneracy
 where each paratope can recognize multiple epitopes, then called a {\itshape mimotopic array}~\cite{cohn2005degeneracy},
 and each epitope can be recognized by multiple paratopes, referred to
 as a {\itshape paratopic clan}~\cite{cohn2005degeneracy}.
 Thus, the TCR repertoire does not need every paratope, i.e., every \CDR sequence, to be present because a paratope can recognize some neighbor epitopes too, as long as they remain sufficiently similar to the \enquote{ideal} epitope. We model the condition of sufficient similarity by requiring the Levenshtein distance between \CDR sequences to be smaller than or equal to a maximum distance $l_{max}$, as explained in the introduction.
 We want to investigate
 how small $l_{max}$ can be to guarantee full coverage,
 which we define as given if the size of the LCC is of the same order of magnitude as the entire TCR repertoire ($LCC \sim N$).
 It should be noted that the Levenshtein distance is a very rough
 measure of the similarity between two \CDR sequences,
 which actually depends on molecular interactions. In some cases a small Levenshtein distance can lead to nonrecognition while in other cases an epitope can be recognized by a  sequence despite a high  Levenshtein distance to the ideal cognate sequence~\cite{wu2002two,aleksic2012different,reiser2003cdr3,wucherpfennig1995molecular,wucherpfennig2007polyspecificity}. Thus, the investigation using the Levenshtein distance
  should be always treated as an approximation.

  We will study the relative size of the LCC, as we did during the robustness analysis (see Fig. \ref{fig:fig_2} (g)), this time in dependence of the threshold value $l_{max}$, as shown in  Fig.  \ref{fig:fig_2} (i).
The relative size
 of the LCC
is plotted in dependence of $l_{max}$ for different parameters
$N\in  \{1\times 10^{3}, 4\times 10^{3},  8\times 10^{3}\}$.
While at $l_{max} = 0$ the LCC covers a minority of nodes,
a transition can be observed for increasing $l_{max}$,
until the LCC covers the entire network and full coverage is given according to the previous definition.
We can observe, that for higher $N$, this transition shifts left, i.e., toward lower thresholds $l_{max}$.
We can interpret the position of this transition as an estimate of how big the disparity between a given \CDR sequence and an epitope's \enquote{ideal} cognate paratopes \CDR sequence can be in order
to still recognize the antigen and therefore enable an immune response.
As larger the total number $N$ of TCRs, as
smaller the requisite tolerance of each single node, thus as smaller the
required degeneracy of each TCR for
provoking a sufficient immune response.
On the other hand, this indicates that smaller organisms may
 need
a wider responsiveness of each TCR to antigens due to their smaller TCR repertoire.
This might be achieved through higher degeneracy.
Alternatively, it might also be possible that smaller
organisms need to cover a less diverse sequence pool,
due to possible lower diversity of cell genes and
possibilities of antigens to react with healthy cells.

\subsection{Distinguishing antigens with network analysis}
\label{virus}
\begin{figure*}[htb]
    \centering
    \includegraphics[scale = 0.7]{ims/virus_colored_3.pdf}
    \caption{{\bf Network topology of different antigens}. The scheme in the upper left illustrates the insertion of antigen cognate \CDR sequences in the TCR network (explained in the text).
Panel (a) shows the degree distribution of
    different antigens ($N=4\times 10^5,  l_{max}=1$). Though they all resemble an algebraic decline, their decay exponents $b$ differ and are shown in panel (b).
Panel (c) shows the monotonically positive
    average neighbor degree
    distribution for the same network.
Panel (d) and (f) show the radarplots of
    assortativity and clustering respectively.
    }
    \label{fig:fig_3}
\end{figure*}

 The paratope is the part of the TCR, that comes into contact with the antigen epitope (presented by an {\itshape antigen presenting cell})
  and
  mainly consists
  of the three CDR regions within the variable domain of the TCR (for a details we refer to chapter 4 in Ref.~\cite{e2003fundamental}.
  Since the diversity within a single individual predominantly arises from the CDR3 region, it suffices to represent a paratope by its CDR3 region, that can then be matched to the specific epitope, called the {\itshape cognate} epitope.
  This of course neglects cross reactivity  (for introduction in cross reactivity see Ref.~\cite{e2003fundamental})
   and hence should be treated as an approximation.
  In a previous work a large database was created~\cite{tickotsky2017mcpas}, matching
  \CDR chains to the corresponding antigen epitopes for different antigen families.
  We use this data in order to represent antigens through their corresponding cognate \CDR sequence
  and place them as a node in our TCR network.
  To do that, we begin with a normal list
  of $N$ \CDR empirical sequences $\{\sigma_i|\ 0 < i \leq N\}=:\{\sigma_i\}$. Moreover, we define the list $\{\sigma_{a}|\ 0 < i \leq N_a, N_a \ll N\}=:\{\sigma_a\}$ as the list of $N_a$ \CDR sequences, that are cognate to at least one antigen from the antigen family $a$ (e.g. Influenza or Tuberculosis).
  We choose the antigen families of the viruses Influenza, HIV, Tuberculosis
  as well as of Alzheimer and Parkinson, where, though not infectious diseases, the immune system plays a major role ~\cite{heppner2015immune, jevtic2017role,  kustrimovic2018parkinson}.
  We now select all $N_a$ \CDR sequences from the
  set $\{\sigma_{a}\}$ and replace the first $N_a$ sequences in $\{\sigma_i\}$ with them. We call the resulting set $\{\sigma_{i,a}\}$ and construct the network from it,  labeling it an {\it antigen specific network} or briefly {\it ASN}.
  In this way we create one network for each antigen and
  can compare the network properties of  antigen cognate nodes, inserted from $\{\sigma_a\}$ (abbreviated with {\itshape ac nodes}).
The network parameters are chosen to be $l_{max} = 1$,
 $N=4\times 10^5$. Moreover, the values $N_A$ are 6051, 1291, 204, 177, 16663 for Influenza, HIV, Alzheimer, Parkinson, Tuberculosis respectively and therefore of a lower order of magnitude than $N$.

Panel (a) of Fig. \ref{fig:fig_3} shows the normalized degree distribution
of ac nodes in different ASNs. Though they share the common behavior of an algebraic decline, as in many other
network types,
their decline exponents differ.
While Alzheimer slowest decline with a
decay exponent of $b = -0.419\pm 0.018$,
the degree distribution of Tuberculosis declines fastest with the exponent of $b = -1.563\pm 0.005$, as shown in Panel (b).  In contrast to the
other investigated antigens, Influenza
has an unusually large number of nodes with high degree,
deviating from the algebraic curve. The given standard deviations should be taken with caution, since the fit relies on the assumption of an algebraic decline. If the degree distribution differs slightly from a perfect algebraic distribution, this would lead to additional deviations.
In panel (c) we show the average neighbor degree distribution, similar to panel (f) of Fig. \ref{fig:fig_2}, but only took ac  nodes into account.
Again, we observe a square root curve, even with similar scaling throughout the different antigens, within
the bounds of uncertainty. Therefore, also the subgraphs of ac nodes are assortative,  if we include edges between ac nodes and other nodes in the network.
Panel (d) and (e) show the assortativity and the clustering coefficient
for different ASNs, only calculated from ac nodes. They deviate significantly throughout the antigens, making them a possible marker to distinguish antigens.

It should be noted that a single antigen often carries
multiple epitopes
and thus can react to multiple TCRs.
More databanks about cross reactivity would offer
 opportunities to perform further network analysis in the future.

\begin{figure}
\centering
\includegraphics[scale = 1.0]{ims/amino.pdf}
\caption{{\bfseries Histogram of the a.a. frequencies.}
We show the distribution of a.a. relative frequencies for a human TCR repertoire of size $N = 10^{6}$ before the thymus selection (\enquote{pre}) and
after thymus selection (\enquote{post}), using SONIA generated \CDR sequences. While the distributions don't differ largely from each other, they both deviate significantly from the a.a. distribution in human plasma (\enquote{HP}). The amino acids are denoted by their standard one-letter abbreviations.
}
\label{fig:aminoHist}
\end{figure}

\section*{Amino acid composition}
TCR receptors consist of amino acid sequences (a.a.s.).
That leads to the question how the a.a.s. are  distributed.
This has been analyzed for various
different tissues over many decades.
Deviations in the a.a. composition can indicate
illnesses~\cite{mally1997changes}.
Ref.~\cite{song2018recent}
reviews the currently applied methods and different tissues that have been
analyzed, including human plasma.
We compare the a.a. composition of the \CDR
sequences with the composition of human plasma~\cite{song2018recent, mcmenamy1960studies} in Fig. \ref{fig:aminoHist}.
It shows the normalized frequency of each a.a. within a sample of $N = 10^6$
 \CDR sequences from TCR$\beta$ chains, generated with SONIA. The distribution before thymus
selection (light blue)
is
similar to the distribution after the thymus
selection (blue), with exception of K, P, R and S (we refer to amino acids, using their 1-letter symbols).
These sequences are significantly enhanced (S) or inhibited (K, P, R)
by the thymus.
Both distributions show much stronger differences, compared to
the average a.a. distribution in human plasma (red),
with exceptions of
A, E, N, P, R and T, where the relative concentrations
  only differ by a relatively small amount.
 It can be concluded that the \CDR sequences do not just resemble
 the average a.a. composition of human plasma but instead have a unique composition, which is
 only changed slightly during thymus selection.

\section{Discussion}
 \label{conclusion}
 We created the package {\scshape TCR-NET}, which improves the performance compared to previous GPU based algorithms for Levenshtein distance calculation~\cite{balhaf2016using} by exploiting the specificities of \CDR sequences in the TCR repertoire.
 Using the SONIA sequence simulation software package~\cite{sethna2020population}, we then
 created a sequence pool in order to analyze the fraction of public and private
 sequences. In coherence with previous analysis~\cite{elhanati2018predicting, quigley2010convergent, robins2009comprehensive, freeman2009profiling, benichou2012rep, six2013past, robins2010overlap, venturi2011mechanism, elhanati2014quantifying, zvyagin2014distinctive,  pogorelyy2017persisting, madi2014t}, it was found that the fraction of public sequences depends on the
 definition of public, which still differs among different scientific publications.
 Moreover, a higher repertoire size increases the share of public  sequences.
  Experimental evidence suggests that the
 final fraction of public sequences lies significantly above the theoretically predicted values resulting from pure stochasting modeling of V(D)J recombination~\cite{soto2020high}.
 We could confirm a difference between generated \CDR sequences and empirical sequences from naive CD8 T cells throughout all investigated definitions of public.
 Moreover, a comparison with a random distribution showed that the public sequences are not just the result of a stochastic overlap,
 which could also occur in an entirely random distribution, but instead must derive from enhancement and suppression effects during V(D)J recombination and thymus selection.

 In the next step, we investigated the network structure of the TCR repertoire.
 It was found that the average neighbor degree of a node increases proportionally to the square root of the degree of the node itself,  showing that the network is assortative. Nodes with many neighbors seem to be connected with each other in a
 very dense center while some isolated nodes can be found far outside.
 The LCC for TCR networks from SONIA-generated \CDR sequences as well as empirical sequences proved to be robust against large damages ($N_r \sim N$), indicating that the TCR repertoire can retain its charactaristic topology.
 Of course the network properties depend significantly on  $l_{max}$. If $l_{max}$ increases for constant $N$, the size of the LCC increases as well until
  covering the entire network. The curve expresses a transition
  behavior, where the transition shifts towards smaller values of $l_{max}$ for larger
  $N$.

  Moreover, we inserted \CDR sequences that correspond to certain antigens~\cite{tickotsky2017mcpas} into TCR networks and analyzed their network properties.
  Especially Influenza exhibited an unusually high number of high degree nodes, deviating from the algebraic decline observed for the degree distributions of all investigated antigens.
  Moreover, the clustering coefficients and assortativities differ between the  antigens.
  The neighbor degree distributions for all antigen specific networks increase monotonically, indicating an assortative network.
  All in all this provides strong evidence, that antigen families can be distinguished by their network topology, especially the degree distribution,  clustering coefficient and assortativity.

  Finally, we investigated the a.a. composition of generated human \CDR sequences  and found that it  deviates significantly from the average a.a. composition in human plasma. Moreover, there are also (though more slight) differences between the
  compositions before and after thymus selection. This might help distinguishing samples and opens up perspectives to empirical  investigations on how
  diseases or special conditions might change the a.a. composition in the TCR repertoire, as it
  was done for other parts of the human body~\cite{engelen2001effects,mally1997changes}.

  In conclusion, we studied the share of public \CDR sequences, the robustness of the immune system, the network structure, the relation to antigen  epitopes as well as the a.a. composition and related it to previous work.
  Especially the application of network theory in TCR repertoire analysis is
  a promising field, which still opens up huge spaces for further investigations and can help understanding the immune systems structure
  and responsiveness in more detail than it was previously possible.
  Especially with the rise of personalized medicine~\cite{hamburg2010path} it will become necessary
  to understand the specific differences in the immune system between individuals. This will require a very detailed analysis of the complex topology of the TCR repertoire. Thus, especially with the rise of high-throughput sequencing and consequently the availability of more and
  more data, the field of network theory is an increasingly helpful tool to
  understand the TCR repertoire.

\acknowledgements
We would like to thank Prof. Dr. Markus Sigrist from ETH Zurich, Prof. Tom Chou from UCLA and Dr. B\"ottcher from University of Frankfurt for their support and helpful discussions.
 \section{Methods}
This section introduces some fundamental concepts and mathematical definitions used throughout the project, ranging from \CDR sequence generation to fundamentals on sets, statistics and network theory to details about our network construction framework, in particular with regard to the GPU implementation.
\subsection{Generation of \CDR sequences}
During V(D)J recombination, the probability that a certain TCR  sequence is formed by a recombination event $E$ can be approximated by the equation~\cite{sethna2020population,murugan2012statistical}
\begin{align}
\begin{split}
P_{\rm recomb}(E) & = P_V(V)P_{DJ}(D, J)\\
&\cdot P_{{\rm del}V}(d_V|V)P_{{\rm del}J}(d_J|J)P_{{\rm del}D}(d_D, d_D'|D)\\
&\cdot P_{{\rm ins}VD}(l_{VD})\prod\limits_{i=1}^{{\rm ins}VD}p_{VD}^{(2)}(x_i|x_{i-1})P({\rm ins}DJ)\\
&\cdot \prod\limits_{i=1}^{{\rm ins}DJ}p_{DJ}^{(2)}(y_i|y_{i-1})
\end{split}
\label{eq:Pscenario}
\end{align}
Here, $V$, $D$, and $J$ are gene sequences from the variable, diversity, and joining regions, respectively. $P_V(V)$ is the probability, that a certain sequence $V$ is chosen from the variable region. $P_{DJ}(D, J)$ is the joint probability, that the genes $D$ and $J$ are chosen from the diversity and joining region. $P_{{\rm del}V}\left(d_V| V\right)$ is the conditional probability that a number $d_V$ of bases is deleted at the end of the $V$ sequence. The probabilities $P_{{\rm del}D}\left( d_D, d_D'|D \right)$ and $P_{{\rm del}J}(d_J| J)$ are defined analogously, where $d_J$ is the number of bases that is deleted at the end of the $J$ section. Since the $D$ region is located at the center of a TCR sequence, $d_D$ and $d_D'$ correspond to the number of deleted bases at the two ends of the sequence. Moreover, $P_{{\rm ins}VD}(l_{VD})$ and $P_{{\rm ins}DJ}\left(l_{DJ}\right)$ are the probabilities that a certain sequence of length $l_{VD}$ or $l_{DJ}$ is inserted at the $V$-$D$ junction or $D$-$J$ junction, respectively.
The conditional probability $p_{VD}^{(2)}(x_i | x_{i-1})$ is the probability that the base at position $i$ has the value $x_i$ given the value $x_{i-1}$ of the base at position $i-1$. This accounts for a nearest-neighbor dependency of the base selection. Analogously, $p_{DJ}^{(2)}(y_{i}|y_{i-1})$ is the probability for the base at position $i$ to have the value $y_i$, while the base at position $i-1$ has the value $y_{i-1}$. The total number of parameters is 2,856~\cite{murugan2012statistical}, and most of them are associated with the deletion probabilities.
\subsection{IGoR, OLGA, SONIA}
\label{sonia}
During T-cell formation, two important processes have to be disentangled. One is the generation of a distribution of T-cell clones due to V(D)J recombination in the thymus. The second process is the selection of clones in the thymus based on their
reactive behavior. Both processes can be described with mathematical models, embedded in software packages. The first process is captured by the IGoR package~\cite{marcou2018high}, which learns the model based on empirical data by adjusting the parameters based on the comparison of generated sequences and experimentally obtained reference sequences. The OLGA algorithm~\cite{10.1093/bioinformatics/btz035} increases the efficiency and thus performance of this task. SONIA is a software tool that additionally incorporates the second process, that is, thymus selection~\cite{sethna2020population}.
 \subsection{Histogram and relative frequency}
 \label{Histogram and relative frequency}
 May $\boldsymbol{x} = \left(x_1, x_2, ...,  x_n\right)$ be a vector of $n$ elements, each chosen from a set $\boldsymbol{\tilde{x}} = \{\tilde{x}_i\ |\ 0<i\leq \tilde{n}\}$ of $\tilde{n}$ different elements  $\tilde{x}_i$.
 Then the  {\itshape histogram function} $m(\tilde{x}\in \boldsymbol{\tilde{x}})$ (or short  {\itshape histogram})  can be defined as
 \begin{equation}
     m(\Tilde{x}) = |\{i\in \mathbb{N}_{|\boldsymbol{x}|}:x_i = \Tilde{x}\}|,
 \end{equation}
 where $|\cdot|$ denotes the {\itshape cardinality} of a set.
 We define the  {\itshape relative frequency $f$} of an element $\Tilde{x}\in \boldsymbol{\tilde{x}}$ as
 \begin{equation}
     f(\tilde{x}) = \frac{ |\{i\in \mathbb{N}_{|x|}:x_i = \Tilde{x}\}|}{|\boldsymbol{x}|} = \frac{m(\tilde{x})}{|\boldsymbol{x}|}.
 \end{equation}
 This normalization assures $\sum\limits_{\tilde{x} \in \boldsymbol{\tilde{x}}} f(\tilde{x}) = 1$. We also refer to the distribution of relative frequencies as the {\it normalized histogram}.
 \subsection{Cartesian product}
 \label{Cartesian product}
 The {\itshape Cartesian product} $C$ between two sets $A$ and $B$ is defined as
 \begin{equation}
     C = A\times B = \{(a, b)|\ a\in A\ \textnormal{and}\ b\in B\},
 \end{equation}
 where $C$ is a new set, comprising tuples $(a, b)$ as elements. For instance,
 the Cartesian product of the two sets $A = \{1, 2, 3\}$ and $B = \{\alpha, \beta, \gamma \}$ is
 \begin{equation}
 \begin{split}
     C = \{(1, \alpha), (1, \beta),  (1&, \gamma), (2, \alpha), (2, \beta),  (2, \gamma), \\
     &(3, \alpha), (3, \beta), (3, \gamma)\}.
  \end{split}
 \end{equation}

 The n-ary  {\itshape Cartesian power} of a set $X$ is defined as
 \begin{equation}
 \begin{split}
     X^{n} &= \underbrace{X\times X\times ... \times X}_{n} \\
     &= \{(x_1, x_2, ...,  x_n)\ |\  x_i\in X\ \forall i: 0<i\leq n\}.
 \end{split}
 \end{equation}

 \subsection{Combination}
 \label{Combination}
 May $\boldsymbol{x}$ denote a set containing $N_{max}$ elements.
 Moreover, $\boldsymbol{y}$ shall be a set of $N$ elements, each chosen from $\boldsymbol{x}$,  where the same element can be chosen multiple times, thus $N$ can be larger than $N_{max}$. Then the function $C_{\boldsymbol{p}}(N, N_{max}, \boldsymbol{x})$ shall be defined as the mapping from $\boldsymbol{x}$ to $\boldsymbol{y}$, where
 $\boldsymbol{p} = (p_1, p_2, ..., p_{N_{max}})$  denotes the vector of probabilities $p_i$ to choose an element $x_i\in \boldsymbol{x}$. In the next section we will write  $C_{p}(N, N_{max})$ or $C_{p}$ for brevity and make the input set $\boldsymbol{x}$ implicit since we are more interested in the mapping function itself than in the concrete elements.

 \subsection{Combination of amino acids}
 \label{Combination of amino acids}
 We represent each amino acid
 with its corresponding 1-letter symbol. Then the 20 different amino acids (a.a.) form the set
\begin{align*}
    \Sigma = \ &\{A, R, N, D, C, Q, E, G, H, I, L, \\
    &K, M, F, P, S, T, W, Y, V\}.
\end{align*}
Thus, a sequence $\boldsymbol{\sigma}$ containing a number $l$ of a.a. is an element of the set  $\Sigma^{l}$ of all possible a.a. sequences with length $l$. The total set of all possible a.a. sequences with arbitrary length then can be defined as $\Sigma_{total} = \bigcup_{l=1}^{\infty}\Sigma^{l}$. In reality, not all of these a.a. sequences are possible during V(D)J recombination. Instead,  the choice reduces to a subset $\Sigma_{N_{max}} \subset \Sigma_{total}$ of $N_{max}$ possible a.a. sequences. Then an array of $N$ a.a. sequences, generated during V(D)J recombination, is an element of the set $\Sigma_{N_{max}, N} := \Sigma_{N_{max}}^{N}$ (the N-ary Cartesian power of $\Sigma_{N_{max}}$ - see section \ref{Cartesian product}).

The combination $C_{\boldsymbol{p}}$ uses $N \in \mathbb{N}$ and $N_{max} \in \mathbb{N}$ as an input in order to generate an array of $N$ a.a. sequences $\{\boldsymbol{\sigma}\}_N:=\{\boldsymbol{\sigma}_i |\ 0<i\leq N\} \in \Sigma_{N_{max}}^{N}$. We view the vector $\boldsymbol{p}$ of weights of all possible a.a. sequences as a parameter instead of an input. Thus, the combination can be defined as

\begin{align*}
C_{\boldsymbol p}: & \ \mathbb{N} \times \mathbb{N} & \rightarrow & \  \bigcup_{0<N<\infty}\ \bigcup_{0<N_{max}<\infty}\Sigma_{N_{max}}^{N} \\
   \ & \left(N, N_{max}\right) & \mapsto & \ \{\boldsymbol{\sigma}\}_N .
\end{align*}

\label{dynamic}

\subsection{Overview of network theory concepts}
\label{complex_networks}

A graph (\ie, network) $G(V, E)$ is an ordered pair $(V,E)$,
where $V$ is the set of nodes and $E\subseteq V\times V$ is the
set of edges that connect pairs of
nodes~\cite{newman2018networks}. In this work, we consider
undirected and unweighted networks. The connectivity properties
of such a network are described by the symmetric {\it adjacency matrix} $A$ with binary entries $0$ and $1$. If an edge connects two nodes $u,v\in V$, the corresponding adjacency matrix entry $A_{uv}=A_{vu}$ is 1 and 0 otherwise. For each node $u \in V$, the {\it degree} is given by the number of edges attached to that node, that is, $d(u)=\sum_{v\in V} A_{u
v}=\sum_{v\in V} A_{v u}$. The {\it degree distribution} $p_d$ is the relative frequency of nodes with degree $d$ and defines the mean degree according to $\langle d \rangle=\sum_{d} d p_d$. The {\it clustersize distribution} is defined as the relative frequency of clusters with size $CC$.
The {\it shortest path length} between a source node $s$ and a target node $t$ is the minimum number of edges that connect these nodes. If the shortest path length is finite for every node pair $(s,t)\in V\times V$, the network is {\it connected}. A network is {\it fully connected} (or complete) if every node is connected to every other node. A {\it subgraph} of $G(V,E)$ is a network $G_{sub}(V',E')$ over a subset $V'\subseteq V$ of nodes and the subset $E'\subseteq E$ of all edges that connect nodes in $V'$. The {\it largest connected component} (LCC) is the largest connected subgraph.
The {\it relative size} of a subgraph $G_{sub}$ of $G$
is defined as the quotient $\frac{|G_{sub}|}{|G|}$, where $|G_{sub}|$ and $|G|$ denote the number
of nodes in the subgraph and graph respectively (also called the {\it size} of a graph).

The {\it degree assortativity} $r$ is defined as~\cite{newman2002assortative}
\begin{equation}
r = \frac{\sum_{kl}kl\left(e_{kl} - q_k q_l\right)}{\sigma_q^2},
\end{equation}
where $q_k=(k+1)p_{k+1}/\langle k \rangle$ is the distribution of {\it excess degree}. The {\it average neighbor degree} of a network is
\begin{equation}
\langle d_{nn} \rangle = \sum\limits_{d'}d'P\left(d'| d\right),
\end{equation}
where $d$ is the degree of a node, $d'$ is the degree of a neighboring node, and $P(d'|d)$ is the conditional degree distribution.
The {\it betweenness centrality} of a graph is defined as
\begin{equation}
g(v) = \sum\limits_{s\neq v\neq t} \frac{\sigma_{st}(v)}{\sigma_{st}},
\end{equation}
where $\sigma_{st}$ denotes the number of shortest paths from a node $s$ to a node $t$ and $\sigma_{st}(v)$ is the number of shortest paths, which cross the node $v$.

\subsection{Percolation transition}
\label{transition}
An interesting phenomenon in many networks, e.g. Erd{\H o}s--R{\'e}nyi graphs, is the percolation phase transition.
We distinguish {\itshape bond percolation} and {\itshape site percolation}. While the former process describes the change of network properties during removal of edges, the latter describes changes under node removal.
Thus, we randomly choose and remove a fraction $p$ of all edges or nodes respectively, while observing how a network property, e.g. the LCC, changes for higher $p$. In the following we focus on bond percolation.
Erd{\H o}s--R{\'e}nyi graphs, as many other graph types, exhibit a transition at a $p_c < 1$ where the relative size of the LCC vanishes~\cite{tishby2018revealing}. This value lies at $p_{c,ER} = \frac{1}{\langle k \rangle}$, where $\langle k\rangle$ refers to the average node degree in the Erd{\H o}s--R{\'e}nyi graph. In contrast, Barabási--Albert graphs do not exhibit this transition for a  $p_c$ below one~\cite{PhysRevLett.85.4626}. Instead, the LCC continuously shrinks and only vanishes when the number of edges goes to 0. Thus, $p_{c,BA} = 1$.

\subsection{Normalization}
\label{normalization}
Throughout the analysis we employed different normalization procedures, which are defined below.
To normalize the sharing number distribution, shown in the panels (a) and (b) of Fig. \ref{fig:fig_2}, the frequency $f(k)$ of \CDR sequences with a certain sharing number $k$  will
be divided by the total number $N$ of sequences:
\begin{equation}
    f_{norm} = \frac{f(k)}{N}
\end{equation}
For the normalization procedure of the betweenness distribution and the degree distribution, we refer to the histogram normalization, described in section \ref{Histogram and relative frequency}.
\subsection{Computing LCC}
\label{computing_LCC}
To normalize the LCC in the robustness analysis in panel (g) of Fig. \ref{fig:fig_2}, we first divide the number of nodes in the LCC (denoted by LCC$(p)$) by the number of
 nodes in the network after node removal, given by $(1-p)N$, i.e.,
\begin{equation}
    \textnormal{LCC}_{norm,1}(p) = \frac{\textnormal{LCC}(p)}{(1-p)N}
\end{equation}
In a second step we divide by the relative size of the LCC at $p=0$, i.e., when
no nodes were removed:
\begin{equation}
\textnormal{LCC}_{norm, 2}(p) = \frac{\textnormal{LCC}_{norm,1}(p)}{\textnormal{LCC}_{norm,1}(0)}.
\end{equation}
In panel (i) of Fig. \ref{fig:fig_2} we don't remove any nodes, only perform the first normalization and therefore plot $\textnormal{LCC}_{norm,1}(p=0) = \frac{\textnormal{LCC}(0)}{N}$.

\subsection{Averaging}
\label{averaging}
In order to reduce fluctuations, we can  average over multiple samples to calculate a distribution $y = f(x)$.
May each sample, referred to via index $i\in \{1, 2, ..., M\}$, where $M$ denotes the number of samples, have its unique distribution $f_i(x)$. Then the average over all distributions is calculated as
\begin{equation}
    \bar{f}(x) = \frac{1}{M}\sum\limits_{i=1}^{M}f_i(x).
\end{equation}

\subsection{Levenshtein distance}
\label{levenshtein_distance}

The {\itshape Levenshtein distance} can be used to compare two strings $a$ and $b$, each having arbitrary length, and is defined as~\cite{levenshtein1966binary}

\begin{align}
\resizebox{.9\hsize}{!}{
$\textnormal{lev}_{a, b}\left( i, j\right) = \left \{
\begin{array}{cc}
\max \left( i-1, j-1 \right)  &\textnormal{if $i=1$}\\ &\textnormal{or $j=1$}\\
\min \left\{
\begin{array}{c}
\textnormal{lev}_{a,b} \left( i-1, j \right) + 1\\
\textnormal{lev}_{a,b} \left(i, j-1\right) + 1\\
\textnormal{lev}_{a,b} \left(i-1, j-1\right) + 1_{\left(a_i \neq b_j\right)}
\end{array}
\right.
& \textnormal{else}
 \end{array}
\right . $
}
\label{eq:levenshtein}
\end{align}

with

\begin{equation}
1_{\left(a_i \neq b_j\right)} = \left\{
\begin{array}{cc}
0&\textnormal{if $a_i = b_j$}\\
1&\textnormal{ otherwise.}
\end{array}
\right.
\end{equation}

$a_i$ denotes the $i$th entry of the string $a$. Likewise, $b_j$  corresponds to the $j$th entry of $b$.
If $a$ and $b$ are equal in length, the Levenshtein distance is also called {\itshape Hamming distance} and corresponds to the number of elements in which both strings differ. E.g. if $a=\text{ACCA}$ and $b=\text{ACGG}$, both strings only differ at positions three and four. Hence, their Hamming distance corresponds to two.
Since \eqref{eq:levenshtein} is recursive, we performed the computation of the Levenshtein distance iteratively, based
on a straight-forward algorithm.
\subsection{TCR network}
\begin{figure*}
\includegraphics[scale=0.9]{example_graph.pdf}
\caption{
In panels (a) to (f)
we draw the largest connected component (LCC) of TCR networks with different sizes $N$ but constant $l_{max} = 6$. In each graph, the color indicates the degree, where the node with maximum degree occurs in yellow and the nodes with lowest degree  in dark blue.
}
\label{fig:graph}
\end{figure*}
We define a TCR network as a graph $G$ where each \CDR sequence $\boldsymbol{\sigma}$ is represented by a node, also called TCR node (see graphs of different sizes in Fig. \ref{fig:graph}). If and only if the Levenshtein distance $l$ between two \CDR sequences does not exceed a threshold value $l_{max}$ (thus $l \leq l_{max}$), the two corresponding nodes are connected via an edge.

 \subsection{Adjacency matrix}
 Let's consider a graph $G$ with $n$ nodes, connected via a set $E$ of edges. Each edge $e\in E$ can be represented via a tuple of two numbers $(n_1, n_2) \in \mathbb{N}_n\times \mathbb{N}_n$, where $n_1$ and $n_2$ denote the indices of the two connected nodes. Since the edges have no direction ($G$ is an {\itshape  undirected} graph), switching the two entries of a tuple yields the same edge, i.e. $(i, j) \cong (j, i)\ \forall {i, j \in \mathbb{N}_n}$.
 The edges can be represented by a symmetric matrix, referred to as the {\itshape  adjacency matrix}, where

 \begin{align*}
     (i, j) \mapsto \begin{cases} 0, \ &\textnormal{if} \ (i, j)\notin E \ \textnormal{and} \ (j, i)\notin E \\
     1, \ &\textnormal{if} \ (i, j)\in E \ \textnormal{or} \ (j, i)\in E.\end{cases}
\end{align*}

 For example, let us consider a  graph $G$ with five nodes and the edges
 \begin{equation*}
     E =  \{(1, 2), (1, 3), (2, 5), (4, 5)\}.
 \end{equation*}
 Then, the adjacency matrix becomes
 \begin{equation}
     \begin{pmatrix}
     0&1&1&0&0\\
     1&0&0&0&1\\
     1&0&0&0&0\\
     0&0&0&0&1\\
     0&1&0&1&0\\
     \end{pmatrix}_{.}
 \end{equation}
 For larger graphs, this matrix increases in scale and the number of matrix entries (and thus the computational effort to construct the matrix) increases quadratically with the $n$.

\subsection{GPU-based network generation}
\label{sec:GPU_algorithm}
\begin{figure*}
    \centering
    \includegraphics[scale=1.0]{ims/performance.pdf}
    \caption{{\bf Performance of TCR-NET.} (a) The runtime to calculate Levenshtein distances between all possible pairs of $N$ \CDR sequences using the CPU-based algorithm {\sc imNet}~\cite{miho2019large} and our GPU-based implementation TCR-NET with $l_{\rm max} = 1$. (b) The computational speedups associated with using TCR-NET for Levenshtein distance calculations. All calculations were performed on one core of an Intel i7-7700 CPU and an Nvidia GTX 1660 GPU. The GPU speedup indicates a lower bound of the number of CPU cores necessary to outperform GPU calculations for the same number of sequences.}
    \label{fig:performance}
\end{figure*}
The advantage of GPUs is their ability to parallelize processes, what
can be utilized to calculate the Levenshtein distance between sequences and with that the entries of the adjacency matrix, more efficiently, as done previously~\cite{balhaf2016using}.
We used a Nvidia GTX 1660, with  $t=1408$  kernels and 6 GB GDDR6.  For large networks its memory does not suffice to store the entire adjacency matrix.
In order to overcome this restriction, we employ the strategy of {\itshape tiling},
 where we divide the total adjacency matrix into blocks,
which are calculated consecutively. To further reduce the memory consumption, each block is saved as a sparse matrix.
The GPU speedup is compromised by the
lower operating frequency of GPU kernels compared to CPU kernels and the
transfer time between the GPU and the CPU RAM.
 All this has to be accounted for while implementing an efficient GPU
 algorithm.

 Knowing that most \CDR sequences are only 8 to 11 a.a. long
 and the relative frequency of sequences declines sharply
 with higher length, we reserved memory space for only 30 a.a. per sequence on the GPU
  in order to limit the computational workload.
  Moreover, there are only 20 different common a.a. Each a.a. is saved as a number. When they are converted, these numbers are not saved as {\sc int} (32 bit) but instead as {\sc byte} (8 bit),
  which occupies only a quarter of the memory on the GPU.
  We use numba to perform distance calculations on the GPU and  determine the sparse adjacency matrix.
  All these processes have
  the goal of improving the computing efficiency.
In order to further reduce the computation time the iterative calculation of the levenshtein distance between two a.a. sequences, will abort, as soon as the smallest intermediate result exceeds the maximum allowed levenshtein distance $l_{max}$. The algorithm to calculate the levenshtein distance is shown as pseudocode in Listing \ref{listing1}.

 In Fig. \ref{fig:performance} we compare the
 performance of {\sc TCR-NET} with the previously published high-performance framework {\sc  imnet}~\cite{miho2019large}. TCR-NET beats {\sc imnet} on a normal computer with a single core of an Intel i7-7700 CPU by approximately one to two orders of magnitude, with a higher speedup for large networks. Of course {\sc imnet}'s performance has an advantage on large CPU clusters, where the high number of kernels  enables an efficient use of the pyspark parallelization. Moreover, TCR-NET only
 calculates efficiently if the number of edges constitutes merely a small fraction of all possible edges, which is usually true for a TCR network, if $l_{max}$ is low.

 We employed {\sc Python} as the main programming language and used  numba and pycuda for GPU optimization. For complex network computation, the package {\sc networkx} was  employed in order to supplement the packages {\sc numpy} for mathematical computation and {\sc matplotlib} for
plotting.
The {\sc imnet} package employs igraph and R
\cite{rcore2016r}
\cite{csardi2006igraph}.
With the {\sc SONIA} package
 a.a. strains could be generated directly from the commmand line~\cite{sethna2020population}.

The full source code and additional information are publicly available on GitHub~\cite{GitHub}.
 \lstset{language=Python}

 \begin{lstlisting}[caption = {Pseudocode for the calculation of the Levenshtein distance. The source code and additional information is publicly available on GitHub~\cite{GitHub}.}, basicstyle=\small, frame = single]
t = len(string_1)
s = len(string_2)

dnew = zeros(s);
dprev = zeros(s);

# info: calculating
#     the Levenshtein distance
for (i = 1; i < len(t); i++):
    dprev[0] = i-1;
    dnew[0] = i;

    for (j = 1; j < len(s); j++):
        if(string1[i-1] == string2[j-1]):
            substitutionCost = 0
        else:
            substitutionCost = 1

        dnew[j] =
            min(dnew[i-1])+1, dprev[i]+1,
                dprev[i-1] +
                substitutionCost)

        if (dnew[j] < min_val):
            min_val = dnew[j]

    for (j = 1; j < len(s); j++):
        dprev[j] = dnew[j];

    if (min_val > max_ld):
        break;
# info: checking, if the Levenshtein
#     distance is in the desired
#     range; c_value represents the
#     entry in the adjacency
#     matrix
if dnew[len(s)-1] < max_ld:
    c_value = 1
else:
    c_value = 0

 \end{lstlisting}
\label{listing1}
\bibliographystyle{apsrev4-1}


\section*{Data availability statement}
The datasets generated during and/or analysed during the current study are available in the Open Science Framework (OSF) repository,

\verb+https://osf.io/jnx8k/?view_only=16426f9a35404264badaaa93162060a7+

\renewcommand{\thefigure}{S\arabic{figure}}
\setcounter{figure}{0}
\renewcommand{\theequation}{S\arabic{equation}}
\setcounter{equation}{0}

\section*{Supplementary Information}
In this Supplementary Information, we provide more details on the behavior of the Hermitian part of the Hamiltonian of our system. We demonstrate how the Hamiltonian parameters are extracted from the experiment, including the non-Hermiticity (linear dichroism). Finally, we demonstrate that the winding of the real effective field constrains the possibility of a non-Hermitian topological transition with annihilation of the exceptional points: The annihilation is only possible if the total winding of the real effective field is zero.

\subsection{Topology of the Hermitian Hamiltonian: Theory}

In this section, we illustrate the topological transitions concerning the Hermitian part of the two-band Hamiltonian (1) of the main text describing the intersection of cavity modes with different numbers. This Hamiltonian is shown below for convenience, together with the definitions:
\begin{equation}
H_{\bf{k}}^{real} = \left( {\begin{array}{*{20}{c}}
{\frac{{E_H^{N+2} + E_V^{N}}}{2} + \frac{{{\hbar ^2}k_x^2}}{{2{m_x}}} + \frac{{{\hbar ^2}k_y^2}}{{2{m_y}}}}&{\Delta  - \beta '{k^2} - \beta {{({k_x} - i{k_y})}^2}}\\
{\Delta  - \beta '{k^2} - \beta {{({k_x} + i{k_y})}^2}}&{\frac{{E_H^{N+2} + E_V^{N}}}{2} + \frac{{{\hbar ^2}k_x^2}}{{2{m_x}}} + \frac{{{\hbar ^2}k_y^2}}{{2{m_y}}}}
\end{array}} \right),
\end{equation}
where $E_H^{N+2}$ and $m_H\sim N+2$ are the energy and mass, of the $N+2$th H-polarised mode and $E_V^{N}$ and $m_V\sim N$, are the energy and mass of the V-polarised mode number $N$. $k_x$, $k_y$ are the 2D wave vector components.
The spin-independent masses $m_x$ and $m_y$ are determined by the birefringence ($n_e$, $n_o$) and the angle of the optical axis $\theta$, according to the following expressions \cite{rechcinska2019engineering}:
\begin{eqnarray}
    \frac{1}{m_x}&=&\frac{1}{2m'_0}\left(\frac{n_o^2\sin^2\theta+n_e^2\cos^2\theta+n_e^2}{n_e^2 n_o^2}\right),\\
    \frac{1}{m_y}&=&\frac{1}{2m'_0}\left(\frac{n_e^2+n_o^2}{n_e^2 n_o^2}\right),
\end{eqnarray}
where $m'_0$ is the mass at the central frequency (between the two modes in question).
$\beta$ is the magnitude of the TE-TM spin orbit coupling. $\beta'=\frac{\hbar^2(m_V-m_H)}{2m_Hm_V}$ and $\Delta=(E_H^{N+2}-E_V^{N})/2$.
This Hermitian Hamiltonian can be written as a linear combination of identity and Pauli matrices which defines a real effective magnetic field $\Omega_r$ acting on the polarization pseudo-spin. The two non-zero components of the field are $\Omega_r^x=\Delta-\beta'k^2-\beta(k_x^2-k_y^2)$ and  $\Omega_r^y=-2 \beta k_x k_y$.

The pseudospin of a given eigenstate is either aligned or anti-aligned with the effective field (so far as the Hamiltonian remains Hermitian). Experimentally, this pseudospin corresponds to the Stokes vector of light, calculated using the following expressions:
\begin{eqnarray}
    S^1&=&\frac{I_V-I_H}{I_H+I_V},\\
    S^2&=&\frac{I_D-I_A}{I_D+I_A},\nonumber\\
    S^3&=&\frac{I_R-I_L}{I_R+I_L},\nonumber
\end{eqnarray}
where the $I$s are the intensities in each of the polarizations, which can be measured experimentally or calculated theoretically as the absolute value squared of the eigenstate wavefunction in the corresponding basis. Each pseudospin component is an average value of the corresponding spin operator represented by a Pauli matrix.

The case $\beta'=0$ has already been considered in \cite{Tercas2014,gianfrate2020measurement}. The two parabolas are split at $k=0$, with different possible signs of this splitting (band detuning).

\begin{figure}
    \centering
    \includegraphics[width=0.95\linewidth]{figS1.pdf}
    \caption{\textbf{Pseudospin textures for the Hermitian Hamiltonian}: \textbf{a}  trivial case $\beta<\beta'$, \textbf{b,c}  topologically non-trivial case $\beta>\beta'$ with positive and negative detunings.  The red marks shows the Dirac points of winding number +1 and the blue marks the Dirac points of winding number -1}
    \label{texture}
\end{figure}

For positive detuning, the TE-TM splitting compensates this splitting along $k_x$ leading to a linear crossing.
On the contrary, the TE-TM splitting is changing sign along $k_y$ and adds to the $k=0$ splitting. As a result, along $k_y$ the two parabola get away the one with respect to the other and do not cross. The resulting spectrum shows two tilted Dirac cones at $\pm k_{0x}$. The winding number of the effective field around these two points is the same (+1), as shown in Fig.~{\ref{texture}b} showing the distribution of the pseudospin (and effective field) texture in the reciprocal space. This winding is inherited from the winding number 2 of the TE and TM modes. The Berry topological charge is $+1/2$. This topological phase remains while $\beta'$ increases up to $\beta'=\beta$ where two new crossing points between the parabola appear at $k_{0y}$ which sets at infinity when the transition occurs. These crossing points correspond to tilted Dirac cones. The winding number of the in-plane pseudo-spin around these two points is (-1) and the Berry topological charge (-1/2), as shown on Figure~\ref{texture}-b. The sum of the topological charges of the 4 Dirac points is zero and the bands are overall topologically trivial which allows the non-Hermitian topological transition we explore in the main text.

For negative detuning (Fig.~\ref{texture}c), the winding of the two Dirac points does not change, but their positions in reciprocal space are rotated by $\pi/2$ and the particular field distribution changes from monopolar to Rashba.

\subsection{Topology of the Hermitian Hamiltonian: Experiment}

\begin{figure}[tbp]
    \centering
    \includegraphics[width=0.8\linewidth]{figS2.pdf}
    \caption{\textbf{4 diabolical points case}. Experimental dispersion relation along \textbf{a} $k_x$ and \textbf{b} $k_y$ directions observed at 1.39\,V applied to the microcavity, corresponding to the case with total winding 0. Energy of the bands visible in horizontal (vertical) polarization is marked by green (purple) dashed line. \textbf{c} Experimental and \textbf{d} simulated by Berreman method pseudospin texture in $S_1$--$S_2$ plane of lower energy band.}
    \label{exp_N2N}
\end{figure}

\begin{figure}[tbp]
    \centering
    \includegraphics[width=0.8\linewidth]{figS3.pdf}
    \caption{\textbf{2 diabolical points case}. Experimental dispersion relation along \textbf{a} $k_x$ and \textbf{b} $k_y$ directions, observed at 11\,V applied to the microcavity, corresponding to the case with winding 2. Energy of the bands visible in horizontal (vertical) polarization is marked by green (purple) dashed line. \textbf{c} Experimental and \textbf{d} simulated by Berreman method pseudospin texture in $S_1$--$S_2$ plane of lower energy band.}
    \label{exp_NN}
\end{figure}

Cavities filled with birefringent nematic liquid crystal allows for realization of both topologically-distinct cases mentioned above: with non-trivial total winding 2 and  trivial total winding 0 (discussed in the main text), depending on voltage applied to the cell. For large enough energy detuning between the two modes of interest, the relative non-Hermiticity is weak enough so that the measured pseudo-spin is simply aligned with the real effective magnetic field (Hermitian limit). In such a case, the winding of the real field in reciprocal space can be determined experimentally through polarization resolved tomography.

Experiments described in this section were preformed on a sample analogous to the one described in the main text, but with wider LC layer of around 3.9\,$\mu$m.

Fig.\,\ref{exp_N2N} presents realization of the case with winding 0, discussed in detail in the main text. Fig.\,\ref{exp_N2N}a,b presents experimental dispersion relation along $k_x$ and $k_y$ directions observed with 1.39\,V applied to ITO electrodes. At this voltage, the numbers of the modes differ by 2, and both modes cross each other along both wave vector directions.  Experimentally determined pseudospin texture in $S_1$--$S_2$ plane of the lower energy band is plotted in Fig.\,\ref{exp_N2N}c. It  matches the Berreman matrix model presented in Fig.\,\ref{exp_N2N}d, showing 0 total winding.

With amplitude of external voltage of 11\,V modes with the same mode numbers are close to degenerate, as shown on experimental dispersion relations in Fig.\,\ref{exp_NN}. In this case modes crosses each other only along $k_y$ direction. Resulting  pseudospin texture in $S_1$--$S_2$ plane shows two 2D monopoles with total winding number of 2 in both experiment Fig.\,\ref{exp_NN}c and Berreman model simulations Fig.\,\ref{exp_NN}d.

The Dirac points appear along the Y axis because of the negative value of the detuning in the experiment, but their winding is that of the TE-TM field, as expected.

\subsection{Hamiltonian parameters}

\begin{figure}[tbp]
\includegraphics[width=0.95\linewidth]{figS4.pdf}
\caption{\textbf{Fit of the dispersions.} Fit of the polarization-resolved dispersions in the two perpendicular directions. The false color map shows the reflectance as a function of wave vector and energy. Fitting allows to obtain the spin-independent masses $m_x$ and $m_y$, as well as the polarization splittings $\beta$ and $\beta'$. \label{figDisp1}}
\end{figure}
The parameters of the Hermitian Hamiltonian (1) are extracted from the fit of the dispersion in the two perpendicular directions. The fits are shown in Fig.~\ref{figDisp1}.
In a given direction, the average of the two parabolicity coefficients (inverse masses) of the two polarizations, provides the value of $1/m_x$ and $1/m_y$. Then, the difference between the two parabolicites in the two polarizations is $\beta+\beta'$ for $Y$ direction and $\beta-\beta'$ for $X$ direction. Fitting both allows to extract $\beta$ and $\beta'$ separately.

\subsection{Non-Hermitian parameter}

In Fig.~\ref{DGamma}, we show an example of a spectrum of reflectance in two polarizations, clearly demonstrating different linewidths. The experimental points for each of the two polarizations are fitted with the Voigt function, as explained in Methods. We stress that while the non-Hermiticity of the Hamiltonian is a single constant, the linewidth at each particular wave vector is determined by the imaginary part of the corresponding eigenstate. Only at the positions of the Dirac points of the Hermitian Hamiltonian does the linewidth of the eigenstates correspond directly to the non-Hermiticity of the Hamiltonian $\delta \Gamma$. It is also the maximal possible imaginary part of the energy that could be observed for any wave vector. This is what we use to extract the non-Hermitian parameter $\delta\Gamma$.

\begin{figure}
    \centering
    \includegraphics[width=0.8\linewidth]{figS5.pdf}
    \caption{\textbf{Polarization-resolved reflectance spectrum.} The reflectance in polarizations H and V, measured at the position of the diabolical point and plotted as a function of energy. The difference of linewidths is approximately $13\%$.}
    \label{DGamma}
\end{figure}

\subsection{Hermitian winding and annihilation of exceptional points}

We consider a two-band system described by a Hamiltonian depending on a two-dimensional wave vector. The Hermitian part of the Hamiltonian is characterized by a presence of diabolical (Dirac) points. The two-band nature allows writing the Hamiltonian as a superposition of Pauli matrices. For many typical situations, such Hamiltonians can be written using only two Pauli matrices that we can choose to be $\sigma_x$ and $\sigma_y$ (massless Dirac Hamiltonian, Rashba and Dresselhaus spin-orbit couplings, TE-TM spin-orbit coupling for photons), whereas the third Pauli matrix $\sigma_z$, responsible for a symmetry breaking, opens the gap at the diabolical points (the mass term in the Dirac Hamiltonian or the Zeeman splitting for electrons and photonic modes). In this case, it is possible to characterize each of the diabolical points and the whole reciprocal space by winding numbers, which, in turn, determine the topology (the Chern number) of the bands once the gap is opened.

We then consider the evolution of the Hamiltonian with an addition of non-Hermiticity, described by a single constant parameter (independent of the wave vector). Without loss of generality, we consider the non-Hermitian contribution to be described by the Pauli matrix $\sigma_x$. In presence of the non-Hermiticity, each of the diabolical points gives rise to two exceptional points.

Thanks to the use of the Pauli matrices, it is possible to map the system's Hamiltonian to a spin-$1/2$ in a magnetic field, using the so-called pseudospin formalism, where the terms of the Hamiltonian are considered as an effective field. In the presence of a non-Hermiticity, this effective field has both real and imaginary components, $\bm{\Omega}=\bm{\Omega}'+i\bm{\Omega}''$. An exceptional point requires the real and imaginary parts of the effective field to be equal to each other in magnitude and perpendicular in direction: $\Omega'=\Omega''$ and $\bm{\Omega}'\perp \bm{\Omega}''$.

\begin{figure}[tbp]
\includegraphics[width=1.05\linewidth]{figS6.pdf}
\caption{\textbf{Hermitian winding and annihilation of exceptional points.} (a) Non-zero total winding of the real field, the exceptional points do not annihilate. (b) Zero total winding of the real field, annihilation of exceptional points. The trajectories of the exceptional points are shown with red and blue arrows (corresponding to the winding of the exceptional points themselves). Crosses mark the annihilation points. \label{Sfig1}}
\end{figure}

Since in our case the direction of the imaginary field is fixed $\bm{\Omega}''=\Omega''\bm{e}_x$, and the real field is in the XY plane, the exceptional points can only be found where the x-component of the real field is zero $\bm{\Omega'}=\Omega'\bm{e}_y\Leftrightarrow\Omega'_x=0$. This condition constrains their evolution in the reciprocal space when the imaginary field is increased $\Omega''\to\infty$ (or the real field is decreased $\Omega'\to 0$). The location of the ensemble of the points $\Omega'_x=0$ is determined by the winding of the real field at the diabolical points. With the increase of $\Omega''$, the exceptional points are moving towards infinity. \textbf{If no points, where the condition for the real part $\Omega'_x=0$ is satisfied, can be found at infinity, then the exceptional points must annihilate for some finite critical value $\Omega''_c$.} Two cases are possible:
\begin{enumerate}
    \item The overall winding of the real field is a non-zero integer. In this case, the orientation angle $\theta_r$ of the real field at infinity changes at least once between $\theta_r=0$ and $\theta_r=2\pi$, necessarily going through $\pi/2$ and $3\pi/2$, where $\Omega'_x=0$. Therefore, the exceptional points can be found for any value of $\Omega''$ up to infinity, and thus they cannot annihilate (or if they annihilate accidentally, they must then reappear when $\Omega''$ is increased).
    \item The overall winding of the real field is zero. In this case, the orientation angle of the real field $\theta_r$ at infinity is constrained to an interval which may not include $\pi/2$ or $3\pi/2$. In this case, the exceptional points cannot be present at infinity, and therefore they must annihilate for a finite value of $\Omega''_c$.
\end{enumerate}

The two situations are illustrated in Fig.~\ref{Sfig1}. Panel (a) shows the case with total winding 2 with a dipolar texture typical for the TE-TM field. The trajectories of the exceptional points are shown with red and blue lines, which are defined by the condition $\Omega'_x=0$ (vertical orientation of the real field). At a large scale, the relative position of the two diabolical points does not play any role and the texture is given by $\theta_r=2\varphi$ (where $\varphi$ is the polar angle of the wave vector). The exceptional points at infinity  are  therefore necessarily located at $\varphi=\pm\pi/4, \pm3\pi/4$ (while initially they move from the diabolical points in the vertical direction).

The second situation is illustrated by Fig.~\ref{Sfig1}(b), where the overall winding of the real field is zero. The angle of the real field at infinity is $\theta_r=0$ and there can be no exceptional points. The condition $\Omega'_x=0$ (vertical real field) gives the trajectories of the exceptional points shown in the figure, with the annihilation points marked with  crosses.

Finally, the winding of the exceptional points themselves is shown with a color of the arrows (red and blue).

The winding of an exceptional point is the winding of the phase of the complex eigenenergies. Indeed, the energy of a second-order exceptional point is given by
\begin{equation}
    E=\pm\sqrt{\alpha q}e^{iw\phi/2}
\end{equation}
where $w=\pm 1$ is precisely the winding number and $q$ is the absolute value of the parameter controlling the deviation from the exceptional point, for example, the wave vector (measured from the exceptional point). This eigenenergy can be obtained as a solution for different Hamiltonians, with different variation of real and imaginary effective fields.
Below, we show that when the imaginary effective field $\bm{\Omega}''$ is constant, the winding of the exceptional point can be easily determined from the texture of the real effective field $\bm{\Omega}'$. The most general form of a Hamiltonian of such type is
\begin{eqnarray}
    \hat{H}&=&\alpha q \left(\sigma_x\cos\varphi+iw\sigma_y\sin\varphi\right)\\
    &+&a(\left(\sigma_x\cos\theta+\sigma_y\sin\theta\right)\nonumber\\
    &+&ia\left(\sigma_x\cos(\theta\pm\pi/2)+\sigma_y\sin(\theta\pm\pi/2)\right)\nonumber
\end{eqnarray}
where the $a$ gives the values of the real and imaginary effective field oriented at an angle $\theta$ and $\theta+\pi/2$ respectively, equal and perpendicular at the exceptional point, $w$ is the winding of the parameter-dependent part of the real effective field and $\alpha$ is its strength. The eigenenergies of this Hamiltonian are given by
\begin{equation}
    E\approx\pm\sqrt{q}\sqrt{2\alpha a}e^{\pm iw\varphi/2}e^{-i\theta/2}
\end{equation}
We see that the winding of an exceptional point is determined by the winding $w$ of the original diabolical point, modified by the relative angle of the real and imaginary fields at each exceptional point $\pm\pi/2$.

\end{document}