###  F.2 Metropolis-Hastings sampler for updates across fluid-types

####  F.2.1 Outline

The index of an unlabeled type $i\in\mathcal{U}$ is selected uniformly at
random. Uniformly choose a subtype $f^{\prime}\in\mathcal{F}_{-i}$ at random,
where $\mathcal{F}_{-i}=\mathcal{F}\setminus y_{i},y_{i}\in\mathcal{F}$. We
then randomly and uniformly select a subtype from the chosen type above. The
potential subtypes include the existing subtypes in fluid-type $f^{\prime}$,
but also include an empty subtype if $|f^{\prime}|<J$. The proposed fluid-type
and subtype for $i$ are accepted or rejected based on the appropriate
Metropolis-Hastings ratio.

####  F.2.2 Details

Let $R\in\Xi^{\,*}_{\mathcal{T},\mathcal{U}}$ be the current state of the
Markov Chain. As remarked above, if $R$ is given, then the corresponding
$\mathcal{V}$ is determined. We choose an unlabeled profile $i\sim
U(\mathcal{U})$ at random, i.e., $i\in\\{N+1,\dots,N+U\\}$. Suppose $i\in
R_{f(k)}$ for some $k\in\\{1,...,K_{f}\\}$. Subsequently, we randomly select
$f^{\prime}\sim U(\mathcal{F}\setminus\\{f\\})$, draw

$k^{\prime}\sim\begin{cases}U\\{1,...,K_{f^{\prime}}+1\\}&\text{for
}K_{f^{\prime}}<J_{f^{\prime}},\\\ U\\{1,...,K_{f^{\prime}}\\}&\text{for
}K_{f^{\prime}}=J_{f^{\prime}}\end{cases}$

and remove $i$ from $R_{f(k)}$ and add it in
$R^{\prime}_{f^{\prime}(k^{\prime})}$. If $k^{\prime}\leq K_{f^{\prime}}$ then
$i$ is added to an existing subtype $R_{f(k^{\prime})}$ of fluid-type
$f^{\prime}$, whereas if $k^{\prime}=K_{f^{\prime}}+1$ then $i$ forms a new
subtype $R_{f(K_{f^{\prime}}+1)}=\\{i\\}$.

If (Case A) $R_{f(k)}=\\{i\\}$, the update removes subtype $R_{f(k)}$ from
fluid-type $f$, so the proposed partition for fluid-type $f$ is
$R^{\prime}_{f}=\cup_{\tilde{k}=1}^{K^{\prime}_{f}}\\{R^{\prime}_{f(\tilde{k})}\\}$,
where

$R^{\prime}_{f(\tilde{k})}=\begin{cases}R_{f(\tilde{k})}&\text{for
}\tilde{k}\in\\{1,...,k-1\\},\\\ R_{f(\tilde{k}+1)}&\text{for
}\tilde{k}\in\\{k,...,K^{\prime}_{f}=K_{f}-1\\}.\end{cases}$

If (Case B) $n_{f(k)}>1$, then $R_{f(k)}$ is not a singleton, so $K_{f}$
remain unchanged and

$R^{\prime}_{f(\tilde{k})}=\begin{cases}R_{f(\tilde{k})}&\text{for
}\tilde{k}\in\\{1,...,K_{f}\\}\setminus k,\\\
R_{f(\tilde{k})}\setminus\\{i\\}&\text{for }\tilde{k}=k.\end{cases}$

For both Cases A and B, if $k^{\prime}=K_{f^{\prime}}+1$, implying
$K_{f^{\prime}}<J_{f^{\prime}}$, then a new subtype is created in the proposed
fluid-type $f^{\prime}$. We obtain
$R^{\prime}_{f^{\prime}}=\cup_{\tilde{k}=1}^{K^{\prime}_{f^{\prime}}}\\{R^{\prime}_{f^{\prime}(\tilde{k})}\\}$
with

$R^{\prime}_{f^{\prime}(\tilde{k})}=\begin{cases}R_{f^{\prime}(\tilde{k})}&\text{for
}\tilde{k}\in\\{1,...,K_{f^{\prime}}\\},\\\ \\{i\\}&\text{for
}\tilde{k}=K^{\prime}_{f^{\prime}}=K_{f^{\prime}}+1.\end{cases}$

On the other hand, if $k^{\prime}\leq K_{f^{\prime}}$, that is, the proposed
subtype assignment already exists, then we have

$R^{\prime}_{f^{\prime}(\tilde{k})}=\begin{cases}R_{f^{\prime}(\tilde{k})}&\text{for
}\tilde{k}\in\\{1,...,K_{f^{\prime}}\\}\setminus k^{\prime},\\\
R_{f^{\prime}(k^{\prime})}\cup\\{i\\}&\text{for
}\tilde{k}=k^{\prime},\end{cases}$

which is permitted for $K_{f^{\prime}}\leq J_{f^{\prime}}$.

Depending on $K_{f^{\prime}}$ only, the proposal probabilities are

$q(R^{\prime}|R)=\begin{cases}(U(F-1)\,(K_{f^{\prime}}+1))^{-1}&\text{if
}K_{f^{\prime}}<J_{f^{\prime}},\\\ (U(F-1)\,K_{f^{\prime}})^{-1}&\text{if
}K_{f^{\prime}}=J_{f^{\prime}}.\end{cases}$

In Case A with $R_{f(k)}=\\{i\\}$, the Hasting’s ratio
$q(R|R^{\prime})/q(R^{\prime}|R)$ is

$H_{A}(R^{\prime}|R)=\begin{cases}\frac{K_{f^{\prime}}+1}{K_{f}}&\text{if
}K_{f^{\prime}}<J_{f^{\prime}},\\\ \frac{K_{f^{\prime}}}{K_{f}}&\text{if
}K_{f^{\prime}}=J_{f^{\prime}}.\end{cases}$

When $n_{f(k)}>1$, for Case B.1 with $K_{f}<J_{f}$, the Hasting’s ratio is

$H_{B.1}(R^{\prime}|R)=\begin{cases}\frac{K_{f^{\prime}}+1}{K_{f}+1}&\text{if
}K_{f^{\prime}}<J_{f^{\prime}},\\\ \frac{K_{f^{\prime}}}{K_{f}+1}&\text{if
}K_{f^{\prime}}=J_{f^{\prime}},\end{cases}$

whereas for Case B.2 with $K_{f}=J_{f}$, the Hastings ratio is simply
$H_{B.2}(R^{\prime}|R)=H_{A}(R^{\prime}|R)$.

###  F.3 Gibbs sampler for updates across fluid-types in the Cut-Model
inference

####  F.3.1 Outline

The Gibbs Sampler presented here is a special case of a Nested MCMC sampler.
Nested MCMC for Cut-Models was suggested in Plummer (2015) and is described
above in Section 6, where some of the following notation is introduced. Our
sampler proposes a new partition by randomly selecting the fluid-type and
subtype for an unlabeled profile $i$ according to the full conditional
probability $p(R^{\prime}|R_{-i})$. Details below explain the sampling
algorithm for classifying a single profile (SPC) and multiple profiles
simultaneously (JPC).

####  F.3.2 Details

After obtaining $\pi(Q|\mathbf{X}_{\mathcal{T}})$ using MCMC and the sampler
in Section F.1, for each $t=1,...,T_{0}$, Gibbs sampling is employed to
simulate a side chain from $\pi(R|\mathbf{X},Q^{(t)})$. For the $t$th side
chain, we first randomly select $R^{(t)}_{0}$, an initial state for the row
partitions of $(\mathcal{T},\mathcal{U})$, conditioned on $Q^{(t)}$ such that
$R^{(t)}_{0}\in\Xi^{\,*}_{\mathcal{T},\mathcal{U}}(Q^{(t)})$. Therefore, the
initial joint partition $R^{(t)}_{0}$ of the labeled and unlabeled data
contains the partition $Q^{(t)}$ of the labeled data as a joint sub-partition.

At step $t^{\prime}=1,\dots,T_{1}$ of this MCMC chain, we randomly select an
unlabeled profile $i\sim U(\mathcal{U})$. Subsequently, we construct a Gibbs
sampler to draw
$R^{(t,t^{\prime})}\in\Xi^{\,*}_{\mathcal{T},\mathcal{U}}(Q^{(t)})$ from the
full conditional distribution

$R^{(t,t^{\prime})}\sim\pi\left(\cdot|R^{(t,t^{\prime}-1)}_{-i}\right),$

which is equivalent to sampling the vector of probabilities with components
given by

$\Pr\left(i\in
R^{(t,t^{\prime})}_{f(k)}|R^{(t,t^{\prime}-1)}_{-i}\right)\propto
p\left(\mathbf{X}|R^{(t,t^{\prime}-1)}_{-i},i\in
R^{(t,t^{\prime})}_{f(k)}\right)\pi\left(R^{(t,t^{\prime})}\right),$

for $f\in\mathcal{F}$ and

$k\in\begin{cases}\\{1,...,K^{(t,t^{\prime})}_{f}+1\\}&\text{if
}K^{(t,t^{\prime})}<J_{f},\\\ \\{1,...,K^{(t,t^{\prime})}_{f}\\}&\text{if
}K^{(t,t^{\prime})}=J_{f},\end{cases}$

where $K^{(t,t^{\prime})}_{f}$ is the number of subtypes in fluid-type
$f\in\mathcal{F}$ in $R^{(t,t^{\prime})}$. When the analyses are performed
with SPC, $T_{1}=1$ is sufficient for all $t\in T_{0}$; in other words, a
single draw by the Gibbs sampler is sufficient for each subchain. However, in
the JPC analyses, for each proposal by the Gibbs sampler, we select an
unlabeled profile $i\sim\text{Uniform}(\mathcal{U})$ at random and draw
$R^{\prime}$ from $p(R^{\prime}|R_{-i})$. In this case, length of each
subchain $T_{1}$ is set such that each unlabeled profile $i\in\mathcal{U}$
would be updated $\geq 10$ times on average.

###  F.4 Runtime, Convergence and Implementation Checks for Metropolis
Hastings

To get an idea of runtimes, in the best case, the Cut-Model setup estimating a
fluid-type for a single profile, we get around 10 Effective Independent
Samples per second (so 600 EIS’/minute) on a standard office desktop computer
in a careful Java implementation, but without parallelization. This value is
computed using the Effective Sample Size (ESS) for the fluid-type of the
unlabeled profile and excludes the one-off cost of sampling subtypes for
training data. In the worst case, Bayesian inference for joint estimation of
fluid-type for our test set of 46 samples, which includes estimating the
subtype structure for the 321 profiles in the training data, we obtain
approximately 1.5 EIS’/minute, equivalently, roughly a half-day to get the
joint fluid-type posterior for 46 profiles if we target an ESS equal 1000 for
each profile. This is just representative, using a central ESS value among
those we observed for sample fluid-type values
$y_{i}=\\{f:i\in\mathcal{V}_{f}\\}$.

For convergence diagnostics, we focus on achieving large ESS values for the
random fluid-type $y_{i}$ of unlabeled samples $i\in\mathcal{U}$. Across all
740 MCMC runs for the different analyses in this paper (the majority being
LOOCV checks on training data) an ESS well above 1000 seems typical. The ESS
values and ranges for all the analyses are given in Table 8.

Table 8: Range of the ESS of the mode type variable across all unlabelled profiles for each analysis. The mode type variable is coded in binary format, where $1=\text{yes}$ and $0=\text{no}.$ The ESS range is calculated excluding any unlabelled profiles with a mode type that has a posterior probability of 1.0. Below SPC and JPC stand for “Single” and “Joint” Profile Classification, respectively. Analysis | ESS Median | Purpose | Section
---|---|---|---
& Range
LOOCV on the training set in Bayesian inference |  4500 (2772, 4500) | Evaluate goodness-of-fit, calibration and | Section 8.3
classification performance with Bayesian
inference, using model and lab-quality data.
LOOCV on the training set in Cut-Model inference |  4500 (2412, 4500) | Evaluate goodness-of-fit, calibration and
classification performance with Cut-Model
inference, using model and lab-quality data,
and compare to the corresponding Bayesian.
analysis.
JPC of profiles in the test set in Bayesian inference |  3955 (1672, 4500) | Evaluate the performance of the trained model | Section 8.4
for jointly classifying casework-like data using
Bayesian inference and lab-quality training data.
JPC of profiles in the test set in Cut-Model inference |  4500 (2402, 4500) | Evaluate the performance of the trained model
for jointly classifying casework-like data using
Cut-Model inference and lab-quality training
data, and compare to the corresponding
Bayesian analysis.
SPC of profiles in the test set in Bayesian inference |  4500 (1313, 4500) | Evaluate the performance of the trained model
for classifying casework-like data one at a time
using Bayesian inference and lab-quality training
data, and compare to the corresponding joint
analysis.
SPC of profiles in the test set in Cut-Model inference |  4500 (936, 4500) | Evaluate the performance of the trained model
for classifying casework-like data one at a time
using Cut-Model inference, and lab-quality
training data, and compare to the corresponding
Bayesian analysis.
JPC of profiles in the test set in Bayesian inference ($J_{f}=1$ & $L_{g}=1$) |  4500 (3117, 4500) |  To investigate the effect of cluster bounds on the classification of casework-like data. | Appendix I
JPC of profiles in the test set in Bayesian inference ($J_{f}=1$ & $L_{g}=M_{g}$) |  4500 (3496, 4500)
JPC of profiles in the test set in Bayesian inference ($J_{f}=10$ & $L_{g}=M_{g}$) |  3703 (1856, 4500)
JPC of profiles in the test set in Bayesian inference ($J_{f}=15$ & $L_{g}=M_{g}$) |  3765 (1277, 4500)

The implementation is checked by verifying our software reproduced the MCRP
distributions when setting the log-likelihood function to return zero. We
hand-checked likelihood calculations for a small number of profiles in a
single subtype. Also, we performed classification on the unlabeled profile
with all missing values and checked that the posterior distribution for its
fluid-type was indeed uniform on $\mathcal{F}$. This is discussed further in
Appendix G. Our LOOCV checks on the training data check correct (but blinded)
classification of training data with known labels.

##  Appendix G Missing data and unbiased fluid-type labelling

There are no missing marker values $x_{i,j},\ i\in\mathcal{N},\
j\in\mathcal{M}$ in the training and test data. However, if a profile did have
missing marker values, then this is easy to handle as explained in Section
5.3.3: in the product over ${(i,j)\in C(f(k),g(l))}$ in Equation 8, we simply
omit from the product any cells $(i,j)$ for which the value of $x_{i,j}$ is
missing.

The prior $\pi_{\mathcal{V}}(\mathcal{V})=F^{-U}$ on the fluid-type assignment
partition $\mathcal{V}$ of unlabeled profiles, presented in Section 5.3.2,
takes the prior for the fluid-type of each unlabeled profile to be independent
of the fluid-types of other unlabeled profiles and uniform on $\mathcal{F}$.
However, it may be of concern that the labeled training data somehow distort
this and changes the prior weighting, perhaps because the number of training
profiles varies from one fluid-type to another, and the number of subtypes is
upper-bounded. We formalize this by considering the posterior distribution of
$\mathcal{V}$ when all the data $x_{i},\ i\in\mathcal{U}$ are missing: when
$\mathbf{X}_{\mathcal{U}}$ provides no information about the fluid-types of
unlabeled profiles we require
$\pi(\mathcal{V}|X_{\mathcal{T}},\mathbf{X}_{\mathcal{U}})=F^{-U}$.

A joint partition $R\in\Xi^{\,*}_{\mathcal{T},\mathcal{U}}$ of all profiles
specifies the assignment $\mathcal{V}$ of unlabeled profiles to fluid-types.
However, it contains more information as each unlabeled profile
$i\in\mathcal{V}_{f}$ is assigned some subtype $k\in\\{1,\dots,K_{f}\\}$, so
$i\in R_{f(k)}$. Let $Q(R)$ be the joint partition on training profiles,
defined in Section 6, which we get by removing all the floating profiles from
$R$.

Consider the likelihood $p(\mathbf{X}|R)$ in Equation 5.3.2 with
$\mathbf{X}=(\mathbf{X}_{\mathcal{T}},\mathbf{X}_{U})$. In this case, with all
$\mathbf{X}_{\mathcal{U}}$ missing, the entries in $\mathbf{X}_{\mathcal{U}}$
do not contribute to $c$ and $s$ in Equation 9. Consequently, the likelihood
is the probability for $\mathbf{X}_{\mathcal{T}}$ only, and this depends only
on the joint partition $Q$ of the labeled profiles. It follows that

$p(\mathbf{X}|R)=p(\mathbf{X}_{\mathcal{T}}|Q),$

where $Q=Q(R)$ as in a Cut-Model.

The outcome $\\{Q,\mathcal{V}\\}$ occurs iff
$R\in\mathbf{X}^{\,*}_{\mathcal{T},\mathcal{V}}(Q)$ occurs, where
$\mathbf{X}^{\,*}_{\mathcal{T},\mathcal{V}}(Q)$, defined in Equation 17, is
the set of joint partitions $R$ that contain $Q$ and assign unlabeled profiles
to fluid-types according to $\mathcal{V}$. We have

$\displaystyle\pi(Q,\mathcal{V}|\mathbf{X})$
$\displaystyle\propto\sum_{R\in\mathbf{X}^{\,*}_{\mathcal{T},\mathcal{V}}(Q)}p(\mathbf{X}|R)\pi_{R}(R)$
$\displaystyle=p(\mathbf{X}_{\mathcal{T}}|Q)\sum_{R\in\mathbf{X}^{\,*}_{\mathcal{T},\mathcal{V}}(Q)}\pi_{R}(R)$
$\displaystyle=p(\mathbf{X}_{\mathcal{T}}|Q)\pi_{Q,\mathcal{V}}(Q,\mathcal{V}),$

with
$\pi_{Q,\mathcal{V}}(Q,\mathcal{V})=\pi_{Q}(Q)\pi_{\mathcal{V}}(\mathcal{V})$,
so $\pi(\mathcal{V}|\mathbf{X})=\pi_{\mathcal{V}}(\mathcal{V})=F^{-U}$. This
result is intuitive: no data, no classification. The same property,
$\pi_{cut}(\mathcal{V}|\mathbf{X})=\pi_{\mathcal{V}}(\mathcal{V})$, holds in
the Cut model case. We used this property to check our code.

The acceptance probabilities for move type $T\in\\{A,B.1,B.2\\}$ are

$\displaystyle\alpha_{T}(R^{\prime}|R)$ $\displaystyle=\min\left\\{1,\
\frac{\pi(R^{\prime}|X)}{\pi(R|X)}H_{T}(R^{\prime}|R)\right\\}.$

There is quite a bit of cancellation in the acceptance probability. Our
implementation always works with the full sum of log-likelihoods across a
fluid-type, but in this sum, only the log-likelihoods of subtypes that change
in the update are themselves updated. This is the time-limiting step and is
handled efficiently.

##  Appendix H Results

In section 8.3, Figure 3 summarises the results of the LOOCV of the Bayesian
analysis on the training set. For LOOCV of the Cut-Model analysis on the
training dataset, the summary in the same style is presented in Figure 8,
which conveys a very similar message as Figure 3, in particular, the confusion
tables (bottom right panel) are identical in both figures. The choice of the
inference framework, Bayes or Cut-Model, does not influence the results of
interest or our conclusion on the classification.

Figure 8: Posterior mode probability and mode of the correct type (for mRNA
profiles in the training set) estimated from the LOOCV analysis on the
training data with Cut-Model inference. This is discussed in Section 8.3.

Confusion tables for the classification of the test data using the fluid-type
of the posterior mode in SPC and JPC analyses (identical for Bayes and Cut-
Model analyses) are given in Table 9. The only differences are observed in the
classification of two mRNA profiles: a CVF sample with an inconclusive
$\log_{10}\text{BF}$ for a CVF label $<1$ across all four analyses; an SMN
profile with $\log_{10}\text{BF}$ for an SMN label equal 1.054 for Bayesian
JPC, but less than one in the other three analyses, so right on the
borderline.

Table 9: The posterior mode estimated across the samples of the test set by SPC and JPC (identical for both Bayesian- and Cut-Model inference). | Posterior Mode
---|---
True Type | One at a Time (SPC) | Joint (JPC)
| CVF | MTB | SLV | BLD | SMN | CVF | MTB | SLV | BLD | SMN
Cervical fluid | 23 | 1 | 0 | 0 | 0 | 22 | 2 | 0 | 0 | 0
Menstrual blood | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0
Saliva | 0 | 0 | 10 | 0 | 0 | 0 | 0 | 10 | 0 | 0
Blood | 0 | 0 | 0 | 2 | 0 | 0 | 0 | 0 | 2 | 0
Semen | 1 | 0 | 0 | 0 | 9 | 0 | 0 | 0 | 0 | 10

##  Appendix I Sensitivity analysis for bounds on partition count.

Figure 9 presents the comparisons in the estimated posterior probabilities of
the correct types for the mRNA profiles in the test set between the default
model ($J_{f}=5,L_{g}=M_{g}$) and other choices of $J_{f}$. Given a
$(J_{f},L_{g})$ combination, the results are very similar between Bayesian and
Cut-Model inferences.

From panels (a)–(d) of Figure 9, there appears to be greater dispersion in the
posterior probabilities of the true types for $(J_{f}=1,L_{g}=1)$ and
$(J_{f}=1,L_{g}=M_{g})$ than the default model. The dispersion appears to be
strong for $(J_{f}=1,L_{g}=1)$. This suggests that the posterior probabilities
are poorly calibrated when we do not accommodate the heterogeneity within a
fluid-type and the variation in the marker amplification within a group
(within a fluid-type). Panels (e)–(h) of Figure 9 show that there is generally
close agreement in the posterior probabilities of the true types between the
default model and $(J_{f}=10,L_{g}=M_{g})$, and between the default model and
$(J_{f}=15,L_{g}=M_{g})$. These results do not depend on the choice of the
inference framework.

Figure 9: Comparing the estimated posterior probabilities of the correct types
of the mRNA profiles in the test set between the default model
($J_{f}=5,L_{g}=M_{g}\,\,\forall f,g\in\mathcal{F}$) and other choices of
$J_{f}$ and $L_{g}$. Each panel in the left columns is a comparison between
two analyses in the Bayesian framework, while each in the right is in the Cut-
Model framework. For legibility, the positions on axes are in logit-space,
while the labels represent the corresponding probability values.  Figure 10:
Posterior distribution of the number of subtypes within each fluid-type. Each
panel in the left columns is a comparison between two analyses in the Bayesian
framework, while each in the right is in the Cut-Model framework.


Without loss of generality, we may assume $\textup{sgn}(\theta^{u}_{k})=1.$
Then,

$\begin{split}\theta^{ver}_{2,j+k}&=\theta^{u}_{k}-\theta^{l}_{k}-\varepsilon\textup{curl}^{d,\varepsilon}((\theta^{hor}_{0,j+k},\theta^{ver}_{0,j+k}))-\varepsilon\textup{curl}^{d,\varepsilon}((\theta^{hor}_{1,j+k},\theta^{ver}_{1,j+k}))\\\
&\leq\frac{14}{6}\pi-\varepsilon\textup{curl}^{d,\varepsilon}((\theta^{hor}_{0,j+k},\theta^{ver}_{0,j+k}))-\varepsilon\textup{curl}^{d,\varepsilon}((\theta^{hor}_{1,j+k},\theta^{ver}_{1,j+k}))\end{split}$
(5.31)

and

$\begin{split}\theta^{ver}_{2,j+k}&=\theta^{u}_{k}-\theta^{l}_{k}-\varepsilon\textup{curl}^{d,\varepsilon}((\theta^{hor}_{0,j+k},\theta^{ver}_{0,j+k}))-\varepsilon\textup{curl}^{d,\varepsilon}((\theta^{hor}_{1,j+k},\theta^{ver}_{1,j+k}))\\\
&\geq\frac{10}{6}\pi-\varepsilon\textup{curl}^{d,\varepsilon}((\theta^{hor}_{0,j+k},\theta^{ver}_{0,j+k}))-\varepsilon\textup{curl}^{d,\varepsilon}((\theta^{hor}_{1,j+k},\theta^{ver}_{1,j+k})).\end{split}$
(5.32)

As $\theta^{ver}_{2,j+k}\in[-\pi,\pi)$ and
$\varepsilon\textup{curl}^{d,\varepsilon}((\theta^{hor}_{0,j+k},\theta^{ver}_{0,j+k})),\varepsilon\textup{curl}^{d,\varepsilon}((\theta^{hor}_{1,j+k},\theta^{ver}_{1,j+k}))\in\\{0,\pm
2\pi\\}$ it follows that
$\varepsilon\textup{curl}^{d,\varepsilon}((\theta^{hor}_{0,j+k},\theta^{ver}_{0,j+k}))+\varepsilon\textup{curl}^{d,\varepsilon}((\theta^{hor}_{1,j+k},\theta^{ver}_{1,j+k}))=2\pi$.
Thus, it holds $|\theta^{ver}_{2,j+k}|\leq\pi/3$. Then we estimate using
(5.30)

$\begin{split}&E_{\varepsilon,\delta}(u,[0,2\varepsilon]\times[j\varepsilon,(j+2)\varepsilon])\\\
&\geq\frac{\varepsilon^{2}}{2}\left(4\cos^{2}\left(\frac{\theta^{ver}_{1,j}+\theta^{ver}_{1,j+1}}{2}\right)-4(1-\delta)\cos(\theta^{ver}_{1,j})-4(1-\delta)\cos(\theta^{ver}_{1,j+1})+4(1-\delta)^{2}\right)\\\
&\geq\frac{\varepsilon^{2}}{2}\left(4\cos^{2}\left(\frac{\pi}{3}\right)-8(1-\delta)\right)\geq\frac{\varepsilon^{2}}{2}\left(1-8(1-\delta)\right)\geq\frac{1}{16}\varepsilon^{2},\end{split}$
(5.33)

where we used $\cos(\pi/3)=1/2.$ This shows (5.26). Eventually, we sum
estimate (5.26) to find

$\displaystyle E_{\varepsilon,\delta}(u,[0,2\varepsilon]\times[0,1])$
$\displaystyle\geq\frac{1}{3}\sum_{j=0}^{\left\lfloor
1/\varepsilon\right\rfloor-2}E_{\varepsilon,\delta}(u,[0,2\varepsilon]\times[j\varepsilon,j\varepsilon+2\varepsilon])$
$\displaystyle\geq\left(\left\lfloor\frac{1}{\varepsilon}\right\rfloor-1\right)\frac{\varepsilon^{2}}{48}\geq\frac{\varepsilon}{192}.$

∎

We continue by discussing the case of $\delta$ being bounded away from one and
away from zero.

######  Lemma 5.5.

There is a constant $c>0$ such that for all
$(\varepsilon,\delta)\in(0,1/2)\times(0,63/64)$ and
$u\in\mathcal{S}_{0}(\Omega,\varepsilon)$ it holds

$E_{\varepsilon,\delta}(u,[0,2\varepsilon]\times[0,1])\geq
c\varepsilon\delta^{2}.$

###### Proof.

Let $\eta:=1/(32\cdot 64^{2})$ and assume that
$E_{\varepsilon,\delta}(u,\Omega)\leq\eta\varepsilon$ (otherwise, there is
nothing to show). We introduce the set of vortices along the left boundary

$\begin{split}V_{0}\coloneqq&\left\\{(0,j)\in\mathbb{Z}^{2}\colon
j\varepsilon\in(0,1-\varepsilon)\cap\varepsilon\mathbb{Z}\textup{ and
}\left|\varepsilon\textup{curl}^{d,\varepsilon}((\theta^{hor}_{0,j},\theta^{ver}_{0,j}))\right|=2\pi\right\\}\\\
\subseteq&\left\\{(0,j)\in\mathbb{Z}^{2}\colon
j\varepsilon\in(0,1-\varepsilon)\cap\varepsilon\mathbb{Z}\textup{ and
}\max\\{|\theta^{hor}_{0,j}|,|\theta^{hor}_{0,j+1}|,|\theta^{ver}_{1,j}|\\}\geq\frac{2}{3}\pi\right\\}=:B_{0}.\end{split}$
(5.34)

Lemma 3.5 and Remark 3.7 for $\beta=2\pi/3$ yield the estimate

$\mathcal{H}^{0}(V)\leq\mathcal{H}^{0}(B_{0})\leq\frac{8E_{\varepsilon,\delta}(u,\Omega)}{\varepsilon^{2}(1-\delta)^{2}}\leq\frac{8\eta\varepsilon}{\varepsilon^{2}(1-\delta)^{2}}\leq\frac{1}{4\varepsilon}.$
(5.35)

Let $\mathcal{I}_{0}$ be the set of cells on the left boundary without
vortices i.e.,

$\mathcal{I}_{0}\coloneqq\\{(0,j)\in\mathbb{Z}^{2}\colon
j\varepsilon\in(0,1)\textup{ and }j\notin B_{0}\\}.$

By (5.35) and $\varepsilon\in(0,1/2)$, we have

$\mathcal{H}^{0}(\mathcal{I}_{0})\geq\frac{3}{4\varepsilon}-1\geq\frac{1}{4\varepsilon}.$
(5.36)

Now, let $(0,j)\in\mathcal{I}_{0}$. Since $V_{0}\subseteq B_{0}$, it holds

$0=\theta^{hor}_{0,j+1}-\theta^{ver}_{1,j}-\theta^{hor}_{0,j}=:\theta_{1}^{j}-\theta_{2}^{j}-\theta_{3}^{j}.$
(5.37)

We claim that there exists $k\in\\{1,2,3\\}$ such that

$|\theta_{k}^{j}|\notin\left[\frac{3}{4}\theta_{opt},\frac{5}{4}\theta_{opt}\right].$
(5.38)

Assuming (5.38) is true for every $j\in\mathcal{I}_{0}$, we estimate the
energy as in (3) in the proof of Lemma 3.5

$\begin{split}&E_{\varepsilon,\delta}(u,\Omega)\\\ \geq&\frac{1}{2\cdot
64^{2}}\varepsilon^{2}\sum_{(0,j)\in\mathcal{I}_{0}}\left[(1-\delta-\cos(\theta^{hor}_{0,j}))^{2}+(1-\delta-\cos(\theta^{hor}_{0,j+1}))^{2}+(1-\delta-\cos(\theta^{ver}_{1,j}))^{2}\right]\\\
\geq&\frac{1}{2\cdot
64^{2}}\varepsilon^{2}\mathcal{H}^{0}(\mathcal{I}_{0})\min\left\\{\left|\cos(\theta_{opt})-\cos\left(\frac{3}{4}\theta_{opt}\right)\right|^{2},\left|\cos(\theta_{opt})-\cos\left(\frac{5}{4}\theta_{opt}\right)\right|^{2}\right\\}\\\
\geq&\frac{\varepsilon}{8\cdot
64^{2}}\min\left\\{\left|\int_{\frac{3}{4}\theta_{opt}}^{\theta_{opt}}\sin(s)\textup{
d}s\right|^{2},\left|\int_{\theta_{opt}}^{\frac{5}{4}\theta_{opt}}\sin(s)\textup{
d}s\right|^{2}\right\\}\\\
\geq&\sin^{2}\left(\frac{3}{4}\theta_{opt}\right)\frac{\varepsilon\theta_{opt}^{2}}{8\cdot
16\cdot 64^{2}}\\\
\geq&\frac{1}{4}\left(\frac{3}{4}\right)^{2}\theta_{opt}^{2}\frac{\varepsilon\delta}{8\cdot
16\cdot 64^{2}}\geq\frac{9}{4\cdot 8\cdot 16^{2}\cdot
64^{2}}\varepsilon\delta^{2},\end{split}$ (5.39)

where we used $\delta^{1/2}\leq\theta_{opt}=\arccos(1-\delta)$. This yields
the lower bound with $c\coloneqq 9/(4\cdot 8\cdot 16^{2}\cdot 64^{2})$.
Hence, it is left to prove that there is $k\in\\{1,\dots,3\\}$ such that
(5.38) holds. For the sake of a contradiction, assume
$|\theta^{j}_{k}|\in[3\theta_{opt}/4,5\theta_{opt}/4]$ for all
$k\in\\{1,\dots,3\\}$. Without loss of generality, we may assume
$\theta_{1}^{j}\geq 0$. Since $(\theta^{hor},\theta^{ver})$ is curl-free, we
have $|\theta_{1}^{j}|=\theta_{1}^{j}=\theta_{2}^{j}+\theta_{3}^{j}$. If
$\textup{sgn}(\theta_{2}^{j})=\textup{sgn}(\theta_{3}^{j})$, we estimate

$\frac{5}{4}\theta_{opt}\geq\theta_{1}^{j}\geq\frac{6}{4}\theta_{opt}>\frac{5}{4}\theta_{opt},$

which is a contradiction. Next, assume
$\textup{sgn}(\theta_{2}^{j})\neq\textup{sgn}(\theta_{3}^{j}).$ Then

$\frac{3}{4}\theta_{opt}\leq\theta_{1}^{j}=\theta_{2}^{j}+\theta_{3}^{j}\leq-\frac{3}{4}\theta_{opt}+\frac{5}{4}\theta_{opt}=\frac{1}{2}\theta_{opt}<\frac{3}{4}\theta_{opt},$

which is again a contradiction. ∎

####  5.2.2 Lower bound for small $\delta\in(0,\overline{\delta})$

In this section we prove the lower bound for small interaction parameters
$\delta\in(0,\overline{\delta})$. In this setting Lemma 3.4 allows us to
exclude the presence of vortices. The remainder of the proof is inspired by
the strategy from [GZ23, Lemma 4] for a related continuum model, see also
[GZ24, GRTZ24].

######  Lemma 5.6.

Let $\eta:=\frac{1}{2\cdot 4\cdot 12^{2}\cdot 256}$. Then there is a constant
$c_{l}>0$ such that for all $\delta\in(0,\eta^{2}/16)$,
$\varepsilon\in(0,\eta/4)$ with $\delta\leq\varepsilon$ and
$u\in\mathcal{S}_{0}(\Omega,\varepsilon)$ it holds

$E_{\varepsilon,\delta}(u,\Omega)\geq
c_{l}\min\left\\{\delta^{2},\varepsilon\delta^{3/2}\left(\left|\ln\left(\frac{\varepsilon}{\delta^{1/2}}\right)\right|+1\right)\right\\}.$

######  Remark 5.7.

Note, that

$\min\left\\{\delta^{2},\varepsilon\delta^{3/2}\left(\left|\ln\left(\frac{\varepsilon}{\delta^{1/2}}\right)\right|+1\right)\right\\}=\begin{cases}\delta^{2}\quad&\textup{
if }\delta^{1/2}\leq\varepsilon,\\\
\varepsilon\delta^{3/2}\left(\left|\ln\left(\frac{\varepsilon}{\delta^{1/2}}\right)\right|+1\right)\quad&\textup{
if }\varepsilon\leq\delta^{1/2}.\end{cases}$

###### Proof.

Let $u\in\mathcal{S}_{0}(\Omega,\varepsilon)$, $c_{1}\coloneqq 1/(32\cdot
12^{2}\cdot 256)$ and $c_{2}\coloneqq\eta$.

Step 1: Preparation. We may assume that $u$ satisfies

$E_{\varepsilon,\delta}(u,\Omega)\leq
c_{1}\min\left\\{\delta^{2},\varepsilon\delta^{3/2}\left(\left|\ln\left(\frac{\varepsilon}{\delta^{1/2}}\right)\right|+1\right)\right\\}.$
(5.40)

Otherwise there is nothing to show. Next, using (2.26) and Lemma 3.4 yields

$\begin{split}E_{\varepsilon,\delta}(u,\Omega)\geq&\varepsilon^{2}\sum_{(i,j)\in\mathcal{I}_{\textup{NN}}^{\textup{ver}}}(1-\delta-\cos(\theta^{ver}_{i,j}))^{2}+2\varepsilon^{2}\sum_{(i,j)\in\mathcal{I}_{\textup{NN}}^{\textup{ver}}}p(\theta^{ver}_{i,j+1},\theta^{ver}_{i,j})\left|\sin\left(\frac{\theta^{ver}_{i+1,j}}{2}\right)-\sin\left(\frac{\theta^{ver}_{i,j}}{2}\right)\right|^{2}\\\
\geq&\varepsilon^{2}\sum_{(i,j)\in\mathcal{I}_{\textup{NN}}^{\textup{ver}}}(1-\delta-\cos(\theta^{ver}_{i,j}))^{2}+\varepsilon^{2}\sum_{(i,j)\in\mathcal{I}_{\textup{NN}}^{\textup{ver}}}\left|\sin\left(\frac{\theta^{ver}_{i,j+1}}{2}\right)-\sin\left(\frac{\theta^{ver}_{i,j}}{2}\right)\right|^{2}\\\
&-2\pi^{2}(1-2c_{p})\varepsilon^{2}\mathcal{H}^{0}(A^{ver}),\end{split}$
(5.41)

where
$A^{ver}=\\{(i,j)\in\mathcal{I}_{\textup{NN}}^{\textup{ver}}\colon\max\\{|\theta^{ver}_{i,j}|,|\theta^{ver}_{i,j+1}|\\}\geq\pi/8\\}$
and
$(\theta^{hor},\theta^{ver})\colon\Omega\cap\varepsilon\mathbb{Z}^{2}\rightarrow[-\pi,\pi)^{2}$
is the associated angular velocity field of $u$. Analogously, the same is true
for the horizontal energy contribution, replacing $A^{ver}$ by the
corresponding set $A^{hor}.$

Step 2: Reduction to estimate on slices. We claim that for all
$i\in\mathbb{N}$ satisfying

$c_{2}\geq
i\varepsilon\geq\frac{c_{2}}{4}\min\left\\{1,\frac{\varepsilon}{\delta^{1/2}}\left(\left|\ln\left(\frac{\varepsilon}{\delta^{1/2}}\right)\right|+1\right)\right\\}=:a(\varepsilon,\delta)$

one of the following holds

$\begin{split}\textup{(i)}&\
\varepsilon\sum_{j\in\mathcal{J}_{N}}(1-\delta-\cos(\theta^{ver}_{i,j}))^{2}\geq
c_{2}\delta^{2}\\\ \textup{(ii)}&\
\varepsilon\sum_{j\in\mathcal{J}_{N}}\max\\{|1-\delta-\cos(\theta^{ver}_{i,j})|,|1-\delta-\cos(\theta^{ver}_{i,j+1})|\\}\left|\sin\left(\frac{\theta^{ver}_{i,j+1}}{2}\right)-\sin\left(\frac{\theta^{ver}_{i,j}}{2}\right)\right|\geq\frac{c_{2}}{i\varepsilon}\varepsilon\delta^{3/2},\end{split}$

where $\mathcal{J}_{N}\coloneqq[0,1)\cap\varepsilon\mathbb{Z}$. We will show
first how to conclude from the claim. The rest of the proof will then be
devoted to proving the claim.

By Young’s inequality, (i) and (ii) from above, and (5.41), we have

$E_{\varepsilon,\delta}(u,\Omega)\geq\varepsilon\sum_{i\varepsilon\geq
a(\varepsilon,\delta)}\min\left\\{c_{2}\delta^{2},\frac{c_{2}}{i\varepsilon}\varepsilon\delta^{3/2}\right\\}-2\pi^{2}(1-2c_{p})\varepsilon^{2}\mathcal{H}^{0}(A^{ver}).$
(5.42)

By Lemma 3.5 we estimate

$\mathcal{H}^{0}(A^{ver})\leq\frac{2E_{\varepsilon,\delta}(u,\Omega)}{(1-\delta)^{2}\varepsilon^{2}(1-\delta-\cos(\frac{\pi}{8}))^{2}}\leq\frac{4\cdot
20^{2}c_{1}s(\varepsilon,\delta)}{\varepsilon^{2}}\leq 4\cdot
20^{2}c_{1}\left(\frac{\delta}{\varepsilon}\right)^{2}\leq 4\cdot
20^{2}c_{1}\leq\frac{1}{2},$

where we used
$1-\delta-\cos(\pi/8)=\cos(\theta_{opt})-\cos(\pi/8)\geq\cos(\pi/16)-\cos(\pi/8)\geq
1/20$, and $\delta\leq\varepsilon$. Consequently,
$\mathcal{H}^{0}(A^{ver})=0$. Similarly, one obtains
$\mathcal{H}^{0}(A^{hor})=0$. Now, we discuss the cases
$\delta^{1/2}\leq\varepsilon$ and $\delta\leq\varepsilon\leq\delta^{1/2}$
separately.

_Step 2.1: $\delta^{1/2}\leq\varepsilon$._ In this case, inequality (5.42)
reduces to

$E_{\varepsilon,\delta}(u,\Omega)\geq\varepsilon\sum_{c_{2}\geq
i\varepsilon\geq
a(\varepsilon,\delta)}\min\left\\{c_{2}\delta^{2},\frac{c_{2}}{i\varepsilon}\varepsilon\delta^{3/2}\right\\}\geq\frac{c_{2}}{2}\delta^{2}\left(c_{2}-\frac{c_{2}}{4}-2\varepsilon\right)\geq\frac{c_{2}^{2}}{8}\delta^{2},$
(5.43)

where we used $a(\varepsilon,\delta)=c_{2}/4$, $c_{2}\leq 1$, and
$\varepsilon\in(0,c_{2}/4).$

_Step 2.2: $\delta\leq\varepsilon\leq\delta^{1/2}$._ By Remark 5.7 and the
assumption $\delta^{1/2}\leq\eta/4=c_{2}/4$, we have

$a(\varepsilon,\delta)=\frac{c_{2}}{4}s(\varepsilon,\delta)=\frac{c_{2}}{4}\frac{\varepsilon}{\delta^{1/2}}\left(\ln\left(\frac{\delta^{1/2}}{\varepsilon}\right)+1\right)\geq\frac{c_{2}}{4}\frac{\varepsilon}{\delta^{1/2}}\geq\varepsilon.$
(5.44)

It follows that

$\displaystyle E_{\varepsilon,\delta}(u,\Omega)\geq$
$\displaystyle\frac{\varepsilon}{2}\sum_{c_{2}\geq i\varepsilon\geq
a(\varepsilon,\delta)}c_{2}\min\left\\{\delta^{2},\frac{\varepsilon\delta^{3/2}}{i\varepsilon}\right\\}$
$\displaystyle\geq$ $\displaystyle\frac{\varepsilon}{2}\sum_{c_{2}\geq
i\varepsilon\geq
a(\varepsilon,\delta)}c_{2}\frac{\varepsilon\delta^{3/2}}{i\varepsilon}$
$\displaystyle\geq$ $\displaystyle
c_{2}\frac{\varepsilon\delta^{3/2}}{4}\int_{a(\varepsilon,\delta)+\varepsilon}^{c_{2}}\frac{1}{t}\,dt$
$\displaystyle=$ $\displaystyle
c_{2}\frac{\varepsilon\delta^{3/2}}{4}\left(\ln(c_{2})-\ln(a(\varepsilon,\delta)+\varepsilon)\right)$
$\displaystyle\geq$ $\displaystyle
c_{2}\frac{\varepsilon\delta^{3/2}}{4}\left(\ln(c_{2})-\ln(2a(\varepsilon,\delta))\right)$
$\displaystyle=$ $\displaystyle
c_{2}\frac{\varepsilon\delta^{3/2}}{4}\left(\ln(2)-\ln\left(\frac{\varepsilon}{\delta^{1/2}}(\ln(\delta^{1/2}/\varepsilon)+1)\right)\right)$
$\displaystyle=$
$\displaystyle\frac{c_{2}}{4}\varepsilon\delta^{3/2}\left(\ln\left(2\right)-\ln\left(\frac{\varepsilon}{\delta^{1/2}}\left(\ln\left(\frac{\delta^{1/2}}{\varepsilon}\right)+1\right)\right)\right)\geq\frac{c_{2}}{4\cdot
50}\varepsilon\delta^{3/2}\left(\left|\ln\left(\frac{\varepsilon}{\delta^{1/2}}\right)\right|+1\right),$

where we used $\ln(|\ln\sigma|+1)\leq\max\\{24\ln(2)/25,3|\ln\sigma|/4\\}$ and
$\ln(2)\geq 1/2$ in the last estimate.

_Step 3: Proof of the claim._ Assume that (i) and (ii) do not hold for
$i\varepsilon$. Define the sets

$\begin{split}\mathcal{M}_{+}&\coloneqq\left\\{j\varepsilon\in(0,1)\cap\varepsilon\mathbb{Z}\colon\theta^{ver}_{i,j}\in\left(\frac{1}{2}\theta_{opt},\frac{3}{2}\theta_{opt}\right)\right\\}\text{\qquad
and}\\\
\mathcal{M}_{-}&\coloneqq\left\\{j\varepsilon\in(0,1)\cap\varepsilon\mathbb{Z}\colon\theta^{ver}_{i,j}\in\left(-\frac{3}{2}\theta_{opt},-\frac{1}{2}\theta_{opt}\right)\right\\}.\end{split}$
(5.45)

Note that for $0\leq\theta<\pi/3$ it holds

$|\arccos(1-\delta)^{2}-\theta^{2}|=\left|\int_{\arccos(1-\delta)}^{\theta}\frac{1}{2}s\,ds\right|\leq\left|\int_{\arccos(1-\delta)}^{\theta}\sin(s)\,ds\right|=|1-\delta-\cos(\theta)|.$
(5.46)

Similarly, the same estimate holds for $-\pi/3<\theta\leq 0$. Together with
$\mathcal{H}^{0}(A^{ver})=0$ and the assumption that (i) does not hold it
follows for the set
$\mathcal{M}^{c}\coloneqq((0,1)\cap\varepsilon\mathbb{Z})\setminus(\mathcal{M}_{+}\cup\mathcal{M}_{-})$
that

$\begin{split}\varepsilon\mathcal{H}^{0}(\mathcal{M}^{c})=&\varepsilon\sum_{j\varepsilon\in\mathcal{M}_{c}}1\leq\frac{16}{9\theta_{opt}^{4}}\varepsilon\sum_{j\varepsilon\in\mathcal{M}^{c}}(\theta_{opt}^{2}-(\theta^{ver}_{i,j})^{2})^{2}\\\
&\leq\frac{4\cdot
16}{9\theta_{opt}^{4}}\varepsilon\sum_{j\varepsilon\in\mathcal{M}^{c}}\left(1-\delta-\cos(\theta^{ver}_{i,j})\right)^{2}\leq\frac{4\cdot
16c_{2}\delta^{2}}{9\theta_{opt}^{4}}\leq\frac{64}{9}c_{2}\leq\frac{1}{3}.\end{split}$
(5.47)

In particular, since $\varepsilon\in(0,1/2)$

$\varepsilon\mathcal{H}^{0}(\mathcal{M}_{+}\cup\mathcal{M}_{-})\geq\frac{2}{3}-\varepsilon\geq\frac{1}{6}.$
(5.48)

Without loss of generality, we may assume
$\varepsilon\mathcal{H}^{0}(\mathcal{M}_{+})\geq 1/12.$ We introduce the
piecewise affine interpolation of $\theta^{ver}_{i,\cdot}$ on $(0,1)$ by (see
Figure 14)

$\theta_{i}(y)\coloneqq\frac{y-j\varepsilon}{\varepsilon}\theta^{ver}_{i,j+1}+\frac{(j+1)\varepsilon-y}{\varepsilon}\theta^{ver}_{i,j}\quad\textup{
for }y\in[j\varepsilon,(j+1)\varepsilon).$

Note, that $\theta_{i}(\cdot)$ is Lipschitz-continuous on $(0,1)$ and
$|\theta_{i}|\leq\pi/3$ since $\mathcal{H}^{0}(A^{ver})=0$.

Figure 14: Sketch of the function $\theta_{i}$. The black dots are the values
of $\theta^{ver}_{i,\cdot}$, the red dots mark the discrete boundary points of
$\mathcal{M}_{+}$, and the blue rectangles indicate where the interpolated
function $\theta_{i}$ leaves the interval $(\theta_{opt}/2,3\theta_{opt}/2)$.

Furthermore, we introduce the function

$F(\theta)\coloneqq\int_{-\frac{\pi}{3}}^{\theta}|\theta_{opt}^{2}-\gamma^{2}|\textup{
d}\gamma\quad\textup{ for }\theta\in[-\frac{\pi}{3},\frac{\pi}{3}].$

By the Lipschitz continuity of $F\circ\theta_{i}$, the coarea formula (see
[EG92, Theorem 1 in Section 3.4.2]), estimate (5.46),
$|s-t|/2\leq|\sin(s)-\sin(t)|$ for all $|s|,|t|\leq\pi/3$ and the fact that we
assume that (ii) does not hold, we estimate

$\begin{split}&\int_{\mathbb{R}}\mathcal{H}^{0}\left(\partial\left\\{F\circ\theta_{i}>s\right\\}\right)+\mathcal{H}^{0}(\partial\\{F\circ\theta_{i}<s\\})\textup{
d}s\\\ \leq&2\int_{0}^{1}|(F\circ\theta_{i})^{\prime}(y)|\,dy\\\
\leq&2\sum_{j}\int_{j\varepsilon}^{(j+1)\varepsilon}|\theta_{opt}^{2}-\theta_{i}^{2}(y)|\left|\frac{\theta^{ver}_{i,j+1}-\theta^{ver}_{i,j}}{\varepsilon}\right|\textup{
d}y\\\
\leq&2\sum_{j}\int_{j\varepsilon}^{(j+1)\varepsilon}\max\\{|\theta_{opt}^{2}-(\theta^{ver}_{i,j})^{2}|,|\theta_{opt}^{2}-(\theta^{ver}_{i,j+1})^{2}|\\}\left|\frac{\theta^{ver}_{i,j+1}-\theta^{ver}_{i,j}}{\varepsilon}\right|\textup{
d}y\\\
\leq&2\varepsilon\sum_{j}\max\\{|1-\delta-\cos(\theta^{ver}_{i,j+1})|,|1-\delta-\cos(\theta^{ver}_{i,j})|\\}\left|\frac{\theta^{ver}_{i,j+1}-\theta^{ver}_{i,j}}{\varepsilon}\right|\textup{
d}y\ \\\
\leq&8\sum_{j}\max\\{|1-\delta-\cos(\theta^{ver}_{i,j+1})|,|1-\delta-\cos(\theta^{ver}_{i,j})|\\}\left|\sin\left(\frac{\theta^{ver}_{i,j+1}}{2}\right)-\sin\left(\frac{\theta^{ver}_{i,j}}{2}\right)\right|\textup{
d}y\\\ \leq&\frac{8c_{2}\delta^{3/2}}{i\varepsilon}.\end{split}$ (5.49)

Next, by Fubini, there exists $\hat{\theta}\in(0,\theta_{opt}/4)$ such that

$\displaystyle\mathcal{H}^{0}\left(\partial\left\\{\theta_{i}>\frac{\theta_{opt}}{2}-\hat{\theta}\right\\}\right)+\mathcal{H}^{0}\left(\partial\left\\{\theta_{i}<\frac{3\theta_{opt}}{2}+\hat{\theta}\right\\}\right)$
(5.50)
$\displaystyle\quad\leq\mathcal{H}^{0}\left(\partial\left\\{F\circ\theta_{i}>F\left(\frac{\theta_{opt}}{2}-\hat{\theta}\right)\right\\}\right)+\mathcal{H}^{0}\left(\partial\left\\{F\circ\theta_{i}<F\left(\frac{3\theta_{opt}}{2}+\hat{\theta}\right)\right\\}\right)$
(5.51)
$\displaystyle\quad\leq\frac{8}{\theta_{opt}}\int_{0}^{\frac{1}{4}\theta_{opt}}\mathcal{H}^{0}\left(\partial\left\\{F\circ\theta_{i}>F\left(\frac{\theta_{opt}}{2}-s\right)\right\\}\right)\mathcal{H}^{0}\left(\partial\left\\{F\circ\theta_{i}<F\left(\frac{3\theta_{opt}}{2}+s\right)\right\\}\right)\textup{
d}s.$ (5.52)

Consider the functions

$g_{1}(s)\coloneqq F\left(\frac{\theta_{opt}}{2}-s\right)\quad\textup{ and
}\quad g_{2}(s)\coloneqq F\left(\frac{3\theta_{opt}}{2}+s\right).\quad$

Then, $\min\\{|g_{1}^{\prime}(s)|,|g_{2}^{\prime}(s)|\\}\geq
3\theta_{opt}^{2}/4$ for all $s\in(0,\theta_{opt}/4).$ Combining (5.49),
(5.51) and the estimate $\delta^{1/2}\leq\arccos(1-\delta)=\theta_{opt}$ for
$\delta\in(0,1)$, we obtain

$\begin{split}&\mathcal{H}^{0}\left(\partial\left\\{\theta_{i}>\frac{\theta_{opt}}{2}-\hat{\theta}\right\\}\right)+\mathcal{H}^{0}\left(\partial\left\\{\theta_{i}<\frac{3\theta_{opt}}{2}+\hat{\theta}\right\\}\right)\\\
\leq&\frac{8}{\theta_{opt}}\int_{0}^{\frac{1}{4}\theta_{opt}}\mathcal{H}^{0}\left(\partial\left\\{F\circ\theta_{i}>g_{1}(s)\right\\}\right)\textup{
d}s+\frac{8}{\theta_{opt}}\int_{0}^{\frac{1}{4}\theta_{opt}}\mathcal{H}^{0}\left(\partial\left\\{F\circ\theta_{i}<g_{2}(s)\right\\}\right)\textup{
d}t\\\
\leq&\frac{8}{\theta_{opt}}\left(\int_{F\left(\frac{1}{4}\theta_{opt}\right)}^{F\left(\frac{1}{2}\theta_{opt}\right)}\frac{4\mathcal{H}^{0}\left(\partial\left\\{F\circ\theta_{i}>s\right\\}\right)}{3\theta_{opt}^{2}}\textup{
d}s+\int_{F\left(\frac{3}{2}\theta_{opt}\right)}^{F\left(\frac{7}{4}\theta_{opt}\right)}\frac{4\mathcal{H}^{0}(\partial\\{F\circ\theta_{i}<s\\})}{3\theta_{opt}^{2}}\textup{
d}s\right)\\\
\leq&\frac{32}{3\theta_{opt}^{3}}\int_{\mathbb{R}}\mathcal{H}^{0}\left(\partial\left\\{F\circ\theta_{i}>s\right\\}\right)+\mathcal{H}^{0}(\partial\\{F\circ\theta_{i}<s\\}\textup{
d}s\leq\frac{256c_{2}}{i\varepsilon}.\end{split}$ (5.53)

Hence, the set

$M_{+}^{\hat{\theta}}\coloneqq\left\\{y\in(0,1)\colon\frac{\theta_{opt}}{2}-\hat{\theta}\leq\theta_{i}(y)\leq\frac{3\theta_{opt}}{2}+\hat{\theta}\right\\}$
(5.54)

has at most $\mathcal{K}(i\varepsilon)\leq 256c_{2}/(i\varepsilon)$ boundary
points in $(0,1)$. We want to show that the number of discrete boundary points
$\partial^{d,\varepsilon}\mathcal{M}_{+}^{\hat{\theta}}$ of the set

$\mathcal{M}_{+}^{\hat{\theta}}\coloneqq
M_{+}^{\hat{\theta}}\cap\varepsilon\mathbb{Z}=\left\\{j\varepsilon\in(0,1)\cap\varepsilon\mathbb{Z}\colon\frac{\theta_{opt}}{2}-\hat{\theta}\leq\theta^{ver}_{i,j}\leq\frac{3\theta_{opt}}{2}+\hat{\theta}\right\\}$
(5.55)

is bounded by $\mathcal{K}(i\varepsilon)$. If $\mathcal{H}^{0}(\partial
M_{+}^{\hat{\theta}})=0$ we either have
$\theta_{i}\in(\theta_{opt}/2-\hat{\theta},\theta_{opt}/2+\hat{\theta})$ or
$\theta_{i}\notin(\theta_{opt}/2-\hat{\theta},\theta_{opt}/2+\hat{\theta})$
which implies the same property for $\theta^{ver}_{i,\cdot}$ in
$(0,1)\cap\varepsilon\mathbb{Z}.$ Thus, we have
$\mathcal{H}^{0}(\partial^{d,\varepsilon}\mathcal{M}_{+}^{\hat{\theta}})=0$.
If $\mathcal{H}^{0}(\partial M_{+}^{\hat{\theta}})=1$, there is $t\in(0,1)$
such that either $(0,t)=M_{+}^{\hat{\theta}}\setminus\\{t\\}$ or
$(t,1)=M_{+}^{\hat{\theta}}\setminus\\{t\\}$. Hence, there is
$l\varepsilon\in(0,1)\cap\varepsilon\mathbb{Z}$ such that either
$(0,l\varepsilon)\cap\varepsilon\mathbb{Z}=M_{+}^{\hat{\theta}}\setminus\\{l\varepsilon\\}$
or
$(l\varepsilon,1)\cap\varepsilon\mathbb{Z}=M_{+}^{\hat{\theta}}\setminus\\{l\varepsilon\\}$,
which yields $\mathcal{H}^{0}(\partial M_{+}^{\hat{\theta}})=1$. Finally, let
us turn to the case $M\coloneqq\mathcal{H}^{0}(\partial
M_{+}^{\hat{\theta}})\geq 2$. Let $0\eqqcolon
y_{0}<y_{1}<\dots<y_{M}<y_{M+1}\coloneqq 1$ be the boundary points of
$M_{+}^{\hat{\theta}}$ and $l=0,\dots,M-1$. By the definition of
$M_{+}^{\hat{\theta}}$ we either have $(y_{l},y_{l+1})\subseteq
M_{+}^{\hat{\theta}}\setminus\\{y_{l},y_{l+1}\\}$ and
$(y_{l+1},y_{l+2})\subseteq(M_{+}^{\hat{\theta}})^{c}\setminus\\{y_{l+1},y_{l+2}\\}$
or vice versa. Therefore,
$(y_{l},y_{l+1})\cap\varepsilon\mathbb{Z}\subseteq\mathcal{M}_{+}^{\hat{\theta}}\setminus\\{y_{l},y_{l+1}\\}$
and
$(y_{l+1},y_{l+2})\cap\varepsilon\mathbb{Z}\subseteq(\mathcal{M}_{+}^{\hat{\theta}})^{c}\setminus\\{y_{l+1},y_{l+2}\\}$
holds. If $(y_{l},y_{l+1}]\cap\mathcal{M}_{+}^{\hat{\theta}}=\emptyset$ we
have
$\partial^{d,\varepsilon}\mathcal{M}_{+}^{\hat{\theta}}\cap(y_{l},y_{l+2})=\emptyset.$
Otherwise, we have
$\mathcal{K}(i\varepsilon)\coloneqq\mathcal{H}^{0}(\partial^{d,\varepsilon}\mathcal{M}_{+}^{\hat{\theta}}\cap(y_{l},y_{l+2}))=1$,
and thus obtain

$\mathcal{H}^{0}(\partial^{d,\varepsilon}\mathcal{M}_{+}^{\hat{\theta}})\leq\mathcal{H}^{0}(\partial
M_{+}^{\hat{\theta}})\leq\mathcal{K}(i\varepsilon)\leq\frac{256c_{2}}{i\varepsilon}.$
(5.56)

Hence, there exist pairwise disjoint, non-empty discrete intervals of the form
$I_{l}\coloneqq[k_{l}\varepsilon,m_{l}\varepsilon]\cap\varepsilon\mathbb{Z}$,
for $l=0,\dots,\mathcal{K}(i\varepsilon)$ and $k_{l},m_{l}\in\mathbb{N}$,
satisfying

$\bigcup_{l=1}^{\mathcal{K}(i\varepsilon)}I_{l}\subseteq\\{j\varepsilon\in(0,1)\cap\varepsilon\mathbb{Z}:\frac{\theta_{opt}}{2}-\hat{\theta}\leq\theta^{ver}_{i,j}\leq\frac{3\theta_{opt}}{2}+\hat{\theta}\\}\quad\textup{
and
}\quad\varepsilon\mathcal{H}^{0}\left(\bigcup_{l=1}^{\mathcal{K}(i\varepsilon)}I_{l}\right)\geq\frac{1}{12}.$

Further, we introduce a discrete potential
$u\colon\Omega\cap\varepsilon\mathbb{Z}^{2}\rightarrow\mathbb{R}$ of the
angular velocity field $(\theta^{hor},\theta^{ver})$ via

$u(i\varepsilon,j\varepsilon)\coloneqq\frac{\varepsilon}{\theta_{opt}}\sum_{k=0}^{i-1}\theta^{hor}_{k,j}\quad\textup{
for
}(i\varepsilon,j\varepsilon)\in[\varepsilon,1)\cap\varepsilon\mathbb{Z}^{2}$

and $u(0,\cdot)=0.$ Since
$\mathcal{H}^{0}(V)\leq\mathcal{H}^{0}(A^{ver})+\mathcal{H}^{0}(A^{hor})=0$
for the set $V$ of vortices of $(\theta^{hor},\theta^{ver})$, we have

$\partial_{1}^{d,\varepsilon}u(i\varepsilon,j\varepsilon)=\frac{1}{\theta_{opt}}\theta^{hor}_{i,j}\quad\textup{
and
}\quad\partial_{2}^{d,\varepsilon}u(i\varepsilon,j\varepsilon)=\frac{1}{\theta_{opt}}\theta^{ver}_{i,j}.$
(5.57)

Furthermore, we have $\partial_{2}^{d,\varepsilon}u(i\varepsilon,\cdot)\geq
1/4$ for all $j\varepsilon\in I_{l}\cap\varepsilon\mathbb{Z}$ and
$l=1,\dots,\mathcal{K}(i\varepsilon).$ Consider the piecewise affine
interpolation $\tilde{u}\colon[0,1)\rightarrow\mathbb{R}$ given by

$\tilde{u}(y)\coloneqq\frac{y-j\varepsilon}{\varepsilon}u(i\varepsilon,(j+1)\varepsilon)+\frac{(j+1)\varepsilon-y}{\varepsilon}u(i\varepsilon,j\varepsilon)\quad\textup{
for }y\in[j\varepsilon,(j+1)\varepsilon)$

for $j\varepsilon\in[0,1-\varepsilon)\cap\varepsilon\mathbb{Z}.$ The function
$\tilde{u}$ satisfies $\partial_{2}\tilde{u}=\partial_{2}^{d,\varepsilon}u\geq
1/4$ for all $y\in\cup_{l=1}^{\mathcal{K}(i\varepsilon)}I_{l}.$ Using Hölder’s
inequality we find for all $l=1,\dots,\mathcal{K}(i\varepsilon)$

$\frac{1}{12\cdot 16}|I_{l}|^{3}\leq\int_{I_{l}}|\tilde{u}(y)|^{2}\textup{
d}y\leq 2\varepsilon\sum_{j\varepsilon\in
I_{l}}\left(\varepsilon\sum_{k=0}^{i-1}|\partial_{1}^{d,\varepsilon}u(i\varepsilon,j\varepsilon)|\right)^{2}\leq
2\varepsilon^{2}\sum_{j\varepsilon\in
I_{l}}\sum_{k=0}^{i-1}\left|\frac{\theta^{hor}_{i,j}}{\theta_{opt}}\right|^{2}\cdot
i\varepsilon.$ (5.58)

Summing over $l$ leads to

$\begin{split}\frac{(i\varepsilon)^{2}}{(256c_{2})^{2}16\cdot
12^{4}}\leq&\frac{1}{(\mathcal{K}(i\varepsilon))^{2}16\cdot 12\cdot 12^{3}}\\\
=&\frac{\mathcal{K}(i\varepsilon)}{16\cdot
12}\left(\frac{1}{\mathcal{K}(i\varepsilon)}\sum_{l=1}^{\mathcal{K}(i\varepsilon)}|I_{l}|\right)^{3}\\\
\leq&\frac{1}{16\cdot 12}\sum_{l=1}^{\mathcal{K}(i\varepsilon)}|I_{l}|^{3}\\\
\leq&\sum_{l=1}^{K}\int_{I_{l}}|\tilde{u}(y)|^{2}\textup{ d}y\\\
\leq&2\varepsilon\sum_{l=1}^{K}\sum_{j\varepsilon\in
I_{l}}\left(\varepsilon\sum_{k=0}^{i-1}\left|\frac{\theta^{hor}_{k,j}}{\theta_{opt}}\right|^{2}\right)i\varepsilon\\\
\leq&\frac{2\varepsilon^{2}}{\theta_{opt}^{2}}\left(\sum_{j\in\mathcal{J}_{NN}}\sum_{k=0}^{i-1}|\theta_{opt}^{2}-(\theta^{hor}_{k,j})^{2}|\cdot
i\varepsilon\right)+2(i\varepsilon)^{2}.\end{split}$ (5.59)

Using (5.46) we find

$\displaystyle\frac{(i\varepsilon)^{2}}{(256c_{2})^{2}16\cdot 12^{4}}$
$\displaystyle\leq\frac{2\varepsilon^{2}}{\theta_{opt}^{2}}\left(\sum_{j\in\mathcal{J}_{NN}}\sum_{k=0}^{i-1}|\theta_{opt}^{2}-(\theta^{hor}_{k,j})^{2}|\cdot
i\varepsilon\right)+2(i\varepsilon)^{2}$ (5.60)
$\displaystyle\leq\frac{2\varepsilon^{2}}{\theta_{opt}^{2}}\sum_{j\in\mathcal{J}_{NN}}\sum_{k\leq
i}|1-\delta-\cos(\theta^{hor}_{k,j})|\cdot i\varepsilon+2(i\varepsilon)^{2}$
(5.61) $\displaystyle\leq
2\left(\frac{E_{\varepsilon,\delta}(u,\Omega)}{\theta_{opt}^{4}}\right)^{1/2}(i\varepsilon)^{3/2}+2(i\varepsilon)^{2}$
(5.62) $\displaystyle\leq
2\left(\frac{E_{\varepsilon,\delta}(u,\Omega)}{\delta^{2}}\right)^{1/2}(i\varepsilon)^{3/2}+2(i\varepsilon)^{2},$
(5.63)

where we used $\delta^{1/2}\leq\arccos(1-\delta)$ for $\delta\in(0,1)$ in the
last line. Since $(256c_{2})^{2}16\cdot 12^{4}=1/4$ we obtain

$\frac{(i\varepsilon)^{2}}{2(256c_{2})^{2}16\cdot
12^{4}}=\frac{(i\varepsilon)^{2}}{(256c_{2})^{2}16\cdot
12^{4}}-2(i\varepsilon)^{2}\leq
2\left(\frac{E_{\varepsilon,\delta}(u,\Omega)}{\delta^{2}}\right)^{1/2}(i\varepsilon)^{3/2}.$
(5.64)

Rearranging (5.64), using $E_{\varepsilon,\delta}(u,\Omega)\leq
c_{1}s(\varepsilon,\delta)$ and $4(4\cdot 256^{2}\cdot 16\cdot
12^{4})^{2}c_{2}^{3}c_{1}=1$ yields that $i\varepsilon\leq
a(\varepsilon,\delta)$. Hence, (i) or (ii) is satisfied for $i\varepsilon\geq
a(\varepsilon,\delta)$. ∎

######  Remark 5.8.

Careful inspection of the proof shows that the lower bound from Lemma 5.6
holds whenever the angular velocity field $(\theta^{hor},\theta^{ver})$ of a
minimizer $u\in\mathcal{S}_{0}(\Omega,\varepsilon)$ is curl-free. Recalling
that the number of vortices of $u$ is bounded by

$\mathcal{H}^{0}(V)\leq
64\frac{E_{\varepsilon,\delta}(u,\Omega)}{\varepsilon^{2}}=64\frac{\min
E_{\varepsilon,\delta}}{\varepsilon^{2}}.$ (5.65)

Hence, if the right-hand side of (5.65) is smaller than $1$ then
$(\theta^{hor},\theta^{ver})$ is curl-free. Since $\frac{\min
E_{\varepsilon,\delta}}{\varepsilon^{2}}\lesssim\min\left\\{\left(\frac{\delta}{\varepsilon}\right)^{2},\frac{\delta^{3/2}}{\varepsilon}\left(\left|\ln\left(\frac{\varepsilon}{\delta^{1/2}}\right)\right|+1\right)\right\\}$
this proves that a respective lower bound holds for a larger range of
parameters than stated in Lemma 5.6.

######  Remark 5.9 (Existence of vortices).

If a minimizer $u\in\mathcal{S}_{0}(\Omega,\varepsilon)$ is curl-free, we can
use the strategy of Lemma 5.6 to obtain the same lower bound

$\min\left\\{\delta^{2},\varepsilon\delta^{3/2}\left(\left|\ln\left(\frac{\varepsilon}{\delta^{1/2}}\right)\right|+1\right)\right\\}\lesssim
E_{\varepsilon,\delta}(u,\Omega)\lesssim
s(\varepsilon,\delta)\leq\varepsilon\delta^{1/2}.$

Fixing $\delta\in(0,1)$ it holds for $\varepsilon\ll\delta$ that
$\varepsilon\delta^{1/2}\ll\min\\{\delta^{2},\varepsilon\delta^{3/2}\left(|\ln\varepsilon/\delta^{1/2}|+1\right)\\}$.
Hence, there must be minimizers that use vortices to lower their energy.

###  5.3 Scaling law for the periodic problem

In this section, we discuss the periodic problem

$\min\\{E_{\varepsilon,\delta}(u,\Omega)\colon
u\in\mathcal{S}_{0}^{|per|}(\Omega,\varepsilon)\\},$ (5.66)

where $\mathcal{S}_{0}^{|per|}(\Omega,\varepsilon)$ consists of all spin
fields $u\in\mathcal{S}_{0}(\Omega,\varepsilon)$ that additionally satisfy the
periodicity condition (2.24). Note that

$\operatorname{argmin}_{\mathcal{S}_{0}^{|per|}(\Omega,\varepsilon)}(E_{\varepsilon,\delta}(\cdot,\Omega))=\operatorname{argmin}_{\mathcal{S}_{0}^{|per|}(\Omega,\varepsilon)}I_{\alpha,\varepsilon}.$

We prove the following partial scaling law.

######  Theorem 5.10 (Scaling Law of the discrete model with periodic boundary
conditions).

Consider the function $s\colon(0,1/2)\times(0,1)\rightarrow(0,\infty)$ given
by

$s(\varepsilon,\delta)\coloneqq\min\left\\{\delta^{2},\varepsilon\delta^{3/2}\left(\left|\ln\left(\frac{\varepsilon}{\delta^{1/2}}\right)\right|+1\right),\varepsilon\delta^{1/2}\right\\}.$

Then following statements are true for $\delta:=\delta^{ver}=\delta^{hor}$:

  * (i)

There is a constant $C>0$ such that for all
$(\varepsilon,\delta)\in(0,1/2)\times(0,1)$ there holds

$\min\left\\{E_{\varepsilon,\delta}(u;\Omega)\colon
u\in\mathcal{S}_{0}^{|per|}(\Omega,\varepsilon)\right\\}\leq
Cs(\varepsilon,\delta).$

  * (ii)

There is $\overline{\delta}\in(0,1)$ and a constant $c>0$, only depending on
$\overline{\delta},$ such that for all
$(\varepsilon,\delta)\in(0,1/2)\times[\overline{\delta},1)$ it holds

$c\varepsilon\delta^{1/2}\leq\min\left\\{E_{\varepsilon,\delta}(u;\Omega)\colon
u\in\mathcal{S}_{0}^{|per|}(\Omega,\varepsilon)\right\\}$

and for all $\delta\in(0,\overline{\delta})$ and $\varepsilon\in[\delta,1/2)$
it holds

$c\min\left\\{\delta^{2},\varepsilon\delta^{3/2}\left(\left|\ln\left(\frac{\varepsilon}{\delta^{1/2}}\right)\right|+1\right)\right\\}\leq\min\left\\{E_{\varepsilon,\delta}(u;\Omega)\colon
u\in\mathcal{S}_{0}^{|per|}(\Omega,\varepsilon)\right\\}.$

###### Proof.

The lower bound follows immediately from the lower bound on
$\mathcal{S}_{0}(\Omega,\varepsilon)$. For the upper bound we need to adjust
the construction in such a way that the associated spin fields satisfy the
periodicity condition (2.24). The constructions are very similar to the ones
used in Proposition 5.1, and we sketch them only for the readers’ convenience.

Ferromagnet. As in the proof of Proposition 5.1 we choose
$u^{(1)}\coloneqq(0,1)^{T}$ on $\Omega\cap\varepsilon\mathbb{Z}^{2}.$ This
spin field satisfies (2.24) and $E_{\varepsilon,\delta}(u^{(1)})\leq
4\delta^{2}.$

Branching-type construction.  We follow the lines of the proof of the upper
bound of Proposition 5.1. Recall Steps 2 and 3 from the branching-type
construction. Let $N\in\mathbb{N}$ and
$(\Psi_{N}^{hor},\Psi^{ver}_{N})\colon\mathbb{R}^{2}\rightarrow[-\pi,\pi)^{2}$
with
$\varepsilon\textup{curl}^{d,\varepsilon}((\Psi_{N}^{hor},\Psi^{ver}_{N}))=0$
on $\mathbb{R}^{2}\cap\varepsilon\mathbb{Z}^{2}.$ Let $\lambda\in(0,2^{-N}]$.
We slightly adjust the previous construction by setting (see Fig. 15 for a
sketch)

$(\Phi_{N}^{hor},\Phi_{N}^{ver})\coloneqq\begin{cases}(\Psi^{hor}_{N},\Psi^{ver}_{N})\quad&\textup{
on }(-\infty,1-2\lambda]\times\mathbb{R},\\\ (0,-\theta_{opt})\quad&\textup{
on }[1-2\lambda,\infty)\times[1/2,\infty),\\\
(0,\phantom{-}\theta_{opt})\quad&\textup{ on
}[1-2\lambda,\infty)\times(-\infty,1/2]\end{cases}$ (5.67)

Note that $(\Phi_{N}^{hor},\Phi_{N}^{ver})$ again satisfies
$\varepsilon\textup{curl}^{d,\varepsilon}((\Phi_{N}^{hor},\Phi_{N}^{ver}))=0$
on $\mathbb{R}^{2}\cap\varepsilon\mathbb{Z}^{2}.$

Figure 15: Sketch of $(\Phi_{N}^{hor},\Phi_{N}^{ver})$. A corresponding spin
field satisfies (2.24).

We continue by mollifying $(\Phi_{N}^{hor},\Phi_{N}^{ver})$ with a suitable
convolution kernel, i.e., let
$p_{\lambda}\colon\mathbb{R}\rightarrow\mathbb{R}$ be a standard mollifier
satisfying (5.8) and consider
$(\Psi_{N,\lambda}^{hor},\Psi_{N,\lambda}^{ver})\colon\mathbb{R}^{2}\rightarrow[-\pi,\pi)^{2}$,
given for $(x,y)\in\mathbb{R}^{2}$ by (c.f. Figure 12)

$\begin{array}[]{rcl}\Psi_{N,\lambda}^{hor}(x,y)&\coloneqq&\int\limits_{\mathbb{R}^{2}}p_{\lambda}(s)p_{\lambda}(t)\Phi_{N}^{hor}(x-\lambda-s,y-t)\textup{
d}s\textup{d}t\\\
\Psi_{N,\lambda}^{ver}(x,y)&\coloneqq&\int\limits_{\mathbb{R}^{2}}p_{\lambda}(s)p_{\lambda}(t)\Phi_{N}^{ver}(x-\lambda-s,y-t)\textup{
d}s\textup{d}t.\end{array}$ (5.68)

Similarly to the proof of Proposition 5.1, we set

$(\theta_{N,\lambda}^{hor},\theta_{N,\lambda}^{ver})\coloneqq(\Psi_{N,\lambda}^{hor},\Psi_{N,\lambda}^{ver})\quad\textup{
on }\Omega\cap\varepsilon\mathbb{Z}^{2}$

and note that $(\theta_{N,\lambda}^{hor},\theta_{N,\lambda}^{ver})$ induces a
spin field $u^{(2)}\in\mathcal{S}_{0}^{|per|}(\Omega,\varepsilon)$. Estimating
the energy, we proceed analogously as for the branching-type construction in
the proof of Proposition 5.1, and take care of the additional energy
contributions from the extra transition layer that near $\\{1\\}\times(0,1)$,
which yields

$E^{hor}_{\varepsilon,\delta}(u^{(2)}_{N,\lambda},\Omega)\lesssim\delta^{2}\cdot
2^{-N}+N\delta^{2}\lambda+N\frac{\varepsilon^{2}\delta}{\lambda}.$ (5.69)

As before one chooses

$N\coloneqq\left\lceil\left|\frac{\ln\left(\frac{\varepsilon}{\delta^{1/2}}\right)}{\ln(2)}\right|\right\rceil\quad\textup{
and
}\quad\lambda\coloneqq\frac{\varepsilon}{2\delta^{1/2}}\in(0,2^{-N})\subseteq(0,1),$

to obtain the desired bound

$E_{\varepsilon,\delta}(u^{(2)}_{N,\lambda},\Omega)\lesssim\varepsilon\delta^{3/2}\left(\left|\ln\left(\frac{\varepsilon}{\delta^{1/2}}\right)\right|+1\right).$
(5.70)

Vortex structure. Again, we proceed as for the non-periodic case in
Proposition 5.1. Recall the definition of
$(\Psi^{hor},\Psi^{ver})\colon\Omega\cap\varepsilon\mathbb{Z}^{2}\rightarrow[-\pi,\pi)^{2}$
in Steps 1 and 2 there (see Figure 13). We define the angular velocity field
$(\theta_{3}^{hor},\theta^{ver}_{3})\colon\Omega\cap\varepsilon\mathbb{Z}^{2}\rightarrow[-\pi,\pi)^{2}$
first on $[0,1/2]^{2}$ and then extend it to
$\Omega\cap\varepsilon\mathbb{Z}^{2}$ by repeated reflections. We again let
$M\coloneqq\lfloor\pi/\theta_{opt}\rfloor\geq 2$ and set

$(\theta_{3}^{hor},\theta^{ver}_{3})\coloneqq\begin{cases}(\Psi^{hor},\Psi^{ver})\quad&\textup{
on }[0,\min\\{M\varepsilon,1/2\\})\times[0,1/2],\\\
(-\theta_{opt},\theta_{opt})\quad&\textup{ on
}[\min\\{M\varepsilon,1/2\\},1/2]\times[0,1/2],\end{cases}$ (5.71)

We proceed by reflecting the construction (horizontally ) over
$\\{\max([0,1/2]\cap\varepsilon\mathbb{Z})\\}\times[0,1/2)$ and afterwards
(vertically) over $[0,1)\times\\{\max([0,1/2]\cap\varepsilon\mathbb{Z})\\}$,
see Figure 16.

Figure 16: Sketch of the construction $(\theta_{3}^{hor},\theta^{ver}_{3})$
satisfying a periodicity condition.

Let $k_{hor}\coloneqq\lfloor 1/(2\varepsilon)\rfloor$. We extend
$(\theta_{3}^{hor},\theta^{ver}_{3})$ to
$([0,1)\times[0,1/2])\cap\varepsilon\mathbb{Z}^{2}$ by setting

$\theta_{3,k_{hor}+i,j}^{hor}\coloneqq-\theta_{3,k_{hor}-i,j}^{hor}\quad\textup{
and
}\quad\theta^{ver}_{3,k_{hor}+i,j}\coloneqq\theta^{ver}_{3,k_{hor}-i,j}\quad\textup{
for }(i\varepsilon,j\varepsilon)\in[0,1/2)\times[0,1/2].$ (5.72)

Note that this extension still satisfies a discrete curl condition on
$([1/2,1)\times[1/2,1))\cap\varepsilon\mathbb{Z}^{2}$ since

$\varepsilon\textup{curl}^{d,\varepsilon}((\theta_{3}^{hor},\theta^{ver}_{3})_{k_{hor}+i,j})=-\varepsilon\textup{curl}^{d,\varepsilon}((\theta_{3}^{hor},\theta^{ver}_{3})_{k_{hor}-i,j})\in\\{2\pi,0,-2\pi\\}$
(5.73)

for $(i\varepsilon,j\varepsilon)\in[0,1/2)^{2}\cap\varepsilon\mathbb{Z}^{2}$.
Analogously, we extend $(\theta_{3}^{hor},\theta^{ver}_{3})$ vertically to the
whole domain $\Omega\cap\varepsilon\mathbb{Z}^{2}$. Since the discrete curl-
condition is satisfied, there is a spin field
$u^{(3)}\colon\Omega\cap\varepsilon\mathbb{Z}^{2}\rightarrow\mathbb{S}^{1}$
whose associated angular velocity field is given by
$(\theta_{3}^{hor},\theta_{3}^{ver})$. The energy can be bounded similarly to
the non-periodic case, including the second vertical layer of vortices on the
left boundary and the additional horizontal and vertical transition layers,
see Figure 16. Hence, we obtain the bound

$E_{\varepsilon,\delta}(u^{(3)},\Omega)\lesssim\varepsilon\delta^{1/2}+\varepsilon\delta\lesssim\varepsilon\delta^{1/2}.$
(5.74)

This concludes the proof of the upper bound of Theorem 5.10. ∎

## Acknowledgement

The authors gratefully acknowledge funding of the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) via project 195170736
- TRR 109. In addition, the authors would like to thank Marco Cicalese and
Marwin Forster for helpful discussions. Furthermore, Melanie Koser would also
like to thank Franz Bethke for illuminating numerical experiments achieved in
a collaboration.

## Appendix

In the following we recall a heuristic argument based on [CS15, CFO19] for the
Gamma-convergence result in Section 4. We adjust the presentation to the
anisotropic model for the interaction parameters
$\delta^{hor},\delta^{ver}\in(0,1).$ For $\varepsilon>0$ small enough one
formally has
$(u_{i+2,j}-2u_{i+1,j}+u_{i,j})/\varepsilon^{2}\approx\partial_{11}u$ and
similarly
$(u_{i,j+2}-2u_{i,j+1}+u_{i,j})/\varepsilon^{2}\approx\partial_{22}u$ for some
smooth enough function $u\colon\Omega\rightarrow S^{1}$ (see below). This
consideration and a formal computation [CS15, CFO19] lead to

$\begin{split}&E_{\varepsilon,\Delta}(u,\Omega)=\frac{\varepsilon^{2}}{2}\sum_{(i,j)}|u_{i,j}-2(1-\delta^{hor})u_{i+1,j}+u_{i+2,j}|^{2}+\frac{\varepsilon^{2}}{2}\sum_{(i,j)}|u_{i,j}-2(1-\delta^{ver})u_{i,j+1}+u_{i,j+2}|^{2}\\\
&\quad=\frac{\varepsilon^{2}}{2}\sum_{(i,j)}\left|\varepsilon^{2}\frac{u_{i+2,j}-2u_{i+1,j}+u_{i,j}}{\varepsilon^{2}}+2\delta^{hor}u_{i+1,j}\right|+\frac{\varepsilon^{2}}{2}\sum_{(i,j)}\left|\varepsilon^{2}\frac{u_{i,j+2}-2u_{i,j+2}+u_{i,j}}{\varepsilon^{2}}+2\delta^{ver}u_{i,j+2}\right|\\\
&\quad\simeq\frac{1}{2}\int|\varepsilon^{2}\partial_{11}u+2\delta^{hor}u|^{2}+\left|\varepsilon^{2}\partial_{22}u+2\delta^{ver}u\right|\operatorname{\textup{
d}\textit{x}}\textup{d}y\\\ &\quad\simeq\frac{1}{2}\int
4(\delta^{hor})^{2}+4\delta^{hor}\varepsilon^{2}u\partial_{11}u+\varepsilon^{4}|\partial_{11}u|^{2}+4(\delta^{ver})^{2}+4\delta^{ver}\varepsilon^{2}u\partial_{22}u+\varepsilon^{4}|\partial_{22}u|^{2}\operatorname{\textup{
d}\textit{x}}\textup{d}y\\\ &\quad\simeq\frac{1}{2}\int
4(\delta^{hor})^{2}-4\delta^{hor}\varepsilon^{2}|\partial_{1}u|^{2}+\varepsilon^{4}|\partial_{11}u|^{2}+4(\delta^{ver})^{2}-4\delta^{ver}\varepsilon^{2}|\partial_{2}u|^{2}+\varepsilon^{4}|\partial_{22}u|^{2}\operatorname{\textup{
d}\textit{x}}\textup{d}y\end{split}$ (5.75)

Using polar coordinates, we write $u$ in the form $u=(\cos(\Psi),\sin(\Psi))$
for some suitable $\Psi\colon\Omega\to[-\pi,\pi)$. Assuming that $\Psi$ if
differentiable and using $|\partial_{i}u|^{2}=|\partial_{i}\Psi|^{2}$ and
$|\partial_{ii}u|^{2}=|\partial_{ii}\Psi|^{2}+|\partial_{i}\Psi|^{4}$ for
$i=1,2$, we formally obtain

$\begin{split}&E_{\varepsilon,\Delta}(u,\Omega)\\\ \simeq&\frac{1}{2}\int
4(\delta^{hor})^{2}-4(\delta^{hor})\varepsilon^{2}|\partial_{1}u|^{2}+\varepsilon^{4}|\partial_{11}u|^{2}+4(\delta^{ver})^{2}-4\delta^{ver}\varepsilon^{2}|\partial_{2}u|^{2}+\varepsilon^{4}|\partial_{22}u|^{2}\operatorname{\textup{
d}\textit{x}}\textup{d}y\\\
&=\frac{1}{2}\int(2\delta^{hor}-\varepsilon^{2}|\partial_{1}\Psi|^{2})^{2}+\varepsilon^{4}|\partial_{11}\Psi|^{2}+(2\delta^{ver}-\varepsilon^{2}|\partial_{2}\Psi|^{2})^{2}+\varepsilon^{4}|\partial_{22}\Psi|^{2}\operatorname{\textup{
d}\textit{x}}\textup{d}y\\\ &=\frac{1}{2}\int
4(\delta^{ver})^{2}\left(\frac{\delta^{hor}}{\delta^{ver}}-\left|\frac{\varepsilon}{\sqrt{2\delta^{ver}}}\partial_{1}\Psi\right|^{2}\right)^{2}+\varepsilon^{4}|\partial_{11}\Psi|^{2}\operatorname{\textup{
d}\textit{x}}\textup{d}y\\\ &\phantom{=}+\frac{1}{2}\int
4(\delta^{ver})^{2}\left(1-\left|\frac{\varepsilon}{\sqrt{2\delta^{ver}}}\partial_{2}\Psi\right|^{2}\right)^{2}+\varepsilon^{4}|\partial_{22}\Psi|^{2}\operatorname{\textup{
d}\textit{x}}\textup{d}y\end{split}$ (5.76)

Substituting $\varphi=\varepsilon\Psi/\sqrt{2\delta^{ver}}$ and
$\gamma^{2}\coloneqq\delta^{hor}/\delta^{ver}$ finally leads to

$\begin{split}&E_{\varepsilon,\Delta}(u,\Omega)\\\ \simeq&\frac{1}{2}\int
4(\delta^{ver})^{2}\left(\gamma^{2}-\left|\partial_{1}\varphi\right|^{2}\right)^{2}+2\delta^{ver}\varepsilon^{2}|\partial_{11}\varphi|^{2}+4(\delta^{ver})^{2}\left(1-\left|\partial_{2}\varphi\right|^{2}\right)^{2}+2\delta^{ver}\varepsilon^{2}|\partial_{22}\varphi|^{2}\operatorname{\textup{
d}\textit{x}}\textup{d}y\\\
=&\sqrt{2}\varepsilon(\delta^{ver})^{3/2}\int\bigg{[}\frac{\sqrt{2\delta^{ver}}}{\varepsilon}\left(\gamma^{2}-\left|\partial_{1}\varphi\right|^{2}\right)^{2}+\frac{\varepsilon}{\sqrt{2\delta^{ver}}}|\partial_{11}\varphi|^{2}+\frac{\sqrt{2\delta^{ver}}}{\varepsilon}\left(1-\left|\partial_{2}\varphi\right|^{2}\right)^{2}\\\
&\hskip
85.35826pt+\frac{\varepsilon}{\sqrt{2\delta^{ver}}}|\partial_{22}\varphi|^{2}\bigg{]}\operatorname{\textup{
d}\textit{x}}\textup{d}y\end{split}$ (5.77)

Assuming small angular changes of the spin field $u$ (see Section 4 for a
rigorous proof) in the horizontal and vertical direction formally yields

$\partial_{1}\varphi=\frac{\varepsilon}{\sqrt{2\delta^{ver}}}\partial_{1}\Psi\simeq\frac{1}{\sqrt{2\delta^{ver}}}\theta^{hor}=\gamma\frac{1}{\sqrt{2\delta^{hor}}}\theta^{hor}\simeq\gamma\sqrt{\frac{2}{\delta^{hor}}}\sin\left(\frac{\theta^{hor}}{2}\right)\eqqcolon\gamma
w^{\Delta}(u)$

and similarly

$\partial_{2}\varphi=\frac{\varepsilon}{\sqrt{2\delta^{ver}}}\partial_{2}\Psi\simeq\sqrt{\frac{2}{\delta^{ver}}}\sin\left(\frac{\theta^{ver}}{2}\right)\eqqcolon
z^{\Delta}(u),$

where $w^{\Delta}(u)$ and $z^{\Delta}(u)$ are the chirality order parameters
introduced in Section 2. Setting
$\sigma\coloneqq\varepsilon/\sqrt{2\delta^{ver}}$ and the formal equality
$\nabla\varphi=(\gamma w^{\Delta}(u),z^{\Delta}(u))$, we obtain the
renormalized energy

$\begin{split}&H_{\sigma,\gamma}(u,\Omega)\\\
\coloneqq&\frac{1}{\sigma}\int(\gamma^{2}-|\partial_{1}\varphi|^{2})^{2}\operatorname{\textup{
d}\textit{x}}\textup{d}y+\sigma\int|\partial_{11}\varphi|^{2}\operatorname{\textup{
d}\textit{x}}\textup{d}y+\frac{1}{\sigma}\int(1-|\partial_{2}\varphi|^{2})^{2}\operatorname{\textup{
d}\textit{x}}\textup{d}y+\sigma\int|\partial_{22}\varphi|^{2}\operatorname{\textup{
d}\textit{x}}\textup{d}y\\\
\simeq&\frac{1}{\sqrt{2}\varepsilon(\delta^{ver})^{3/2}}E_{\varepsilon,\Delta}(u,\Omega).\end{split}$
(5.78)

## References

  * [ABC+22] R. Alicandro, A. Braides, M. Cicalese, L. De Luca, and A. Piatnitski.  Topological singularities in periodic media: Ginzburg–Landau and core-radius approaches.  Arch. Rat. Mech. Anal., 243(2):559–609, 2022.
  * [ABCS23] R. Alicandro, A. Braides, M. Cicalese, and M. Solci.  Discrete Variational Problems with Interfaces.  Cambridge Monographs on Applied and Computational Mathematics. Cambridge University Press, 2023.
  * [AC09] R. Alicandro and M. Cicalese.  Variational analysis of the asymptotics of the ${{XY}}$ model.  Arch. Rat. Mech. Anal., 192(3):501–36, 2009.
  * [ACP11] R. Alicandro, M. Cicalese, and M. Ponsiglione.  Variational equivalence between Ginzburg-Landau, ${{XY}}$ spin systems and screw dislocation energies.  Indiana Univ. Math. J., 60(1):171–208, 2011.
  * [ACR15] R. Alicandro, M. Cicalese, and M. Ruf.  Domain formation in magnetic polymer composites: an approach via stochastic homogenization.  Arch. Rat. Mech. Anal., 218(2):945–984, 2015.
  * [ADGP16] R. Alicandro, L. De Luca, A. Garroni, and M. Ponsiglione.  Dynamics of discrete screw dislocations on glide directions.  J. Mech. Phys. Solids, 92:87–104, 2016.
  * [ADLGP14] R. Alicandro, L. De Luca, A. Garroni, and M. Ponsiglione.  Metastability and dynamics of discrete topological singularities in two dimensions: A $\Gamma$-convergence approach.  Arch. Rat. Mech. Anal., 214(1):269–330, 2014.
  * [AGZ24] L. Abel, J. Ginster, and B. Zwicknagl.  A scaling law for a model of epitaxially strained elastic films with dislocations, 2024.  arXiv:2403.13646.
  * [BCDLP18] R. Badal, M. Cicalese, L. De Luca, and M. Ponsiglione.  $\Gamma$-convergence analysis of a generalized xy model: Fractional vortices and string defects.  Comm. Math. Phys., 358(2):705–739, 2018.
  * [BCDM02] H. Ben Belgacem, S. Conti, A. DeSimone, and S. Müller.  Energy scaling of compressed elastic films – three-dimensional elasticity and reduced theories.  Arch. Rat. Mech. Anal., 164:1–37, 2002.
  * [BCL05] H. Brézis, P. G. Ciarlet, and J. L. Lions.  Analyse fonctionnelle: Théorie et applications.  Sciences sup. Dunod, 2005.
  * [BCM17] D. Bourne, S. Conti, and S. Müller.  Energy bounds for a compressed elastic film on a substrate.  J. Nonlinear Science, 27:453–494, 2017.
  * [BG15] P. Bella and M. Goldman.  Nucleation barriers at corners for cubic-to-tetragonal phase transformation.  Proc. Roy. Soc. Edinburgh A, 145:715–724, 2015.
  * [BK14] P. Bella and R. V. Kohn.  Wrinkles as the result of compressive stresses in an annular thin film.  Comm. Pure Appl. Math., 67(5):693–747, 2014.
  * [BK24] F. Bethke and M. Koser.  Towards the numerical optimization of frustrated spin systems at the helimagnetic/ferromagnetic transition point.  In preperation, 2024.
  * [BT08] A. Braides and L. Truskinovsky.  Asymptotic expansions by $\Gamma$-convergence.  Continuum Mech. Thermodyn., pages 21–62, 2008.
  * [CC15] A. Chan and S. Conti.  Energy scaling and branched microstructures in a model for shape-memory alloys with SO(2) invariance.  Mathematical Models Methods App. Sci., 25:1091–1124, 2015.
  * [CDKZ21] S. Conti, J. Diermeier, M. Koser, and B. Zwicknagl.  Asymptotic self-similarity of minimizers and local bounds in a model of shape-memory alloys.  J. Elast., 147:149–200, 2021.
  * [CDMZ20] S. Conti, J. Diermeier, D. Melching, and B. Zwicknagl.  Energy scaling laws for geometrically linear elasticity models for microstructures in shape memory alloys.  ESAIM Control Optim. Calc. Var., 26:115, 2020.
  * [CFO19] M. Cicalese, M. Forster, and G. Orlando.  Variational analysis of a two-dimensional frustrated spin system: emergence and rigidity of chirality transitions.  SIAM J. Math. Anal., 51(6):4848–4893, 2019.
  * [CFO22] M. Cicalese, M. Forster, and G. Orlando.  Variational analysis of the $J_{1}-J_{2}-J_{3}$ model: A non-linear lattice version of the aviles–giga functional.  Arch. Rat. Mech. Anal., 245(2):1059–1133, 2022.
  * [CGOS18] S. Conti, M. Goldman, F. Otto, and S. Serfaty.  A branched transport limit of the Ginzburg-Landau functional.  J. Éc. polytech. Math., 5:317–375, 2018.
  * [CK98] R. Choksi and R. V. Kohn.  Bounds on the micromagnetic energy of a uniaxial ferromagnet.  Comm. Pure Appl. Math., 51(3):259–289, 1998.
  * [CKO99] R. Choksi, R. V. Kohn, and F. Otto.  Domain branching in uniaxial ferromagnets: a scaling law for the minimum energy.  Comm. Math. Phys., 201(1):61–79, 1999.
  * [CKO04] R. Choksi, R. V. Kohn, and F. Otto.  Energy minimization and flux domain structure in the intermediate state of a type-I superconductor.  J. Nonlinear Sci., 14:119–171, 2004.
  * [Con00] S. Conti.  Branched microstructures: Scaling and asymptotic self-similarity.  Comm. Pure Appl. Math., 53(11):1448–1474, 2000.
  * [Con06] S. Conti.  A lower bound for a variational model for pattern formation in shape-memory alloys.  Contin. Mech. Thermodyn., 17:469–476, 2006.
  * [COR23] M. Cicalese, G. Orlando, and M. Ruf.  A classical $S^{2}$ spin system with discrete out-of-plane anisotropy: Variational analysis at surface and vortex scalings.  Nonlinear Analysis, 231:112929, 2023.
  * [COS16] S. Conti, F. Otto, and S. Serfaty.  Branched microstructures in the Ginzburg-Landau model of type-I superconductors.  SIAM J. Math. Anal., 48:2994–3034, 2016.
  * [CS15] M. Cicalese and F. Solombrino.  Frustrated ferromagnetic spin chains: a variational approach to chirality transitions.  J. Nonlinear Sci., 25(2):291–313, 2015.
  * [Die13] H. T. Diep, editor.  Frustrated spin systems.  World Scientific, 2nd edition, 2013.
  * [DKMO02] A. DeSimone, R. V. Kohn, S. Müller, and F. Otto.  A reduced theory for thin-film micromagnetics.  Commun. Pure Appl. Math., 55:1408–1460, 2002.
  * [DVJ19] V. Dabade, R. Venkatraman, and R.D. James.  Micromagnetics of Galfenol.  J. Nonlinear Sci., 29:415–460, 2019.
  * [EG92] L. C. Evans and R. F. Gariepy.  Measure theory and fine properties of functions.  Studies in advanced mathematics. CRC Press, 1992.
  * [For18] M. Forster.  Discrete-to-Continuum Gamma-Limit of a Planar Chirality Model.  Master thesis, 2018.
  * [Gin23] J. Ginster.  On the formation of microstructure for singularly perturbed problems with 2,3 or 4 preferred gradients, 2023.  CVGMT preprint: https://cvgmt.sns.it/paper/6022/.
  * [Gol20] M. Goldman.  Self-similar minimizers of a branched transport functional.  Indiana Univ. Math. J., 69(4):1073–1104, 2020.
  * [GRTZ24] J. Ginster, A. Rüland, A. Tribuzio, and B. Zwicknagl.  On the effect of geometry on scaling laws for a class of martensitic phase transformations, 2024.  arXiv:2405.05927.
  * [GZ23] J. Ginster and B. Zwicknagl.  Energy scaling law for a singularly perturbed four-gradient problem in helimagnetism.  J. Nonlinear Sci., 33(1):20, 2023.
  * [GZ24] J. Ginster and B. Zwicknagl.  Energy scaling laws for microstructures: from helimagnets to martensites.  Calc. Var. Partial Differential Equations, 63(1):8, 2024.
  * [JS01] W. Jin and P. Sternberg.  Energy estimates of the von Kármán model of thin-film blistering.  J. Math. Phys., 42:192–199, 2001.
  * [KKO13] H. Knüpfer, R. V. Kohn, and F. Otto.  Nucleation barriers for the cubic-to-tetragonal phase transformation.  Comm. Pure Appl. Math., 66(6):867–904, 2013.
  * [KL23] A. Kubin and L. Lamberti.  Variational analysis in one and two dimensions of a frustrated spin system: chirality and magnetic anisotropy transitions.  Math. Eng., 5(6):1–37, 2023.
  * [KM92] R.V. Kohn and S. Müller.  Branching of twins near an austenite-twinned martensite interface.  Phil. Mag. A, 66:697–715, 1992.
  * [KM94] R. V. Kohn and S. Müller.  Surface energy and microstructure in coherent phase transitions.  Comm. Pure Appl. Math., 47(4):405–435, 1994.
  * [KM11] H. Knüpfer and C.B. Muratov.  Domain structure of bulk ferromagnetic crystals in applied fields near saturation.  J. Nonlin. Sc., 21:921–962, 2011.
  * [KO19] H. Knüpfer and F. Otto.  Nucleation barriers for the cubic-to-tetragonal phase transformation in the absence of self-accommodation.  Z. Angew. Math. Mech., 99:e201800179, 2019.
  * [Koh06] R.V. Kohn.  Energy-driven pattern formation.  In International Congress of Mathematicians, ICM 2006, volume 1, pages 359–383. EMS Publishing House, 2006.
  * [Kos24] M. Koser.  On singularly-perturbed variational problems for pattern formation in helimagnets and martensites.  PhD thesis, Humboldt-Universität zu Berlin, submitted 2024.
  * [MBJ+09] S. Mühlbauer, B. Binz, F. Jonietz, C. Pfleiderer, A. Rosch, A. Neubauer, R. Georgii, and P. Böni.  Skyrmion lattice in a chiral magnet.  Science, 323(5916):915–919, 2009.
  * [MNA+10] W. Münzer, A. Neubauer, T. Adams, S. Mühlbauer, C. Franz, F. Jonietz, R. Georgii, P. Böni, B. Pedersen, M. Schmidt, A. Rosch, and C. Pfleiderer.  Skyrmion lattice in the doped semiconductor ${\text{fe}}_{1-x}{\text{co}}_{x}\text{Si}$.  Phys. Rev. B, 81:041203, 2010.
  * [OV10] F. Otto and T. Viehmann.  Domain branching in uniaxial ferromagnets: asymptotic behavior of the energy.  Calc. Var. Partial Differential Equations, 38(1-2):135–181, 2010\.
  * [RR23] T. Ried and C. Román.  Domain branching in micromagnetism: scaling law for the global and local energies, 2023.  arXiv:2308.07239.
  * [RT22] A. Rüland and A. Tribuzio.  On the energy scaling behaviour of a singularly perturbed tartar square.  Arch. Rational Mech. Anal., 243:401–431, 2022.
  * [Sim21a] T. M. Simon.  Quantitative aspects of the rigidity of branching microstructures in shape memory alloys via h-measures.  SIAM J. Math. Anal., 53(4):4537–4567, 2021.
  * [Sim21b] T. M. Simon.  Rigidity of branching microstructures in shape memory alloys.  Arch. Ration. Mech. Anal., 241:1707–1783, 2021.
  * [SPD+20] H. Seiner, P. Plucinsky, V. Dabade, B. Benešová, and R.D. James.  Branching of twins in shape memory alloys revisited.  J. Mech. Phys. Solids, 141:103961, 2020.
  * [SV18] G. Scilla and V. Valocchia.  Chirality transitions in frustrated ferromagnetic spin chains: A link with the gradient theory of phase transitions.  J. Elasticity, 132:271–293, 2018.
  * [Tri78] H. Triebel.  Interpolation Theory, Function Spaces, Differential Operators, volume 18 of North-Holland Mathematical Library.  Elsevier, 1978.
  * [UOMT06] M. Uchida, Y. Onose, Y. Matsui, and Y. Tokura.  Real-space observation of helical spin order.  Science, 311(5759):359–361, 2006.
  * [VDJ20] R. Venkatraman, V. Dabade, and R. D. James.  Bounds on the energy of a soft cubic ferromagnet with large magnetostriction.  J. Nonlinear Sci., 30(6):3367–3388, 2020.
  * [YKO+11] X. Z. Yu, N. Kanazawa, Y. Onose, K. Kimoto, W. Z. Zhang, S. Ishiwata, Y. Matsui, and Y. Tokura.  Near room-temperature formation of a skyrmion crystal in thin-films of the helimagnet fege.  Nat. Mater., 10(2):106–109, 2011.
  * [Zie89] W. P. Ziemer.  Weakly Differentiable Functions.  Graduate Texts in Mathematics, 120. Springer, New York, NY, 1989.
  * [Zwi14] B. Zwicknagl.  Microstructures in low-hysteresis shape memory alloys: Scaling regimes and optimal needle shapes.  Arch. Rat. Mech. Anal., 213(2):355–421, 2014.


# High Performance of Gradient Boosting in Binding Affinity Prediction

Dmitrii Gavrilev  Nurlybek Amangeldiuly  Sergei Ivanov  Evgeny Burnaev

###### Abstract

Prediction of protein-ligand (PL) binding affinity remains the key to drug
discovery. Popular approaches in recent years involve graph neural networks
(GNNs), which are used to learn the topology and geometry of PL complexes.
However, GNNs are computationally heavy and have poor scalability to graph
sizes. On the other hand, traditional machine learning (ML) approaches, such
as gradient-boosted decision trees (GBDTs), are lightweight yet extremely
efficient for tabular data. We propose to use PL interaction features along
with PL graph-level features in GBDT. We show that this combination
outperforms the existing solutions.

Machine Learning, ICML

##  1 Introduction

Binding affinity is one of the major determinants to identify whether a
molecule binds to a desired macromolecular target. It helps speculate if the
molecule can induce a biological response in an organism, which is the central
mechanism in drug action. Reliable and fast prediction of such parameter is an
important and challenging task in cheminformatics (Wan et al., 2020). To build
an accurate binding affinity prediction model PL complex 3D structures are
needed (Son & Kim, 2021), and are readily accessible from steadily growing
open source databases (Jumper et al., 2021). While docking is an established
approach to predict affinity (Kitchen et al., 2004), recent ML-based
algorithms consistently outperform existing classical methods (Meli et al.,
2021). Considered a natural way of processing chemical structures, recently
GNNs are proposed to predict binding affinity. In a typical prediction
pipeline for cheminformatics, a GNN generalizes to a labeled and weighted
graph, nodes of which represent molecular atoms, and edges may represent their
bonds by including various interaction features such as bond angles and
distances (Jones et al., 2021; Son & Kim, 2021). This application domain of
GNN can lead to a promising results (Kearnes et al., 2016). The main
bottleneck of GNNs is, however, their computationally complex iterative
convergence (Gupta et al., 2020; Tiezzi et al., 2020).

In contrast, GBDTs provide faster convergence and perform out-of-the-box,
which makes them a popular tool for ML contests and industrial applications
(Ivanov & Prokhorenkova, 2021; Gorishniy et al., 2021). Notwithstanding the
advantages, powerful model of GBDTs often overlooked (Qin et al., 2021).

In this work, we leverage the advantages of GBDTs to graph-level structural
and interaction feature representations of PL complexes for binding affinity
prediction. Assessment of the model performance on external test sets showed
that the GBDTs can perform equally well or outperform the recent state-of-the-
art GNN-based approach, SIGN (Li et al., 2021).

Table 1: Metrics on PDBbind core and CSAR-HiQ sets. Models | PDBbind core set | CSAR-HiQ set
---|---|---
RMSE $\downarrow$ |  MAE $\downarrow$ |  SD $\downarrow$ |  R $\uparrow$ |  RMSE $\downarrow$ |  MAE $\downarrow$ |  SD $\downarrow$ |  R $\uparrow$
LR | 1.608 (0.000) | 1.285 (0.000) | 1.576 (0.000) | 0.690 (0.000) | 1.994 (0.000) | 1.555 (0.000) | 1.909 (0.000) | 0.650 (0.000)
MLP | 1.456 (0.032) | 1.158 (0.025) | 1.454 (0.033) | 0.744 (0.014) | 2.195 (0.203) | 1.727 (0.178) | 2.046 (0.116) | 0.575 (0.066)
SIGN | 1.316 (0.031) | 1.027 (0.025) | 1.312 (0.035) | 0.797 (0.012) | 1.735 (0.031) | 1.327 (0.040) | 1.709 (0.044) | 0.754 (0.014)
Ours: CatBoost | 1.321 (0.008) | 1.045 (0.011) | 1.270 (0.012) | .812 (0.004) | 1.798 (0.031) | 1.391 (0.040) | 1.679 (0.014) | 0.744 (0.005)
Ours: LightGBM | 1.316 (0.010) | 1.040 (0.007) | 1.279 (0.015) | 0.809 (0.005) | 1.725 (0.038) | 1.305 (0.046) | 1.660 (0.049) | 0.751 (0.017)

##  2 Related Work

The early ML methods of binding affinity prediction rely on careful feature
engineering. Namely, RF-Score proposed training a Random Forest model on
intermolecular interaction features (Ballester & Mitchell, 2010). For sets of
atom types $P$ in protein and $L$ in ligand, the authors count the number of
atom type pairs in the complex. More formally, the resulting interaction
features vector $x$ is $|P|\times|L|$ dimensional and

$x_{ji}=\sum_{a\in protein}\sum_{b\in
ligand}\delta_{aj}\delta_{bi}\Theta\left(d_{cutoff}-d_{ab}\right),$

where $\delta$ and $\Theta$ are the Kronecker and Heavyside functions
respectively, $d_{cutoff}$ is the interaction cutoff distance and $d_{ab}$ is
the distance between atoms $a$ and $b$. Particularly, the authors set
$d_{cutoff}=12$. As for $P$ and $L$, the following sets are recommended:

$\displaystyle P=\left\\{C,N,O,S\right\\},$ $\displaystyle
L=\left\\{C,N,O,F,P,S,Cl,Br,I\right\\}.$

Therefore, the interaction features can be described by a 36-dimensional
vector. Designing representations for the molecules has also received a great
deal of attention in the context of graph kernels (Kriege et al., 2020; Ivanov
& Burnaev, 2018) which have been applied successfully for various applications
on graphs (Sharaev et al., 2018; Ivanov et al., 2018). On the other hand,
recent trends in binding affinity prediction involve extensive usage of GNNs
to obtain spatial representations. Such approach also allows characterizing
intermolecular and intramolecular interactions. To enrich the representation,
the efforts are being directed towards integrating atomic distances (Lim et
al., 2019), bond angles and physics-derived descriptors (Moon et al., 2022).
In such work (Jiang et al., 2021), authors construct protein graphs, ligand
graphs and PL graphs separately to represent PL structure and interaction by
also integrating distance, angle, and various 2D node and edge features.
Another approach, graphDelta (Karlov et al., 2020) characterizes binding site
of the complex by computing Behler-Parrinello symmetry functions (Behler &
Parrinello, 2007) per each common PL heavy atom types in the complex and
obtains complex interaction features to train a message-passing neural
network. Very recent work, SIGN exploits the interaction features similar to
RF-Score. Moreover, the authors of SIGN propose to use not only topological
information but also angles and distances. Despite outperforming ML and NN-
based models, following drawbacks inherent to GNN models are not solved in
SIGN:

  * •

Time-consuming construction of graphs from raw 3D coordinates (Wang & Xia,
2012);

  * •

Requires large computational resources for training (Liu et al., 2020).

We provide time and memory consumption details of SIGN in Section 4.

##  3 Our Approach

In this work, we incorporate PL interaction features, as introduced in RF-
Score, alongside graph-level features. To construct graph-level features, the
atomic features are pooled together. For each atomic feature, we sum the
values over all atoms in the complex. Additionally, we calculate the standard
deviation of atomic features to characterize atomic distribution statistics.
Categorical features were encoded to a one-hot numeric array: atom types were
transformed into a binary array, where each bit encodes a type. Then, the
graph-level feature constructed from these arrays represents the distribution
of atom types in a complex. Below we list the features that we used to train
our models:

  * •

RF-Score PL interaction features;

  * •

Summarized graph-level atomic features;

  * •

Standard deviation of atomic features;

After obtaining features, we pass them to a GBDT model. GBDT (Friedman, 2001)
is a widely used algorithm that builds ensemble of decision trees. The idea of
GBDT is to iteratively add weak decision trees and aggregate their predictions
to obtain an overall strong model. At each iteration the model is updated by
building a new tree that minimizes the loss function w.r.t. previous
predictions of the ensemble. In particular, at each iteration $t$ the model
$f(\mathbf{x})$ is updated in an additive manner:

$f^{t}(\mathbf{x})=f^{t-1}(\mathbf{x})+\epsilon\,h^{t}(\mathbf{x}),$ (1)

where $f^{t-1}$ is previously built trees, $h^{t}$ is a weak learner that is
chosen from some family of functions $\mathcal{H}$ (e.g. decision trees), and
$\epsilon$ is a learning rate. The weak learner $h^{t}\in\mathcal{H}$
approximates the direction of the negative gradient of a loss function $L$
w.r.t. the current model’s predictions:

$h^{t}=\operatorname*{arg\,min}_{h\in\mathcal{H}}\sum_{i}\left(-\frac{\partial
L(f^{t-1}(\mathbf{x}_{i}),y_{i})}{\partial
f^{t-1}(\mathbf{x}_{i})}-h(\mathbf{x}_{i})\right)^{2}.$ (2)

As such one can think of GBDT’s training as functional gradient descent, where
at each iteration we move in the direction that minimizes the loss. This
allows flexibly training GBDT on various tabular continuous/discrete, noisy
and heterogeneous data.

##  4 Experiments

### Data.

For evaluating the performance of the models, we use the same dataset and
splits as in (Li et al., 2021). Specifically, all experiments were conducted
on PDBbind v2016 (Liu et al., 2015) and CSAR-HiQ (Dunbar Jr et al., 2011),
which contain experimental PL binding affinity data. The core set of PDBbind
(290 complexes) was used as a benchmark for model comparison (Su et al.,
2018). The models were trained on the difference between the refined and core
sets (3767 complexes). For tuning hyperparameters, 377 out of 3767 complexes
were used as a validation set. After training and tuning, the models were
additionally evaluated on CSAR-HiQ dataset. Since CSAR-HiQ and the refined set
of PDBbind overlap, we filter out these overlapping complexes for evaluation.
The filtered dataset consists of 135 complexes.

Figure 1: Feature importance (interaction features).

We follow the same feature extraction procedure as the authors of SIGN. Each
atom is represented by a 36-dimensional vector, 18 of which describe protein
features and the rest describe ligand features. Those features include atom
type, hybridization, partial charge, aromaticity, etc. For this set of
features, we run the pipeline presented in Section 3. As a result, we obtain
108 features for each complex — 36 from the interaction matrix, 36 by summing
atomic features, and the other 36 by calculating standard deviation.

### Models.

Our work focuses mainly on GBDT models since they have shown great performance
for tabular data in both time and quality. We compare two popular GBDT models,
CatBoost (Dorogush et al., 2018) and LightGBM (Ke et al., 2017), with SIGN.
First, we train the models on PDBbind refined set and tune their
hyperparameters on validation set. For tuning, we utilize Optuna (Akiba et
al., 2019), automatic hyperparameter search framework. While validating, we
aim to minimize root mean squared error (RMSE). The best hyperparameters for
CatBoost are following: depth = 10, grow policy = symmetric tree, L2 leaf
regularization = 73.47. As for LightGBM, the optimal hyperparameters are:
bagging fraction = 1, bagging frequency = 0, feature fraction = 0.8, L1
regularization = 0, L2 regularization = 0, minimal child samples = 20, number
of leaves = 244.

Figure 2: Feature importance (sum of atomic features).

Figure 3: Feature importance (standard deviation of atomic features).

Next, we evaluate the models on the PDBbind core set. Table 1 summarizes our
main findings. The following metrics are reported: RMSE (root mean squared
error), MAE (mean absolute error), SD (standard deviation), R (Pearson’s
correlation). The metrics are averaged across 5 different runs (their standard
deviation are presented in parentheses). The best values are denoted in bold.
On the core set, LightGBM has the same RMSE as SIGN on average but it deviates
less, while CatBoost performs better in terms of SD and R. The only metric on
which SIGN outperforms GBDT model is MAE. Besides GBDT models, we also provide
the results for linear regression and three-layer perceptron.

To assess the generalization ability, the models are further evaluated on
CSAR-HiQ set. As can be seen from Table 1, classical GBDT models perform on
par with SIGN method. LightGBM demonstrates slight improvement over SIGN in
terms of RMSE and is able to outperform all presented models. Nevertheless,
SIGN shows better R score. It should be noted that GBDTs achieved these
results easily, owing to efficiency and scalability of their architecture.

### Feature importance.

Finally, we obtain the feature importance to analyze their contribution.
Figures 2 and 3 show the degree of contribution of atomic features (their sum
and standard deviation). Green and red bars correspond to ligand and protein
features respectively. Figure 1 illustrates the importance of interaction
features. It is notable that mainly uncommon interactions in the complexes,
such as $S-F$, $N-F$ have the highest importance in decision making of the
model. It is also important to note that the sum of atomic features has a more
even distribution of importance than those of interaction features while
having the same peak in magnitude. This might suggest that pooling atomic-
level features into graph-level features can significantly boost the
performance of RF-Score.

### Training time and memory consumption.

We also report the benefits of using GBDTs in terms of speed and memory
consumption. All experiments were conducted on 22 CPU cores and 1 NVIDIA A10
GPU (GPU was used for training SIGN). While training SIGN took 11 hours (300
epochs with batch size set to 16), training CatBoost took only 3 minutes.
Moreover, LightGBM has been trained only under 4 seconds in total. As for
memory consumption, SIGN allocates 19 GB of RAM and 13 GB of VRAM, whereas
GBDT models required less than 1 GB of RAM.

##  5 Discussion and Future Work

In this work, we looked back at traditional ML methods for predicting PL
binding affinity and proposed a way to improve their performance. We showed
that the inclusion of graph-level features with their standard deviation
results in better metrics. To our surprise, GBDT models can be comparable to
heavy GNN models without much effort. We hope that our findings would bring
attention back to simpler models, as there might be more undiscovered
potential. Next, we plan to integrate GBDT and GNN into one end-to-end
architecture to utilize the advantages of heterogeneous learning and
representation learning, as recently reported in our work, BGNN (Ivanov &
Prokhorenkova, 2021). The approach demonstrated state-of-the-art performance
and efficiency results on node-level prediction tasks, and adapting the
architecture to process the graph-level molecular features is a promising
direction.

## References

  * Akiba et al. (2019) Akiba, T., Sano, S., Yanase, T., Ohta, T., and Koyama, M.  Optuna: A next-generation hyperparameter optimization framework.  In _Proceedings of the 25rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , 2019.
  * Ballester & Mitchell (2010) Ballester, P. J. and Mitchell, J. B.  A machine learning approach to predicting protein–ligand binding affinity with applications to molecular docking.  _Bioinformatics_ , 26(9), 2010.
  * Behler & Parrinello (2007) Behler, J. and Parrinello, M.  Generalized neural-network representation of high-dimensional potential-energy surfaces.  _Phys. Rev. Lett._ , 98, 2007.
  * Dorogush et al. (2018) Dorogush, A. V., Ershov, V., and Gulin, A.  Catboost: gradient boosting with categorical features support.  _arXiv preprint arXiv:1810.11363_ , 2018.
  * Dunbar Jr et al. (2011) Dunbar Jr, J. B., Smith, R. D., Yang, C.-Y., Ung, P. M.-U., Lexa, K. W., Khazanov, N. A., Stuckey, J. A., Wang, S., and Carlson, H. A.  Csar benchmark exercise of 2010: selection of the protein–ligand complexes.  _Journal of chemical information and modeling_ , 51(9), 2011\.
  * Friedman (2001) Friedman, J. H.  Greedy function approximation: a gradient boosting machine.  _Annals of statistics_ , 2001.
  * Gorishniy et al. (2021) Gorishniy, Y., Rubachev, I., Khrulkov, V., and Babenko, A.  Revisiting deep learning models for tabular data.  _Advances in Neural Information Processing Systems_ , 34, 2021.
  * Gupta et al. (2020) Gupta, P., Gasse, M., Khalil, E., Mudigonda, P., Lodi, A., and Bengio, Y.  Hybrid models for learning to branch.  _Advances in neural information processing systems_ , 33, 2020.
  * Ivanov & Burnaev (2018) Ivanov, S. and Burnaev, E.  Anonymous walk embeddings.  In _International conference on machine learning_. PMLR, 2018.
  * Ivanov & Prokhorenkova (2021) Ivanov, S. and Prokhorenkova, L.  Boost then convolve: Gradient boosting meets graph neural networks.  _arXiv preprint arXiv:2101.08543_ , 2021.
  * Ivanov et al. (2018) Ivanov, S., Durasov, N., and Burnaev, E.  Learning node embeddings for influence set completion.  In _2018 IEEE International Conference on Data Mining Workshops (ICDMW)_. IEEE, 2018.
  * Jiang et al. (2021) Jiang, D., Hsieh, C.-Y., Wu, Z., Kang, Y., Wang, J., Wang, E., Liao, B., Shen, C., Xu, L., Wu, J., et al.  Interactiongraphnet: A novel and efficient deep graph representation learning framework for accurate protein–ligand interaction predictions.  _Journal of medicinal chemistry_ , 64(24), 2021.
  * Jones et al. (2021) Jones, D., Kim, H., Zhang, X., Zemla, A., Stevenson, G., Bennett, W. D., Kirshner, D., Wong, S. E., Lightstone, F. C., and Allen, J. E.  Improved protein–ligand binding affinity prediction with structure-based deep fusion inference.  _Journal of chemical information and modeling_ , 61(4), 2021\.
  * Jumper et al. (2021) Jumper, J., Evans, R., Pritzel, A., Green, T., Figurnov, M., Ronneberger, O., Tunyasuvunakool, K., Bates, R., Žídek, A., Potapenko, A., et al.  Highly accurate protein structure prediction with alphafold.  _Nature_ , 596(7873), 2021.
  * Karlov et al. (2020) Karlov, D. S., Sosnin, S., Fedorov, M. V., and Popov, P.  graphdelta: Mpnn scoring function for the affinity prediction of protein–ligand complexes.  _ACS omega_ , 5(10), 2020.
  * Ke et al. (2017) Ke, G., Meng, Q., Finley, T., Wang, T., Chen, W., Ma, W., Ye, Q., and Liu, T.-Y.  Lightgbm: A highly efficient gradient boosting decision tree.  _Advances in neural information processing systems_ , 30, 2017.
  * Kearnes et al. (2016) Kearnes, S., McCloskey, K., Berndl, M., Pande, V., and Riley, P.  Molecular graph convolutions: moving beyond fingerprints.  _Journal of computer-aided molecular design_ , 30(8), 2016\.
  * Kitchen et al. (2004) Kitchen, D. B., Decornez, H., Furr, J. R., and Bajorath, J.  Docking and scoring in virtual screening for drug discovery: methods and applications.  _Nature Reviews Drug Discovery_ , 3(11), 2004.
  * Kriege et al. (2020) Kriege, N. M., Johansson, F. D., and Morris, C.  A survey on graph kernels.  _Applied Network Science_ , 5(1), 2020.
  * Li et al. (2021) Li, S., Zhou, J., Xu, T., Huang, L., Wang, F., Xiong, H., Huang, W., Dou, D., and Xiong, H.  Structure-aware interactive graph neural networks for the prediction of protein-ligand binding affinity.  In _Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery & Data Mining_, 2021.
  * Lim et al. (2019) Lim, J., Ryu, S., Park, K., Choe, Y. J., Ham, J., and Kim, W. Y.  Predicting drug–target interaction using a novel graph neural network with 3d structure-embedded graph representation.  _Journal of Chemical Information and Modeling_ , 59(9), 2019\.
  * Liu et al. (2020) Liu, S., Park, J. H., and Yoo, S.  Efficient and effective graph convolution networks.  In _Proceedings of the 2020 SIAM International Conference on Data Mining_. SIAM, 2020.
  * Liu et al. (2015) Liu, Z., Li, Y., Han, L., Li, J., Liu, J., Zhao, Z., Nie, W., Liu, Y., and Wang, R.  Pdb-wide collection of binding data: current status of the pdbbind database.  _Bioinformatics_ , 31(3), 2015.
  * Meli et al. (2021) Meli, R., Anighoro, A., Bodkin, M. J., Morris, G. M., and Biggin, P. C.  Learning protein-ligand binding affinity with atomic environment vectors.  _Journal of Cheminformatics_ , 13(1), 2021.
  * Moon et al. (2022) Moon, S., Zhung, W., Yang, S., Lim, J., and Kim, W. Y.  Pignet: a physics-informed deep learning model toward generalized drug–target interaction predictions.  _Chemical Science_ , 13(13), 2022.
  * Qin et al. (2021) Qin, Z., Yan, L., Zhuang, H., Tay, Y., Pasumarthi, R. K., Wang, X., Bendersky, M., and Najork, M.  Are neural rankers still outperformed by gradient boosted decision trees?  In _International Conference on Learning Representations (ICLR)_ , 2021.
  * Sharaev et al. (2018) Sharaev, M., Artemov, A., Kondrateva, E., Ivanov, S., Sushchinskaya, S., Bernstein, A., Cichocki, A., and Burnaev, E.  Learning connectivity patterns via graph kernels for fmri-based depression diagnostics.  In _2018 IEEE International Conference on Data Mining Workshops (ICDMW)_. IEEE, 2018.
  * Son & Kim (2021) Son, J. and Kim, D.  Development of a graph convolutional neural network model for efficient prediction of protein-ligand binding affinities.  _PLoS ONE_ , 16(4), 2021.
  * Su et al. (2018) Su, M., Yang, Q., Du, Y., Feng, G., Liu, Z., Li, Y., and Wang, R.  Comparative assessment of scoring functions: the casf-2016 update.  _Journal of chemical information and modeling_ , 59(2), 2018\.
  * Tiezzi et al. (2020) Tiezzi, M., Marra, G., Melacci, S., Maggini, M., and Gori, M.  A lagrangian approach to information propagation in graph neural networks.  _arXiv preprint arXiv:2002.07684_ , 2020.
  * Wan et al. (2020) Wan, S., Bhati, A. P., Zasada, S. J., and Coveney, P. V.  Rapid, accurate, precise and reproducible ligand–protein binding free energy prediction.  _Interface Focus_ , 10(20200007), 2020.
  * Wang & Xia (2012) Wang, J. and Xia, Y.  Fast graph construction using auction algorithm.  _arXiv preprint arXiv:1210.4917_ , 2012.