# Non-convexity of extremal length

Nathaniel Sagman  Nathaniel Sagman: University of Luxembourg, 2 Av. de
l’Universite, 4365 Esch-sur-Alzette, Luxembourg<EMAIL_ADDRESS>
###### Abstract.

With respect to every Riemannian metric, the Teichmüller metric, and the
Thurston metric on Teichmüller space, we show that there exist measured
foliations on surfaces whose extremal length functions are not convex. The
construction uses harmonic maps to $\mathbb{R}$-trees and minimal surfaces in
$\mathbb{R}^{n}.$

##  1\. Introduction

Let $\Sigma_{g}$ be a closed oriented surface of genus $g\geq 2,$ and let
$\mathbf{T}_{g}$ be the Teichmüller space of marked Riemann surface structures
on $\Sigma_{g}.$ To any measured foliation $\mathcal{F}$ on $\Sigma_{g}$ we
can associate the extremal length function
$\mathbf{EL}_{\mathcal{F}}:\mathbf{T}_{g}\to(0,\infty).$ Extremal length
functions play a large role in Teichmüller theory. See, for instance,
Kerckhoff’s formula [8, Theorem 4] and the Gardiner-Masur compactification
[6].

Liu-Su proved that $\mathbf{EL}_{\mathcal{F}}$ is plurisubharmonic, and
Miyachi proved the stronger result that it is log-plurisubharmonic (see [11]
and [17]). Note that convexity with respect to a Riemannian metric implies
plurisubharmonicity. Rafi-Lenhzen proved that, on Teichmüller geodesics,
extremal length is $K$-quasi-convex, but they also constructed a Teichmüller
geodesic along which the extremal length is not convex [10]. Continuing in
this direction, Bourque-Rafi proved that the Teichmüller metric admits non-
convex balls by finding foliations and geodesics where the extremal length is
not convex under any reparametrization [3] (see especially Lemma 1.2 in [3]).

In this note, we extend the non-convexity result of Rafi-Lenhzen [10]. Let
$\mathcal{C}$ denote the class of (possibly asymmetric) Finsler metrics on
$\mathbf{T}_{g}$ such that for every point $S$ in $\mathbf{T}_{g}$ and every
tangent vector $\mu$ at that point, there is a $C^{2}$ geodesic starting at
$S$ and tangent to $\mu$ at time zero. $\mathcal{C}$ includes every Riemanian
metric, notably the Weil-Petersson metric, but also the Teichmüller metric and
the Thurston metric. Rafi-Lenhzen build an explicit foliation and a
Teichmüller ray that has pieces along which the slope of the extremal length
function decreases. In contrast, we show that convexity fails at an
infinitesimal level.

######  Theorem 1.1.

For all $g\geq 2$ and $m\in\mathcal{C}$, there exists a measured foliation
$\mathcal{F}$ on $\Sigma_{g}$ with real analytic extremal length function and
a geodesic $t\mapsto S_{t}$ for $m$ with the property that

$\frac{d^{2}}{dt^{2}}|_{t=0}\mathbf{EL}_{\mathcal{F}}(S_{t})<0.$

In particular, $\mathbf{EL}_{\mathcal{F}}$ is not convex with respect to $m$.

As noted in [13], it follows from the main result of [13] that with respect to
every Riemannian metric on $\mathbf{T}_{g}$, the energy functional for
harmonic maps associated with a Fuchsian representation can be non-convex. By
the paper [23], the same result holds for (non-Fuchsian) Hitchin
representations. We prove Theorem 1.1 by interpreting extremal length as an
energy. Drawing from recent work on minimal surfaces (see
[12],[13],[14],[15],[23]), we establish a link between non-convexity of
extremal length and instability of minimal surfaces in $\mathbb{R}^{n}$.

One of the main takeaways of the proof is that a destabilizing variation of an
equivariant minimal surface in $\mathbb{R}^{n}$ produces a foliation (or even
a number of foliations) whose extremal length can be lowered to second order.
And although it’s probably difficult in practice, if one has the explicit
minimal surface data, then one could compute quantities associated with the
extremal length (see Remark 3.13).

It would require some care, but one could try to use minimal surfaces to
construct a foliation and a geodesic (for some metric) such that, in
restriction to the geodesic, the extremal length has a local maximum at time
zero. This would imply that the extremal length is not convexoidal for the
metric. In fact, it is conjectured in [2] that the extremal length systole
attains a local maximum at the regular octahedron punctured at its vertices,
which would imply that Voronoï’s criterion fails for the extremal length
systole, and moreover that extremal length is not convexoidal for any metric
in $\mathcal{C}$ (see [1, Definition 1.4 and Proposition 1.5] for definitions
and justification).

Finally, let us remark that a number of questions remain open related to
convexity in Teichmüller geometry. It is not known if the Teichmüller metric
convex hull of $3$ points in $\mathbf{T}_{g}$ can be all of $\mathbf{T}_{g}$.
While sufficiently small Teichmüller balls are always convex (the analogous
fact holds for any Finsler metric), it is unclear if sets of the form
$\\{S\in\mathbf{T}_{g}:\mathbf{EL}_{\mathcal{F}}(S)<\alpha\\}$, referred to as
horoballs in [3], are convex for $\alpha$ small. Our proof of Theorem 1.1
suggests a new way to probe the convexity question for such horoballs.

###  1.1. Acknowledgements

I’d like to thank Kasra Rafi and Maxime Fortier Bourque for discussion on this
topic. I am funded by the FNR grant O20/14766753, Convex Surfaces in
Hyperbolic Geometry.

##  2\. Preliminaries

###  2.1. Measured foliations

Let $\mathcal{S}$ be the set of non-trivial homotopy classes of simple closed
curves on $\Sigma_{g}$, and $\mathbb{R}^{\mathcal{S}}$ the product space with
the weak topology. Any $\gamma\in\mathcal{S}$ determines a point in
$\mathbb{R}^{\mathcal{S}}$ through the intersection number,

(1) $\gamma\mapsto(i(\gamma,\alpha))_{\alpha\in\mathcal{S}}.$

A weighted multicurve is a formal positive linear combination of classes in
$\mathcal{S}$. We extend the intersection number to the space of weighted
multicurves $\mathcal{WS}$ by

$i\Big{(}\sum_{j=1}^{n}a_{j}\gamma_{j},\sum_{k=1}^{m}b_{k}\alpha_{k}\Big{)}=\sum_{j=1}^{n}\sum_{k=1}^{m}a_{j}b_{k}i(\gamma_{j},\alpha_{k}),$

which as above yields an embedding from $\mathcal{WS}$ into
$\mathbb{R}^{\mathcal{S}}$ via the same map (1). To us, the space of measured
foliations $\mathcal{MF}$ is the closure of $\mathcal{WS}$ in
$\mathbb{R}^{\mathcal{S}}$. Note that the intersection number extends
continuously to $\mathcal{MF}$ [5].

Alternatively, a measured foliation $\mathcal{F}$ is a singular foliation on
$\Sigma_{g}$, the singularities being $k$-prongs, $k\geq 3$, equipped with a
transverse measure: an absolutely continuous measure defined on arcs
transverse to the foliation and which is invariant under leaf-preserving
isotopy. Two measured foliations are measure equivalent if they differ by a
leaf-preserving isotopy and Whitehead moves. See [5, Exposé 5] for the precise
definitions. The intersection function is defined on simple closed curves by
integration against the transverse measure.

Let $S$ be a Riemann surface structure on $\Sigma_{g}$. The vertical (resp.
horizontal) foliation of a holomorphic quadratic differential $\phi$ on $S$ is
the singular foliation whose leaves are the integral curves of the line field
on $S\backslash\phi^{-1}(0)$ on which $\phi$ is a negative (resp. positive)
real number. The singularities are indeed prongs at the zeros, with a zero of
order $k$ corresponding to a prong with $k+2$ segments. Both foliations come
with transverse measures determined by $|\textrm{Re}\sqrt{\phi}|$ and
$|\textrm{Im}\sqrt{\phi}|$ respectively. In this paper, we will always use the
vertical foliation.

The Hubbard-Masur theorem asserts that on the given Riemann surface $S$, every
measured foliation $\mathcal{F}$ is measure equivalent to one arising from the
construction above [7]. We refer to the corresponding differential $\phi$ as
the Hubbard-Masur differential.

###  2.2. Extremal length (and its regularity)

Let $S$ be a Riemann surface structure on $\Sigma_{g},$ and $A\subset S$ a
doubly connected domain, conformally equivalent to an annulus
$\\{z\in\operatorname{\mathbb{C}}:1<|z|<R\\}.$ The modulus of $A$ is the
quantity

$\textrm{Mod}(A)=\frac{1}{2\pi}\log R.$

######  Definition 2.1.

The extremal length of a homotopically non-trivial simple closed curve
$\gamma$ with respect to $S$ is

$\textrm{EL}(S,\gamma)=\inf_{A}\textrm{Mod}(A),$

where the infimum is taken over all doubly connected domains $A$ homotopic to
$\gamma.$

Given a weighted multicurve $\gamma=\sum_{j=1}^{n}a_{j}\gamma_{j},$ we define

$\textrm{EL}(S,\gamma)=\sum_{j=1}^{n}a_{j}^{2}\textrm{EL}(S,\gamma_{j}).$

Kerckhoff showed that the map $\textrm{EL}(S,\cdot)$ extends continuously to
all measured foliations, defining a map
$\textrm{EL}(S,\cdot):\mathcal{MF}\to(0,\infty)$ [8].

Fix a measured foliation $\mathcal{F}$ on $\Sigma_{g}.$ We define the extremal
length function on Teichmüller space,
$\mathbf{EL}_{\mathcal{F}}:\mathbf{T}_{g}\to(0,\infty)$, by

$\mathbf{EL}_{\mathcal{F}}(S)=\textrm{EL}(S,\mathcal{F}).$

In terms of the Hubbard-Masur differential $\phi$, the extremal length is the
$L^{1}$ norm:

$\mathbf{EL}_{\mathcal{F}}(S)=\int_{S}|\phi|.$

Recall that the tangent space of $\mathbf{T}_{g}$ at a surface $S$ identifies
with the vector space of harmonic Beltrami forms on $S$. By direct
computation, $\mathbf{EL}_{\mathcal{F}}$ is $C^{1}$, and the derivative is
given by

$d(\mathbf{EL}_{\mathcal{F}})_{S}(\mu)=-4\textrm{Re}\int\phi\mu.$

From our understanding, it is unknown if $\mathbf{EL}_{\mathcal{F}}$ is
$C^{2}$, and we’ll have to address this point in the main proof. Royden’s
computation in [22, Lemma 1] seems relevant to this problem, and the papers
[19] and [20] suggest that it is at most $C^{2}.$ Around a point in
$\mathbf{T}_{g}$ where all zeros of the Hubbard-Masur differential are simple,
$\mathbf{EL}_{\mathcal{F}}$ is real analytic (see [16]). This condition is
generic, and guaranteed when $\mathcal{F}$ has only $3$-pronged singularities
and no saddle connections. We do not pursue the general regularity question in
the current paper.

###  2.3. Harmonic maps

We plan to interpret extremal length in terms of harmonic maps to
$\mathbb{R}$-trees. As above, let $S$ be a Riemann surface structure on
$\Sigma_{g}$ and $\nu$ a smooth metric that is conformal with respect to the
complex structure. Let $(M,d)$ be a complete and non-positively curved (NPC)
length space equipped with an action
$\rho:\pi_{1}(\Sigma_{g})\to\textrm{Isom}(M,d).$ Let $\tilde{S}$ be the
universal cover and $h:\tilde{S}\to(M,d)$ a $\rho$-equivariant and Lipschitz
map. Korevaar-Schoen [9, Theorem 2.3.2] associate a locally $L^{1}$ measurable
metric $g=g(h)$, defined on pairs of Lipschitz vector fields. If $h$ is a
$C^{1}$ map to a smooth Riemannian manifold $(M,\sigma)$, and the distance $d$
is induced by a Riemannian metric $\sigma$, then $g(h)$ is represented by the
pullback metric $h^{*}\sigma$. Since $\rho$ is acting by isometries, the
tensor $g(h)$ descends to $S$. Henceforth we consider it a function on $S$.
The energy density is the locally $L^{1}$ function

$e(h)=\frac{1}{2}\textrm{trace}_{\nu}g(h).$

The total energy is

$\mathcal{E}(S,h)=\int_{S}e(h)dA,$

where $dA$ is the area form of $\nu$. The measurable $2$-form $e(h)dA$ does
not depend on the choice of compatible metric $\nu$, but only on the complex
structure.

######  Definition 2.2.

$h$ is harmonic if it is a critical point for the energy
$h\mapsto\mathcal{E}(S,h)$.

Let $g_{ij}(h)$ be the components of $g(h)$ in a holomorphic local coordinate
$z=x_{1}+ix_{2}$. The Hopf differential of $h$ is the measurable tensor on $S$
given in the local coordinate by

(2) $\phi(h)(z)=\frac{1}{4}(g_{11}(h)(z)-g_{22}(h)(z)-2ig_{12}(h)(z))dz^{2}.$

In the Riemannian setting, (2) is

$\phi(h)(z)=h^{*}\sigma\Big{(}\frac{\partial}{\partial
z},\frac{\partial}{\partial z}\Big{)}(z)dz^{2}.$

When $h$ is harmonic, even in the metric space setting, the Hopf differential
is represented by a holomorphic quadratic differential.

Assume that $\rho$ has the following property: for any Riemann surface $S$
representing a point in $\mathbf{T}_{g}$, there is a unique $\rho$-equivariant
harmonic map $h:\tilde{S}\to(M,d)$. The energy functional on Teichmüller space
$\mathbf{E}_{\rho}:\mathbf{T}_{g}\to[0,\infty)$ is defined by

$\mathbf{E}_{\rho}(S)=\mathcal{E}(S,h).$

When $\mathbf{E}_{\rho}$ is $C^{1}$ and the associated harmonic map has Hopf
differential $\phi$, the derivative in the direction of a harmonic Beltrami
form $\mu$ is

(3) $d(\mathbf{E}_{\rho})_{S}(\mu)=-4\textrm{Re}\int\phi\mu.$

See [25] for the proof in the metric space context.

###  2.4. $\mathbb{R}$-trees dual to foliations

######  Definition 2.3.

An $\mathbb{R}$-tree is a length space $(T,d)$ such that any two points are
connected by a unique arc, and every arc is a geodesic, isometric to a segment
in $\mathbb{R}$.

It is a basic exercise to show that $\mathbb{R}$-trees are complete and NPC
length spaces. We concern ourselves with a particular class of actions on
$\mathbb{R}$-trees, obtained as follows. Let $\mathcal{F}$ be a measured
foliation on $S$ with Hubbard-Masur differential $\phi$. Lifting $\mathcal{F}$
to the universal cover $\tilde{S}$, we define an equivalence relation on
$\tilde{S}$ by $x\sim y$ if $x$ and $y$ lie on the same leaf. The quotient
space $\tilde{S}/\sim$ is denoted $T$. Pushing the transverse measure down via
the projection $\pi:\tilde{S}\to T$ yields a distance function $d$ that turns
$(T,d)$ into an $\mathbb{R}$-tree, with an induced action
$\rho:\pi_{1}(\Sigma_{g})\to\textrm{Isom}(T,d).$ Under this distance, the
projection map $\pi:\tilde{S}\to(T,d)$ is $\rho$-equivariant and harmonic, and
the Hopf differential is exactly $\phi/4$ (see [26, Section 3]).

The energy density of $\pi$ can be described explicitly: at a point
$p\in\tilde{S}$ on which $\phi(p)\neq 0$, $\pi$ locally isometrically factors
through a segment in $\mathbb{R}$. In a small neighbourhood around that point,
$g(h)$ is represented by the pullback metric of the locally defined map to
$\mathbb{R}$. Therefore, it can be computed that

(4) $e(\pi)=\nu^{-1}|\phi|/2.$

Similarly, this provides one way to compute $\phi(\pi)=\phi/4$. In view of
(4), we will always rescale the metric on $T$ from $(T,d)$ to $(T,2d).$ In
this normalization, the total energy is

(5) $\mathcal{E}(S,\pi)=\int_{S}|\phi|.$

Keeping $\rho$ and varying the source Riemann surface, $\rho$-equivariant
harmonic maps always exist and are unique [27], and hence there is an energy
functional $\mathbf{E}_{\rho}$. From the formula (5), we deduce the following.

######  Proposition 2.4.

Let $\mathcal{F}$ be a measured foliation with Hubbard-Masur differential
$\phi$, and $\rho$ the action on the $\mathbb{R}$-tree dual to $\mathcal{F}$.
As functions on $\mathbf{T}_{g},$
$\mathbf{E}_{\rho}=\mathbf{EL}_{\mathcal{F}}.$

Accordingly, the same discussion on regularity from Section 2.2 applies to
$\mathbf{E}_{\rho}.$ From this point on, we will think about extremal length
solely in terms of harmonic maps to $\mathbb{R}$-trees.

##  3\. Non-convexity

Let $m$ be a metric distance function on $\mathbf{T}_{g}$ in which
$(\mathbf{T}_{g},m)$ is a length space. We say that a function
$F:\mathbf{T}_{g}\to\mathbb{R}$ is convex with respect to $m$ if for all
geodesics $c:[0,1]\to\mathbf{T}_{g}$, the function $F\circ
c:[0,1]\to\mathbb{R}$ is convex. If $c$ is $C^{2}$ and $F$ is $C^{2}$ around
the image of $c$, then $F\circ c$ is convex if and only if the second
derivative is non-negative at all points.

After discussing metrics on $\mathbf{T}_{g}$ in 3.1, we recall constructions
from [15] relating variations of minimal surfaces in $\mathbb{R}^{n}$ to
variations of minimal maps to products of $\mathbb{R}$-trees. We will then use
minimal surfaces in $\mathbb{R}^{n}$ to find a harmonic map to an
$\mathbb{R}$-tree (or, a measured foliation) whose energy (extremal length)
can be lowered to second order.

###  3.1. Metrics on $\mathbf{T}_{g}$

Recall the class of metrics $\mathcal{C}$ from the introduction. Let’s briefly
justify that the Teichmüller metric $d_{T}$ and the Thurston metric $d_{Th}$
are contained in $\mathcal{C}.$ This may follow from a general theory of
asymmetric Finsler metrics with certain properties, but we couldn’t find a
source and we prefer to be hands-on.

######  Definition 3.1.

The Teichmüller metric is defined by $d_{T}(S,S^{\prime})=\inf_{g}\log K(g),$
where $K(g)$ is the maximum quasiconformal dilatation of a quasiconformal map
$g:S\to S^{\prime}$.

Recall that at a Riemann surface $S$, $T_{S}\mathbf{T}_{g}$ identifies with
the space of harmonic Beltrami forms on $S$. If $\mu$ is any such Beltrami
form, away from the finite zero set of $\mu$ we can locally choose a
coordinate $z=x+iy$ in which $\mu=d\overline{z}/dz.$ The Teichmüller mapping
in the direction of $\mu$ at scale $K$ is defined in such a coordinate by

(6) $f_{\mu,K}(x,y)=K^{1/2}x+K^{-1/2}y$

We define $f_{\mu,K}$ globally by doing (6) over the local patches, and
extending to all of $S$ by continuity. The Teichmüller ray $K\mapsto
f_{\mu,K}$ is a geodesic for $d_{T}$ tangent to $\mu$ at $K=0.$

######  Definition 3.2.

The Thurston metric is defined by
$d_{Th}(S,S^{\prime})=\inf_{g}\log\textrm{Lip}(g),$ where $\textrm{Lip}(g)$ is
the Lipschitz constant of a Lipschitz map $g$ taking $S$ to $S^{\prime}$.

We couldn’t find a clean statement in the literature about the existence of
Thurston geodesics in a given tangent direction. The result can probably be
established through the constructions of Thurston’s original paper [24], but
for the sake of brevity we’ll instead cite the recent work of Pan and Wolf
[18].

For any Riemann surface $S$ and projective measured lamination on $S$, Pan and
Wolf construct a “harmonic stretch line,” which is a Thurston geodesic that in
some sense solves an energy-minimization problem. Every unit tangent vector to
$\mathbf{T}_{g}$ at $S$ is tangent to a harmonic stretch line [18, Remark
1.12]. Moreover, harmonic stretch lines are special examples of “piecewise
harmonic stretch lines,” which are, as stated in Theorem 1.7 of [18], real
analytic paths in Teichmüller space.

###  3.2. Harmonic functions

Let $S$ be a Riemann surface structure on $\Sigma_{g}$ and $\phi$ a
holomorphic quadratic differential. Let $(T,2d)$ be the dual $\mathbb{R}$-tree
with action $\rho$ and harmonic projection map $\pi:\tilde{S}\to(T,2d)$.
Assume that $\phi$ is the square of an abelian differential $\alpha$. The
cohomology class of the harmonic $1$-form $\textrm{Re}\alpha$ determines a
representation $\chi:\pi_{1}(\Sigma_{g})\to(\mathbb{R},+),$ and integrating
from a basepoint $p\in\tilde{S}$ yields a $\chi$-equivariant harmonic function

$h:\tilde{S}\to\mathbb{R},\hskip
2.84526ptf_{i}(z)=\int_{p}^{z}\textrm{Re}\alpha.$

We can compute directly that $\phi(h)=\phi$ and that for any choice of
conformal metric on $S$, $e(h)=e(\pi)$. Geometrically, $h$ is related to $\pi$
through the folding map $p$, which is a map $p:(T,d)\to\mathbb{R}$ satisfying
$h=p\circ\pi$ and restricting to an isometry on geodesic segments of $(T,d)$
(see [14, Section 4]). For any Riemann surface, $\chi$-equivariant harmonic
functions exist, but they are unique only up to translations of
$\mathbb{R}^{n}$. Nevertheless, the energy density is independent of the
choice of harmonic function, and it is possible to choose the harmonic
functions locally to vary real analytically with the choice of Riemann surface
(see [4, Section 5]). Thus, as in 2.3, we may define a (real analytic) energy
functional $\mathbf{E}_{\chi}$ on $\mathbf{T}_{g}.$

In general, there is a degree $2$ branched covering $\tau:C\to S$ in which
$\phi$ lifts to a square $\tilde{\phi}=\alpha^{2}$, disconnected if $\phi$ is
already a square, and which has the universal property that any other branched
cover on which $\phi$ lifts to a square must factor through $\tau$. We can
repeat the above construction on $C,$ obtaining an equivariant harmonic map to
an $\mathbb{R}$-tree $\pi^{\prime}:\tilde{C}\to(T^{\prime},2d^{\prime})$ that
folds onto an equivariant harmonic function $h$ from $\tilde{C}$ to
$\mathbb{R}$. $C$ comes with a holomorphic involution that negates $\alpha$.
The involution leaves the energy densities of $\pi^{\prime}$ and $f$
invariant, so that they descend all the way to $S$, where they agree with that
of the map to the original $\mathbb{R}$-tree.

###  3.3. Minimal maps

Classically, a minimal map to a Riemannian manifold is a harmonic and
conformal immersion. An immersion is conformal precisely when the Hopf
differential vanishes identically. For an NPC space $(M,d)$, we make the
following definition.

######  Definition 3.3.

$h:\tilde{S}\to(M,d)$ is minimal if it is harmonic and $\phi(h)=0.$

In the presence of a $C^{1}$ energy functional $\mathbf{E}_{\rho},$ by (3),
$S$ is minimal if and only if $S$ is a critical point of $\mathbf{E}_{\rho}.$
For equivariant maps to $\mathbb{R}^{n}$, minimal maps are also critical
points of the area functional

$f\mapsto A(f)=\int_{\Sigma_{g}}dA_{f},$

where $dA_{f}$ is the area form of the pullback of the Euclidean metric by
$f$. We record the following consequence of the definitions. Let $(X,d)$ be a
product of NPC spaces $(M_{i},d_{i}).$

######  Proposition 3.4.

$h:\tilde{S}\to(X,d)$ is harmonic if and only if every component map
$h_{i}:\tilde{S}\to(M_{i},d_{i})$ is harmonic. Moreover,
$\phi(h)=\sum_{i=1}^{n}\phi(h_{i}).$

Let $\phi_{1},\dots,\phi_{n}$ be holomorphic quadratic differentials on $S$
and let $(M,d)$ be the product of the dual $\mathbb{R}$-trees, with product
action $\rho$ and product of projection maps
$\pi=(\pi_{1},\dots,\pi_{n}):\tilde{S}\to(M,d).$ The energy functional
$\mathbf{E}_{\rho}$ is the sum of the component energy functionals. Similar to
the previous subsection, there is a degree $2^{n}$ branched covering
$\tau:C\to S$ on which each $\phi_{i}$ lifts to a square, with the analogous
universal property, and which comes with $n$ commuting holomorphic involutions
that each negate a $1$-form. There is a product representation
$\chi=(\chi_{1},\dots,\chi_{n}):\pi_{1}(C)\to(\mathbb{R}^{n},+)$ with energy
functional $\mathbf{E}_{\chi}$ and equivariant harmonic function
$h=(h_{1},\dots,h_{n}):\tilde{C}\to\mathbb{R}^{n}$. Proposition 3.4 gives the
observation below.

######  Proposition 3.5.

$h$ is minimal if and only if $\pi$ is minimal if and only if
$\sum_{i=1}^{n}\phi_{i}=0.$

Let us thus assume the $\phi_{i}$’s sum to zero. The main input toward Theorem
1.1 is Proposition 3.6 below, which is used to turn variations of $h$ into
variations of $\pi$. We will need to restrict to a class of variations. Set
$\textrm{Var}_{\tau}(h)$ to be the space of $C^{\infty}$ functions
$\dot{h}:\tilde{C}\to\mathbb{R}^{n}$ that are invariant under $\pi_{1}(C)$ and
the lifts to $\tilde{C}$ of the $n$ holomorphic involutions. Of course, such
$\dot{h}$ is equivalent to a function on $S$.

######  Proposition 3.6 (Propositions 4.4, 4.6, and 5.1 of [15]).

Let $\dot{h}\in\textrm{Var}_{\tau}(h)$. For every $\epsilon>0$, there exists a
$C^{\infty}$ path of Riemann surfaces $t\mapsto C_{t}$ and $C^{\infty}$ paths
of $C^{\infty}$ maps $t\mapsto f_{i}^{t}:C\to C_{t}$ starting at the identity
such that

(7)
$\frac{d^{2}}{dt^{2}}|_{t=0}\sum_{i=1}^{n}\mathbf{E}_{\chi_{i}}(C_{t},h_{i}\circ(\tilde{f}_{i}^{t})^{-1})\leq\frac{d^{2}}{dt^{2}}|_{t=0}A(h_{t})+\epsilon,$

where $\tilde{f}_{i}^{t}$ is the lift to $\tilde{C}$. The Riemann surfaces
$C_{t}$ descend through the branched cover $\tau$ to Riemann surface
structures $S_{t}$ on $S$. Similarly, the $f_{i}^{t}$’s descend to
$f_{i}^{t}:S\to S_{t}.$

Very briefly: without perturbing the area too much, one can modify $\dot{h}$
to be zero in a neighbourhood of the zeros of the $\phi_{i}$’s. There is a
canonical way to pull such a variation back to $n$ vector fields on the
surface $\tilde{C}$, which then generate flows $t\mapsto f_{i}^{t}$. With
respect to the conformal structure of the image of
$h_{t}=(h_{1}\circ(f_{1}^{t})^{-1},\dots,h_{n}\circ(f_{n}^{t})^{-1})$, which
we label $C_{t},$ the energy of $h_{t}$ is equal to its area. By invariance
properties of $\dot{h},$ everything can be chosen to descend to $\Sigma_{g}$.
Such a self-maps variation gives a variation of $\pi$,

(8)
$\pi_{t}=(\pi_{1}\circ(f_{1}^{t})^{-1},\dots,\pi_{n}\circ(f_{n}^{t})^{-1}).$

If $\tilde{f}_{i}^{t}$ also denotes the lift to $\tilde{S}$, then we can see
by the local isometric factoring described in 2.4, or the folding map of 3.2,

(9)
$e(\pi_{i}\circ(\tilde{f}_{i}^{t})^{-1})=e(h_{i}\circ(\tilde{f}_{i}^{t})^{-1})$

(note that both densities descend to $S$, which is where the equation (9) is
defined).

Finally, we will want destabilizing variations.

######  Definition 3.7.

$h$ is $\tau$-unstable if there exists $\dot{h}\in\textrm{Var}_{\tau}(h)$ such
that

$\frac{d^{2}}{dt^{2}}|_{t=0}A(h+t\dot{h})<0.$

For all $g\geq 3$ and $n\geq 3$, we can take $C=S$: one can find abelian
differentials on $S$ whose squares sum to $0$ that give an unstable minimal
map. One way to produce such a map is to lift an unstable minimal surface in
the $3$-torus to the universal covers. For $g=2$, it turns out that any
equivariant minimal surface is stable. However, we proved

######  Theorem 3.8 (Section 5.3 in [15]).

There exists a Riemann surface $S$ of genus $2$ with
$\phi_{1},\dots,\phi_{n},$ $n\geq 3,$ summing to $0$ and which give a non-
trivial branched cover $\tau:C\to S$ and a $\tau$-unstable minimal map
$h:\tilde{C}\to\mathbb{R}^{n}.$

###  3.4. Proof of Theorem 1.1

Resuming the setup from Section 3.3, assume that $h$ is destabilized by some
$\dot{h}\in\textrm{Var}_{\tau}(h)$. Choose $\epsilon>0$ small enough so that

$\frac{d^{2}}{dt^{2}}|_{t=0}A(h+t\dot{h})+\epsilon<0.$

Applying Proposition 3.6, defining $\pi_{t}$ as in (8) and using (9), we
arrive at the following.

######  Proposition 3.9.

There exist $C^{\infty}$ paths of $C^{\infty}$ maps to Riemann surfaces
$t\mapsto f_{i}^{t}:S\to S^{t}$ starting at the identity such that

(10)
$\frac{d^{2}}{dt^{2}}|_{t=0}\sum_{i=1}^{n}\mathcal{E}(S_{t},\pi_{i}\circ(f_{i}^{t})^{-1})\leq
2^{-n}\frac{d^{2}}{dt^{2}}|_{t=0}A(h+t\dot{h})+\epsilon<0.$

All recorded examples of $\tau$-unstable maps, in particular examples from
Theorem 3.8, come from differentials with even order zeros. Recalling our
regularity concerns from Section 2.2, we need the Proposition below. Say that
a holomorphic quadratic differential $\phi$ is generic if it has only simple
zeros. Generic quadratic differentials on $S$ form an open and dense subset.

######  Proposition 3.10.

For all $g\geq 2,$ $n\geq 3,$ we can choose generic holomorphic quadratic
differentials $\phi_{i}$ that give maps to $\mathbb{R}$-trees $\pi_{i},$ and
$C^{\infty}$ paths $t\mapsto f_{i}^{t}:S\to S_{t}$ starting at the identity
such that (10) holds:

$\frac{d^{2}}{dt^{2}}|_{t=0}\sum_{i=1}^{n}\mathcal{E}(S_{t},\pi_{i}\circ(f_{i}^{t})^{-1})<0.$

Morally, we’re using that instability in the sense of (10) is an open
property. To formalize the argument, we borrow a formula from Reich-Strebel.

######  Proposition 3.11 (Equation 1.1 in [21] and Proposition 3.1 in [15]).

Let $\pi:\tilde{S}\to(T,2d)$ be any equivariant map to an $\mathbb{R}$-tree
with Hopf differential $\phi$ and let $f:S\to S^{\prime}$ be any
quasiconformal map to another Riemann surface $S^{\prime}$ with lift
$\tilde{f}$ to $\tilde{S}$ and Beltrami form $\mu$.

(11)
$\mathcal{E}(S,\pi\circ\tilde{f}^{-1})-\mathcal{E}(S,\pi)=-4\textrm{Re}\int_{S}\phi\cdot\frac{\mu}{1-|\mu|^{2}}dxdy+4\int_{S}|\phi_{i}|\cdot\frac{|\mu|^{2}}{1-|\mu|^{2}}dxdy.$

###### Proof of Proposition 3.10.

Begin with the data $\phi_{i},f_{i}^{t}$ from Proposition 3.9. The
$\phi_{i}$’s may not be generic, but we know that (10) holds. Let
$\mu_{i}^{t}$ be the Beltrami form of $f_{i}^{t}$, and $\alpha_{i}$ the
$C^{\infty}$ $(1,-1)$-form and $\beta_{i}$ the $C^{\infty}$ function on $S$
described by

$\alpha_{i}(z)=\frac{d^{2}}{dt^{2}}\bigg{|}_{t=0}\frac{\mu_{i}^{t}(z)}{1-|\mu_{i}^{t}(z)|^{2}},\hskip
2.84526pt\beta_{i}(z)=\frac{d^{2}}{dt^{2}}\bigg{|}_{t=0}\frac{|\mu_{i}^{t}(z)|^{2}}{1-|\mu_{i}^{t}(z)|^{2}}.$

By the formula (11),

$\displaystyle\frac{d^{2}}{dt^{2}}\bigg{|}_{t=0}\sum_{i=1}^{n}\mathcal{E}(S_{t},\pi_{i}\circ(f_{i}^{t})^{-1})$
$\displaystyle=\frac{d^{2}}{dt^{2}}\bigg{|}_{t=0}\sum_{i=1}^{n}\Big{(}\mathcal{E}(S_{t},\pi_{i}\circ(f_{i}^{t})^{-1})-\mathcal{E}(S,\pi_{i})\Big{)}$
$\displaystyle=\frac{d^{2}}{dt^{2}}\bigg{|}_{t=0}4\sum_{i=1}^{n}\Big{(}-\textrm{Re}\int_{S}\phi_{i}\cdot\frac{\mu_{i}^{t}}{1-|\mu_{i}^{t}|^{2}}dxdy+\int_{S}|\phi_{i}|\cdot\frac{|\mu_{i}^{t}|^{2}}{1-|\mu_{i}^{t}|^{2}}dxdy\Big{)}.$

Taking the derivative into the integral,

(12)
$0>\frac{d^{2}}{dt^{2}}\bigg{|}_{t=0}\sum_{i=1}^{n}\mathcal{E}(S_{t},\pi_{i}\circ(f_{i}^{t})^{-1})=4\textrm{Re}\sum_{i=1}^{n}\Big{(}-\int_{S}\phi_{i}\cdot\alpha_{i}+\int_{S}|\phi_{i}|\cdot\beta_{i}dA\Big{)}.$

We perturb $\phi_{1},\dots,\phi_{n-1}$ ever so slightly to be generic, and
then redefine

$\phi_{n}:=-\phi_{1}-\dots-\phi_{n-1},$

which is, of course, very close to the original $\phi_{n}$. If $\phi_{n}$ is
not generic, then we perturb it to be. By openness, if our perturbation of
$\phi_{n}$ is sufficiently small, then if we redefine $\phi_{n-1}$ so that the
sum of the $\phi_{i}$’s is again zero, $\phi_{n-1}$ will still be generic.
Since $\alpha_{i}$ and $\beta_{i}$ are uniformly bounded, it is clear from the
right hand side of (12) that if the perturbations are small enough, then both
sides of (12) remain negative. So, we can take these new $\phi_{i}$’s to be
our holomorphic quadratic differentials, and keep the same Riemann surfaces
$S_{t}$ and paths of $C^{\infty}$ maps $t\mapsto f_{i}^{t}:S\to S_{t}.$ ∎

###### Proof of Theorem 1.1.

We find an $\mathbb{R}$-tree such that energy can be decreased to second
order. We will then invoke Proposition 2.4 to say that the same happens for
the extremal length of the associated foliation.

Our starting point is Theorem 3.8: we fix a $\tau$-unstable minimal map from a
branched cover of a Riemann surface $S$ of genus $g$. By Proposition 3.10, we
can adjust the Hopf differentials of the component maps to obtain generic
quadratic differentials $\phi_{1},\dots,\phi_{n}$ that yield an action
$\rho=(\rho_{1},\dots,\rho_{n})$ on a product of $\mathbb{R}$-trees and an
equivariant minimal map $\pi=(\pi_{1},\dots,\pi_{n})$, as well as paths of
$C^{\infty}$ maps $t\mapsto f_{1}^{t},\dots,f_{n}^{t}:S\to S_{t}$ starting at
the identity such that

(13)
$\frac{d^{2}}{dt^{2}}\bigg{|}_{t=0}\sum_{i=1}^{n}\mathcal{E}(S_{t},\pi_{i}\circ(f_{i}^{t})^{-1})<0.$

Since the $\phi_{i}$’s are generic, their energy functionals on
$\mathbf{T}_{g}$ are all real analytic. By the definition of minimality and
(3),

(14) $\frac{d}{dt}|_{t=0}\mathbf{E}_{\rho}(S_{t})=0\textrm{ and
}\frac{d}{dt}\bigg{|}_{t=0}\sum_{i=1}^{n}\mathcal{E}(S_{t},\pi_{i}\circ(f_{i}^{t})^{-1})=0.$

Since harmonic maps minimize energy, for all $t$,

(15)
$\mathbf{E}_{\rho}(S_{t})=\sum_{i=1}^{n}\mathbf{E}_{\rho_{i}}(S_{t})\leq\sum_{i=1}^{n}\mathcal{E}(S_{t},\pi_{i}\circ(f_{i}^{t})^{-1}).$

It follows from (13), (14), and (15) that

(16)
$\frac{d^{2}}{dt^{2}}|_{t=0}\mathbf{E}_{\rho}(S_{t})\leq\frac{d^{2}}{dt^{2}}\bigg{|}_{t=0}\sum_{i=1}^{n}\mathcal{E}(S_{t},\pi_{i}\circ(f_{i}^{t})^{-1})<0.$

By the first equation in (14), the left hand side of (16) does not depend on
the specific path in Teichmüller space, but only on the initial tangent vector
$\mu$. Thus, if we fix a metric $m$ in $\mathcal{C}$, we can replace the path
with the geodesic for $m$ starting at $S$ and tangent to $\mu$ at time zero,
say $t\mapsto S_{t}^{\prime}$. Finally, since the energy splits into the
energies of the component maps, for (16) to hold along our path there must be
at least one component representation $\rho_{i}$ such that

$\frac{d^{2}}{dt^{2}}|_{t=0}\mathbf{E}_{\rho_{i}}(S_{t}^{\prime})<0.$

If $\mathcal{F}$ is the measured foliation corresponding to $\rho_{i},$ then
by Proposition 2.4,

$\frac{d^{2}}{dt^{2}}|_{t=0}\mathbf{EL}_{\mathcal{F}}(S_{t}^{\prime})<0.$

This completes the proof. ∎

######  Remark 3.12.

One could instead ask about convexity with respect to a connection: a function
is convex with respect to a connection if the restriction to every geodesic
for that connection is a convex function. This setting is considered in [1].
Our proof shows that extremal length is not convex for any connection that
admits $C^{2}$ geodesics through every tangent vector of $\mathbf{T}_{g}.$

######  Remark 3.13.

In principle, one can write down the tangent vectors explicitly. One begins
with a $\tau$-unstable minimal surface and a destabilizing variation. For
example, one could work in genus $3$, and take $C=S$ and any non-planar
equivariant minimal map from $\tilde{S}\to\mathbb{R}^{3}$ with its
destabilizing unit normal variation (see [15, Section 5.3]). The proofs of
Propositions 4.4, 4.6, and 5.1 in [15] explain how to build the flows
$f_{1}^{t},f_{2}^{t},f_{3}^{t}.$ One can then try to compute the corresponding
path $t\mapsto S_{t}$ and take the derivative at time zero, although the
computation may be involved.

## References

  * [1] Christophe Bavard.  Théorie de Voronoĭ géométrique. Propriétés de finitude pour les familles de réseaux et analogues.  Bull. Soc. Math. France, 133(2):205–257, 2005.
  * [2] Maxime Fortier Bourque, Dídac Martínez-Granado, and Franco Vargas Pallete.  The extremal length systole of the bolza surface, 2021.
  * [3] Maxime Fortier Bourque and Kasra Rafi.  Non-convex balls in the Teichmüller metric.  J. Differential Geom., 110(3):379–412, 2018.
  * [4] J. Eells and L. Lemaire.  Deformations of metrics and associated harmonic maps.  Proc. Indian Acad. Sci. Math. Sci., 90(1):33–45, 1981.
  * [5] Albert Fathi, François Laudenbach, and Valentin Poénaru.  Thurston’s work on surfaces, volume 48 of Mathematical Notes.  Princeton University Press, Princeton, NJ, 2012.  Translated from the 1979 French original by Djun M. Kim and Dan Margalit.
  * [6] Frederick P. Gardiner and Howard Masur.  Extremal length geometry of Teichmüller space.  Complex Variables Theory Appl., 16(2-3):209–237, 1991.
  * [7] John Hubbard and Howard Masur.  Quadratic differentials and foliations.  Acta Math., 142(3-4):221–274, 1979.
  * [8] Steven P. Kerckhoff.  The asymptotic geometry of Teichmüller space.  Topology, 19(1):23–41, 1980.
  * [9] Nicholas J. Korevaar and Richard M. Schoen.  Sobolev spaces and harmonic maps for metric space targets.  Comm. Anal. Geom., 1(3-4):561–659, 1993.
  * [10] Anna Lenzhen and Kasra Rafi.  Length of a curve is quasi-convex along a Teichmüller geodesic.  J. Differential Geom., 88(2):267–295, 2011.
  * [11] Lixin Liu and Weixu Su.  Variation of extremal length functions on Teichmüller space.  Int. Math. Res. Not. IMRN, (21):6411–6443, 2017.
  * [12] Vladimir Marković.  Uniqueness of minimal diffeomorphisms between surfaces.  Bull. Lond. Math. Soc., 53(4):1196–1204, 2021.
  * [13] Vladimir Marković.  Non-uniqueness of minimal surfaces in a product of closed Riemann surfaces.  Geom. Funct. Anal., 32(1):31–52, 2022.
  * [14] Vladimir Markovic and Nathaniel Sagman.  Minimal surfaces and the new main inequality.  arXiv preprint arXiv:2301.00249, 2022.
  * [15] Vladimir Markovic, Nathaniel Sagman, and Peter Smillie.  Unstable minimal surfaces in $\mathbb{R}^{n}$ and in products of hyperbolic surfaces.  arXiv preprint arXiv:2206.02938, 2022.
  * [16] Howard Masur.  The Teichmüller flow is Hamiltonian.  Proc. Amer. Math. Soc., 123(12):3739–3747, 1995.
  * [17] Hideki Miyachi.  Extremal length functions are log-plurisubharmonic.  In In the tradition of Ahlfors-Bers. VII, volume 696 of Contemp. Math., pages 225–250. Amer. Math. Soc., Providence, RI, 2017.
  * [18] Huiping Pan and Michael Wolf.  Ray structures on teichmüller space, 2022.
  * [19] Mary Rees.  Teichmüller distance for analytically finite surfaces is $C^{2}$.  Proc. London Math. Soc. (3), 85(3):686–716, 2002.
  * [20] Mary Rees.  Teichmüller distance is not $C^{2+\epsilon}$.  Proc. London Math. Soc. (3), 88(1):114–134, 2004.
  * [21] Edgar Reich and Kurt Strebel.  On the Gerstenhaber-Rauch principle.  Israel J. Math., 57(1):89–100, 1987.
  * [22] H. L. Royden.  Automorphisms and isometries of Teichmüller space.  In Advances in the Theory of Riemann Surfaces (Proc. Conf., Stony Brook, N.Y., 1969), Ann. of Math. Studies, No. 66, pages 369–383. Princeton Univ. Press, Princeton, N.J., 1971.
  * [23] Nathaniel Sagman and Peter Smillie.  Unstable minimal surfaces in symmetric spaces of non-compact type.  arXiv preprint arXiv:2208.04885, 2022.
  * [24] William P. Thurston.  Minimal stretch maps between hyperbolic surfaces, 1998.
  * [25] Richard A. Wentworth.  Energy of harmonic maps and Gardiner’s formula.  In In the tradition of Ahlfors-Bers. IV, volume 432 of Contemp. Math., pages 221–229. Amer. Math. Soc., Providence, RI, 2007.
  * [26] Michael Wolf.  Harmonic maps from surfaces to $\mathbb{R}$-trees.  Math. Z., 218(4):577–593, 1995.
  * [27] Michael Wolf.  On realizing measured foliations via quadratic differentials of harmonic maps to $\mathbf{R}$-trees.  J. Anal. Math., 68:107–120, 1996.


# Integrating LSTMs and GNNs for COVID-19 Forecasting

Nathan Sesti  Juan Jose Garau-Luis  Edward Crawley  Bruce Cameron

###### Abstract

The spread of COVID-19 has coincided with the rise of Graph Neural Networks
(GNNs), leading to several studies proposing their use to better forecast the
evolution of the pandemic. Many such models also include Long Short Term
Memory (LSTM) networks, a common tool for time series forecasting. In this
work, we further investigate the integration of these two methods by
implementing GNNs within the gates of an LSTM and exploiting spatial
information. In addition, we introduce a skip connection which proves critical
to jointly capture the spatial and temporal patterns in the data. We validate
our daily COVID-19 new cases forecast model on data of 37 European nations for
the last 472 days and show superior performance compared to state-of-the-art
graph time series models based on mean absolute scaled error (MASE). This area
of research has important applications to policy-making and we analyze its
potential for pandemic resource control.

COVID-19, Graph Neural Networks, LSTM, policy-making

##  1 Introduction

Since its outbreak in late 2019, the COVID-19 virus (Fauci et al., 2020;
Velavan & Meyer, 2020) has devastated the world, causing over 3 million deaths
(World Health Organization, 2021) and economic losses estimated up to $10
trillion (United Nations, 2020). Crucial to curbing such damage is prompt
response and effective policy-making, and, even with vaccination numbers
increasing, this ability will be crucial to develop as the possibility of a
future pandemic looms. Due to the exponential nature of epidemiological
transmission (Li et al., 2020), even slight improvements on early intervention
can have outsized impact, making preemptiveness one of our best tools.
Therefore, knowledge of future spread is of critical importance.

Given the urgency of the pandemic, the Machine Learning (ML) community has
stepped in to help in multiple capacities (Alimadadi et al., 2020; Buckee et
al., 2020), from assisting in early COVID-19 diagnosis based on CT scans
(Barstugan et al., 2020; Wang et al., 2021), to predicting new cases and
hospitalizations using different time series models (Alazab et al., 2020).
Among the latter efforts, some authors have taken advantage of both the
temporal correlations and the inherent network structure of epidemiological
data and thus combined temporal models such as Long Short Term Memory (LSTM)
networks (Hochreiter & Schmidhuber, 1997) with recent advancements in Graph
Neural Networks (GNNs) (Scarselli et al., 2008). This has proven to be one of
the most successful approaches.

In this work we aim to improve upon their accuracy, presenting a new time
series forecasting model based on GNNs and LSTMs. We introduce a general
method that further integrates both networks and apply it to the problem of
forecasting new COVID-19 cases. Our specific use case consists of 37 countries
in Europe, with data from Jan. 2020 to May 2021. We show that our model is
able to forecast new cases out of the training distribution, showing better
results than four other GNN+LSTM models and lag, with a prediction error of
10% when forecasting 7 days ahead. We also analyze how the model performance
translates into real-world policy-making by means of the fraction of missed
cases. Finally, our findings suggest that introducing a skip connection in our
model is key to reducing underfitting.

Figure 1: Overview of the problem and our proposed method. We consider an area
with $N$ different regions (e.g., different countries, states), each
displaying a different COVID-19 case curve. The task consists of predicting,
$M$ timesteps in advance (we use days in our paper), the number of cases in
each region. To that end, a graph with $N$ nodes is created based on the
geographical structure of the regions and a feature vector is associated to
each node. These are fed to our model, consisting of three layers of GraphSAGE
convolutions, an LSTM with graph convolution gates, three additional GraphSAGE
layers, and two final fully-connected layers. We use a skip connection between
the first batch of GraphSAGE layers and the GraphLSTM.

##  2 Related Work

Following the outbreak of COVID-19, the ML community has given the problem of
epidemiological modeling a refreshed solution. The majority of the studies
have made LSTMs the tool of choice. In addition, many researchers have
answered the call for use of mobility data (Buckee et al., 2020; Santosh,
2020) and used graphs to model spatial information, some converging upon
different combinations of LSTMs and GNNs to solve the problem, mirroring their
use in other applications.

Temporal-Only Models Since early in the pandemic, researchers have realized
the value of LSTM networks due to the salient temporal patterns of
transmission (Chimmula & Zhang, 2020; Shahid et al., 2020; Arora et al.,
2020). Their predictions displayed low error, competing well with standard
epidemiological models, such as susceptible-infected-recovered-dead (SIRD)
models (Bailey et al., 1975).

GNNs for COVID-19 Capitalizing on the call for use of mobility data, numerous
papers have applied GNNs to the task. Some have leveraged classic
epidemiological models, one feeding the concatenation of outputs from a GNN
and an epidemiological model into a final distributional regression layer (La
Gatta et al., 2021). Alternatively, (Fritz et al., 2021) uses a combination of
a GNN and LSTM to predict the parameters of a SIRD model. Our work will differ
in that it is not dependent upon any explicit epidemiological model, instead
being able to construct its own parametrization of input features.

While many studies take advantage of both the spatial and temporal
dependencies of the problem, (Cao et al., 2021) pays particular attention to
modeling them jointly — i.e., leveraging the correlations that depend on both
other regions and timesteps rather than taking these separately. While they
achieve this using Fourier Transforms, we instead draw from (Seo et al.,
2018), which replaces the typical linear transformation within an LSTM with a
Chebyshev Spectral CNN (Defferrard et al., 2016), although we switch it out
with GraphSAGE (Hamilton et al., 2017). Our contribution also includes a novel
skip connection, which concatenates the output of the aforementioned LSTM
module with a purely spatial GNN in order to improve upon the challenges
common to joint modeling of spatial and temporal dimensions.

Other Applications of GNN + LSTM At the intersection of these two network
models has been an explosion of research aiming to tackle other graph-
structured time series problems. Specifically, DCRNN (Li et al., 2017) was
applied to traffic prediction, GCLSTM (Chen et al., 2018) was applied to
predicting emails and other contact between people, and GConvGRU (Seo et al.,
2018) was able to predict moving MNIST data and model natural language.

##  3 Problem Statement

Our objective in this paper is to predict some numerical feature of each node
in a graph given a sequence of previous graph snapshots. In the case of
COVID-19, we want to predict new cases across different regions, where each
region constitutes a different node in the graph. The data can be
characterized as a graph $\mathcal{G}=(X,A,W)$ where
$X=\\{x_{ijt}\\}\in\mathbb{R}^{K_{X}\times N\times T}$ indicates the $K_{X}$
input features for each of $N$ nodes at $T$ different timesteps, with
$X_{t}\in\mathbb{R}^{K_{X}\times N}$ indicating the values at timestep $t$.
Then, $A$ is the adjacency matrix —fixed over time— and
$W\in\mathbb{R}^{K_{W}\times N\times N}$ represents each of the $K_{W}$ static
edge features used in our model between each pair of countries.

The task at timestep $t$ is, given the $L$ most recent timesteps
$\\{X_{t-L+1},\mathellipsis,X_{t}\\}$, to predict the number of new cases $M$
timesteps into the future, $\hat{X}_{t+M}\in\mathbb{R}^{N}$. When each graph
snapshot corresponds to a unique day, the task becomes predicting new cases
$M$ days ahead, which benefits effective policy-making and early response if
$M$ is large enough.

##  4 Methods

Our method consists of a combination of GraphSAGE and LSTM layers that jointly
exploit spatio-temporal relationships in the data. This section covers the
main details of the model, wich is depicted in Figure 1. Our approach
considers edge features in the GraphSAGE layers and replaces linear operations
in the LSTM by graph convolutions.

GraphSAGE Our model makes extensive use of GraphSAGE (Hamilton et al., 2017),
a spatial GNN that aggregates the features of local nodes to generate an
embedding for each node in a graph, $\mathcal{G}=(\mathcal{V},\mathcal{E})$.
Let $\boldsymbol{x}_{j}=\boldsymbol{h}_{j}^{0}$ represent the feature vector
of node $j$, which is also the embedding for that node before any iterations.
In each iteration $k\in\\{1,...,n\\}$, each node $v\in\mathcal{V}$ aggregates
the embeddings of all immediately adjacent nodes
$\\{\boldsymbol{h}_{u}^{k-1},\ \forall u\in\mathcal{N}(v)\\}$, usually by
taking the mean of these vectors. This aggregated vector,
$\boldsymbol{h}_{\mathcal{N}(v)}^{k-1}$, is then concatenated to the current
embedding $\boldsymbol{h}_{v}^{k-1}$ and fed into a single neural network
layer with a sigmoid activation function to compute $\boldsymbol{h}_{v}^{k}$.
We slightly alter this scheme to allow for edge features $e_{vu}$ between any
two nodes $u$ and $v$, taking a weighted mean,
$\boldsymbol{h}_{\mathcal{N}(v)}^{k-1}=\frac{1}{|\mathcal{N}(v)|}\sum\limits_{u\in\mathcal{N}(v)}{e_{vu}\cdot
h_{u}^{k-1}}$. In order to use multiple edge features, we learn to calculate
the scalar $e_{vu}$ from a feature vector $\boldsymbol{e}_{vu}$ using a
single-layer perceptron.

GraphLSTM LSTMs (Hochreiter & Schmidhuber, 1997) are a class of recurrent
neural architectures which use a sequence of gates to keep track of both
short- and long-term dependencies in the data. Each of these gates usually
applies weights via matrix multiplication; however, we instead replace each of
these linear transformations with a graph convolution operation, specifically
GraphSAGE, in order to maintain the spatial structure of our inputs and
jointly model spatio-temporal dependencies. We refer to this modification of
the LSTM architecture as GraphLSTM.

Skip Connection The overall structure of our model (see Figure 1) consists of
an initial GraphSAGE layer with $n=3$, feeding into a GraphLSTM with embedded
GraphSAGE, followed by a final GraphSAGE layer also using $n=3$. A multi-layer
perceptron with ReLU activation is applied to finally generate a prediction.
Rather than being fed only the output of the GraphLSTM, the final GraphSAGE
layer recieves a concatenation of the output of the first GraphSAGE layer and
the LSTM. This skip connection is an important contribution of our work, as it
is found to speed stabilization of the model and discourage underfitting.

Refinements In addition to the presented structure, in both GraphSAGE layers,
we make use of common ML techniques including: 1) dropout, which randomly
disables nodes in the graph in order to discourage overfitting and increase
model robustness (Srivastava et al., 2014), and 2) residual connections, which
concatenate raw inputs to posterior layers in order to address the vanishing
gradient problem (He et al., 2016).

Loss Function The specific loss function we use is mean absolute scaled error
(MASE). To train the model we look at the per-person error, i.e., we aggregate
all geographical regions and consider the absolute number of COVID-19 cases.
Specifically,

$\mathcal{L}_{t}=\frac{|\sum_{i=0}^{N-1}{\hat{X}_{ti}-X_{ti}}|}{\sum_{i=0}^{N-1}{X_{ti}}}$
(1)

where $N$ is the total number of geographical regions considered.

Table 1: Per-person and per-country MASE on the test data for our model and
each of the 5 compared approaches, including Lag prediction.
Model | Per-person | Per-country
---|---|---
Our Model | 0.10 | 0.27
GConvLSTM | 0.83 | 2.15
GConvGRU | 0.78 | 1.59
DCRNN | 0.83 | 2.15
GCLSTM | 0.83 | 1.66
Lag | 0.13 | 0.30
Figure 2: Predictions and labels for the test data, consisting of new COVID-19
cases in 37 European countries for the period 22 Mar. 2021 - 9 May 2021 (48
days).

##  5 Results

We now apply our model111A repository with our code can be found at:
https://github.com/jjgarau/GNND to COVID-19 forecasting; we make use of
(Rozemberczki et al., 2021). We start with the dataset and training procedure,
and then discuss the results.

Dataset and Training Motivated by the dynamics of the virus in Europe (Dye et
al., 2020; Saglietto et al., 2020), we select a use case consisting of daily
new COVID-19 cases for 37 European nations ($N=37$) with more than 100k
inhabitants. We use new COVID-19 cases data from (Dong et al., 2020) for the
period 24 Jan. 2020 - 9 May 2021. We use the first 377 days as training data,
the following 47 days as validation data, and the last 48 days as test data.
Our task consists of predicting the smoothed number222Moving average over one
week of new cases one week ahead ($M=7$). We train our model by individually
passing the data from each of 21 consecutive days ($L=21$), beginning 28 days
prior to our prediction target, and use the output of the 21${}^{\text{st}}$
day as our prediction. The edges connecting countries are chosen based on
proximity, each node being assigned to its 3 nearest neighbors, as determined
by geodesic distance between landmass centroids. We use a single edge feature
($K_{W}=1$), consisting of the social connectivity score between two incident
nations according to (Facebook, 2021).

Baseline Comparison The results of running our model on the aforementioned
data are shown in Figure 2. We compare our approach with four state-of-the-art
models which have been shown effective in similar spatial time series
forecasting problems: GConvLSTM (Seo et al., 2018), GConvGRU (Seo et al.,
2018), DCRNN (Li et al., 2017), and GCLSTM (Chen et al., 2018). We find that
ours displays significant improvements for the task of COVID-19 prediction. We
also compare against lag, which uses the last input date as its prediction.
This mechanism has strong results, and though our model learns a similar
lagged behaviour, it yields a 20% reduction in loss over pure lag. Specific
numbers are presented in Table 1 for both our loss function, per-person MASE,
and per-country MASE, which is instead calculated such that each nation
weights equally regardless of its total number of cases.

Application to Policy-making Both of the aforementioned metrics are useful for
policy-making: per-person MASE in cases of international policy-making, such
as that done by the EU; and per-country MASE in the case where countries must
put national policies in place. To better understand the real impact of our
predictions, we also present an illustration of the fraction of missed cases
(i.e., undercounting) in each country (Figure 3) when using our model on the
48 days of test data. Although time series models are usually trained with
symmetric loss functions such as MASE or MSE, their real-world impact might be
asymmetric, as in the case of COVID-19, where underprediction might hamper
preparedness. In our case, the fraction of missed cases ranges from 0 to 0.16,
and stays below 0.1 for the majority of the countries.

Figure 3: Fraction of missed cases on the 48 days of test data.

Skip Connection Ablation Study To validate the efficacy of our proposed skip
connection, we perform an ablation test, finding that the altered model
achieves an inferior result of 0.76 for per-person MASE (1.68 per-country
MASE). We conclude that the skip connection plays a critical role in learning
the spatio-temporal patterns present in the data.

##  6 Conclusion and Future Work

In this paper, we presented a model which built upon existing approaches to
COVID-19 prediction by further integrating LSTMs and GNNs. The outputs of our
model could prove useful to policy-makers attempting to take preemptive action
by giving them improved knowledge of future evolution of the pandemic. Our
model embeds the GraphSAGE graph convolution operator in place of linear
transformations within the gates of an LSTM to create GraphLSTM, a module
capable of jointly capturing spatio-temporal patterns. We also propose a skip
connection to solve common challenges of time series modeling, which is shown
to be an important addition to the model. Our work also presents a solution
useful beyond the application of epidemiological prediction. Numerous problems
similarly rely on a graph structure with node features changing in time.
Future work may also further capitalize on choice of data—while we used only
single edge and node features, there are certainly other relevant factors
potentially including hospital capacity, poverty rates, and age demographics.

## References

  * Alazab et al. (2020) Alazab, M., Awajan, A., Mesleh, A., Abraham, A., Jatana, V., and Alhyari, S.  Covid-19 prediction and detection using deep learning.  _International Journal of Computer Information Systems and Industrial Management Applications_ , 12:168–181, 2020.
  * Alimadadi et al. (2020) Alimadadi, A., Aryal, S., Manandhar, I., Munroe, P. B., Joe, B., and Cheng, X.  Artificial intelligence and machine learning to fight covid-19.  _Physiological Genomics_ , 52(4):200–202, 2020\.
  * Arora et al. (2020) Arora, P., Kumar, H., and Panigrahi, B. K.  Prediction and analysis of covid-19 positive cases using deep learning models: A descriptive case study of india.  _Chaos, Solitons & Fractals_, 139:110017, 2020.  ISSN 0960-0779.
  * Bailey et al. (1975) Bailey, N. T. et al.  _The mathematical theory of infectious diseases and its applications_.  Charles Griffin & Company Ltd, 5a Crendon Street, High Wycombe, Bucks HP13 6LE., 1975.
  * Barstugan et al. (2020) Barstugan, M., Ozkaya, U., and Ozturk, S.  Coronavirus (COVID-19) Classification using CT Images by Machine Learning Methods.  _arXiv preprint arXiv:2003.09424_ , mar 2020.
  * Buckee et al. (2020) Buckee, C. O., Balsari, S., Chan, J., Crosas, M., Dominici, F., Gasser, U., Grad, Y. H., Grenfell, B., Halloran, M. E., Kraemer, M. U. G., Lipsitch, M., Metcalf, C. J. E., Meyers, L. A., Perkins, T. A., Santillana, M., Scarpino, S. V., Viboud, C., Wesolowski, A., and Schroeder, A.  Aggregated mobility data could help fight covid-19.  _Science_ , 368(6487):145–146, 2020.  ISSN 0036-8075.
  * Cao et al. (2021) Cao, D., Wang, Y., Duan, J., Zhang, C., Zhu, X., Huang, C., Tong, Y., Xu, B., Bai, J., Tong, J., et al.  Spectral temporal graph neural network for multivariate time-series forecasting.  _arXiv preprint arXiv:2103.07719_ , 2021.
  * Chen et al. (2018) Chen, J., Xu, X., Wu, Y., and Zheng, H.  Gc-lstm: Graph convolution embedded lstm for dynamic link prediction.  _arXiv preprint arXiv:1812.04206_ , 2018.
  * Chimmula & Zhang (2020) Chimmula, V. K. R. and Zhang, L.  Time series forecasting of covid-19 transmission in canada using lstm networks.  _Chaos, Solitons & Fractals_, 135:109864, 2020.  ISSN 0960-0779.
  * Defferrard et al. (2016) Defferrard, M., Bresson, X., and Vandergheynst, P.  Convolutional neural networks on graphs with fast localized spectral filtering.  _arXiv preprint arXiv:1606.09375_ , 2016.
  * Dong et al. (2020) Dong, E., Du, H., and Gardner, L.  An interactive web-based dashboard to track covid-19 in real time.  _The Lancet infectious diseases_ , 20(5):533–534, 2020.
  * Dye et al. (2020) Dye, C., Cheng, R. C., Dagpunar, J. S., and Williams, B. G.  The scale and dynamics of covid-19 epidemics across europe.  _Royal Society open science_ , 7(11):201726, 2020\.
  * Facebook (2021) Facebook.  Social Connectedness Index, 2021.  URL https://dataforgood.fb.com/tools/social-connectedness-index/.  Accessed: 8 May 2021.
  * Fauci et al. (2020) Fauci, A. S., Lane, H. C., and Redfield, R. R.  Covid-19 — navigating the uncharted.  _New England Journal of Medicine_ , 382(13):1268–1269, 2020.
  * Fritz et al. (2021) Fritz, C., Dorigatti, E., and Rügamer, D.  Combining graph neural networks and spatio-temporal disease models to predict covid-19 cases in germany.  _arXiv preprint arXiv:2101.00661_ , 2021.
  * Hamilton et al. (2017) Hamilton, W. L., Ying, R., and Leskovec, J.  Inductive representation learning on large graphs.  _arXiv preprint arXiv:1706.02216_ , 2017.
  * He et al. (2016) He, K., Zhang, X., Ren, S., and Sun, J.  Deep residual learning for image recognition.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pp. 770–778, 2016.
  * Hochreiter & Schmidhuber (1997) Hochreiter, S. and Schmidhuber, J.  Long short-term memory.  _Neural computation_ , 9(8):1735–1780, 1997.
  * La Gatta et al. (2021) La Gatta, V., Moscato, V., Postiglione, M., and Sperlí, G.  An epidemiological neural network exploiting dynamic graph structured data applied to the covid-19 outbreak.  _IEEE Transactions on Big Data_ , 7(1):45–55, 2021.
  * Li et al. (2020) Li, Q., Guan, X., Wu, P., Wang, X., Zhou, L., Tong, Y., Ren, R., Leung, K. S., Lau, E. H., Wong, J. Y., Xing, X., Xiang, N., Wu, Y., Li, C., Chen, Q., Li, D., Liu, T., Zhao, J., Liu, M., Tu, W., Chen, C., Jin, L., Yang, R., Wang, Q., Zhou, S., Wang, R., Liu, H., Luo, Y., Liu, Y., Shao, G., Li, H., Tao, Z., Yang, Y., Deng, Z., Liu, B., Ma, Z., Zhang, Y., Shi, G., Lam, T. T., Wu, J. T., Gao, G. F., Cowling, B. J., Yang, B., Leung, G. M., and Feng, Z.  Early transmission dynamics in wuhan, china, of novel coronavirus–infected pneumonia.  _New England Journal of Medicine_ , 382(13):1199–1207, 2020.
  * Li et al. (2017) Li, Y., Yu, R., Shahabi, C., and Liu, Y.  Diffusion convolutional recurrent neural network: Data-driven traffic forecasting.  _arXiv preprint arXiv:1707.01926_ , 2017.
  * Rozemberczki et al. (2021) Rozemberczki, B., Scherer, P., He, Y., Panagopoulos, G., Riedel, A., Astefanoaei, M., Kiss, O., Beres, F., Lopez, G., Collignon, N., and Sarkar, R.  PyTorch Geometric Temporal: Spatiotemporal Signal Processing with Neural Machine Learning Models, 2021.
  * Saglietto et al. (2020) Saglietto, A., D’Ascenzo, F., Zoccai, G. B., and De Ferrari, G. M.  Covid-19 in europe: the italian lesson.  _Lancet_ , 395(10230):1110–1111, 2020.
  * Santosh (2020) Santosh, K.  Covid-19 prediction models and unexploited data.  _Journal of medical systems_ , 44(9):1–4, 2020\.
  * Scarselli et al. (2008) Scarselli, F., Gori, M., Tsoi, A. C., Hagenbuchner, M., and Monfardini, G.  The graph neural network model.  _IEEE transactions on neural networks_ , 20(1):61–80, 2008.
  * Seo et al. (2018) Seo, Y., Defferrard, M., Vandergheynst, P., and Bresson, X.  Structured sequence modeling with graph convolutional recurrent networks.  In _International Conference on Neural Information Processing_ , pp. 362–373. Springer, 2018.
  * Shahid et al. (2020) Shahid, F., Zameer, A., and Muneeb, M.  Predictions for covid-19 with deep learning models of lstm, gru and bi-lstm.  _Chaos, Solitons & Fractals_, 140:110212, 2020.  ISSN 0960-0779.
  * Srivastava et al. (2014) Srivastava, N., Hinton, G., Krizhevsky, A., Sutskever, I., and Salakhutdinov, R.  Dropout: a simple way to prevent neural networks from overfitting.  _The journal of machine learning research_ , 15(1):1929–1958, 2014.
  * United Nations (2020) United Nations.  COVID-19 to slash global economic output by $8.5 trillion over next two years, 2020.  URL https://www.un.org/development/desa/en/news/policy/wesp-mid-2020-report.html.  Accessed: 8 May 2021.
  * Velavan & Meyer (2020) Velavan, T. P. and Meyer, C. G.  The covid-19 epidemic.  _Tropical medicine & international health_, 25(3):278, 2020.
  * Wang et al. (2021) Wang, S.-H., Govindaraj, V. V., Górriz, J. M., Zhang, X., and Zhang, Y.-D.  Covid-19 classification by fgcnet with deep feature fusion from graph convolutional network and convolutional neural network.  _Information Fusion_ , 67:208–229, 2021.  ISSN 1566-2535.
  * World Health Organization (2021) World Health Organization.  WHO Coronavirus (COVID-19) Dashboard, 2021.  URL https://covid19.who.int/.  Accessed: 8 May 2021.