EUROPEAN ORGANIZATION FOR NUCLEAR RESEARCH

Search for a feebly interacting particle $X$
in the decay $K^{+}\rightarrow\pi^{+}X$
The NA62 Collaboration

###### Abstract

A search for the $K^{+}\rightarrow\pi^{+}X$ decay, where $X$ is a long-lived
feebly interacting particle, is performed through an interpretation of the
$K^{+}\rightarrow\pi^{+}\nu\bar{\nu}$ analysis of data collected in 2017 by
the NA62 experiment at CERN. Two ranges of $X$ masses,
$0$–$110\,\text{MeV}/c^{2}$ and $154$–$260\,\text{MeV}/c^{2}$, and lifetimes
above $100\,\text{ps}$ are considered. The limits set on the branching ratio,
$\text{BR}(K^{+}\rightarrow\pi^{+}X)$, are competitive with previously
reported searches in the first mass range, and improve on current limits in
the second mass range by more than an order of magnitude.

Accepted for publication by JHEP

The NA62 Collaboration 111Corresponding authors: J. Swallow, R. Volpe, email:
<EMAIL_ADDRESS><EMAIL_ADDRESS>

Université Catholique de Louvain, Louvain-La-Neuve, Belgium
E. Cortina Gil, A. Kleimenova, E. Minucci 111Corresponding authors: J.
Swallow, R. Volpe, email<EMAIL_ADDRESS>roberta.volpe@cern.ch${}^{,}\,$222Deceased, S. Padolski 333Present address:
Brookhaven National Laboratory, Upton, NY 11973, USA, P. Petrov, A. Shaikhiev
444Also at Institute for Nuclear Research of the Russian Academy of Sciences,
117312 Moscow, Russia, R. Volpe$^{,}$111Corresponding authors: J. Swallow, R.
Volpe, email<EMAIL_ADDRESS><EMAIL_ADDRESS>555Present address: Faculty of Mathematics, Physics and Informatics, Comenius
University, 842 48, Bratislava, Slovakia

TRIUMF, Vancouver, British Columbia, Canada
T. Numao, B. Velghe

University of British Columbia, Vancouver, British Columbia, Canada
D. Bryman 666Also at TRIUMF, Vancouver, British Columbia, V6T 2A3, Canada, J.
Fu 777Present address: UCLA Physics and Biology in Medicine, Los Angeles, CA
90095, USA

Charles University, Prague, Czech Republic
T. Husek 888Present address: Department of Astronomy and Theoretical Physics,
Lund University, Lund, SE 223-62, Sweden, J. Jerhot 999Present address:
Université Catholique de Louvain, B-1348 Louvain-La-Neuve, Belgium, K. Kampf,
M. Zamkovsky

Institut für Physik and PRISMA Cluster of excellence, Universität Mainz,
Mainz, Germany
R. Aliberti 101010Present address: Institut für Kernphysik and Helmholtz
Institute Mainz, Universität Mainz, Mainz, D-55099, Germany, G. Khoriauli
111111Present address: Universität Würzburg, D-97070 Würzburg, Germany, J.
Kunze, D. Lomidze 121212Present address: European XFEL GmbH, D-22761 Hamburg,
Germany, R. Marchevski 131313Present address: CERN, European Organization for
Nuclear Research, CH-1211 Geneva 23, Switzerland, L. Peruzzo, M. Vormstein, R.
Wanke

Dipartimento di Fisica e Scienze della Terra dell’Università e INFN, Sezione
di Ferrara, Ferrara, Italy
P. Dalpiaz, M. Fiorini, I. Neri, A. Norton, F. Petrucci, H. Wahl

INFN, Sezione di Ferrara, Ferrara, Italy
A. Cotta Ramusino, A. Gianoli

Dipartimento di Fisica e Astronomia dell’Università e INFN, Sezione di
Firenze, Sesto Fiorentino, Italy
E. Iacopini, G. Latino, M. Lenti, A. Parenti

INFN, Sezione di Firenze, Sesto Fiorentino, Italy
A. Bizzeti 141414Also at Dipartimento di Fisica, Università di Modena e Reggio
Emilia, I-41125 Modena, Italy, F. Bucci

Laboratori Nazionali di Frascati, Frascati, Italy
A. Antonelli, G. Georgiev 151515Also at Faculty of Physics, University of
Sofia, BG-1164 Sofia, Bulgaria, V. Kozhuharov 151515Also at Faculty of
Physics, University of Sofia, BG-1164 Sofia, Bulgaria, G. Lanfranchi, S.
Martellotti, M. Moulson, T. Spadaro

Dipartimento di Fisica “Ettore Pancini” e INFN, Sezione di Napoli, Napoli,
Italy
F. Ambrosino, T. Capussela, M. Corvino 131313Present address: CERN, European
Organization for Nuclear Research, CH-1211 Geneva 23, Switzerland, D. Di
Filippo, P. Massarotti, M. Mirra, M. Napolitano, G. Saracino

Dipartimento di Fisica e Geologia dell’Università e INFN, Sezione di Perugia,
Perugia, Italy
G. Anzivino, F. Brizioli, E. Imbergamo, R. Lollini, R. Piandani 161616Present
address: Institut für Experimentelle Teilchenphysik (KIT), D-76131 Karlsruhe,
Germany, C. Santoni

INFN, Sezione di Perugia, Perugia, Italy
M. Barbanera 171717Present address: School of Physics and Astronomy,
University of Birmingham, Birmingham, B15 2TT, UK, P. Cenci, B. Checcucci, P.
Lubrano, M. Lupi 181818Present address: Institut am Fachbereich Informatik und
Mathematik, Goethe Universität, D-60323 Frankfurt am Main, Germany, M. Pepe,
M. Piccini

Dipartimento di Fisica dell’Università e INFN, Sezione di Pisa, Pisa, Italy
F. Costantini, L. Di Lella, N. Doble, M. Giorgi, S. Giudici, G. Lamanna, E.
Lari, E. Pedreschi, M. Sozzi

INFN, Sezione di Pisa, Pisa, Italy
C. Cerri, R. Fantechi, L. Pontisso, F. Spinella

Scuola Normale Superiore e INFN, Sezione di Pisa, Pisa, Italy
I. Mannelli

Dipartimento di Fisica, Sapienza Università di Roma e INFN, Sezione di Roma I,
Roma, Italy
G. D’Agostini, M. Raggi

INFN, Sezione di Roma I, Roma, Italy
A. Biagioni, E. Leonardi, A. Lonardo, P. Valente, P. Vicini

INFN, Sezione di Roma Tor Vergata, Roma, Italy
R. Ammendola, V. Bonaiuto 191919Also at Department of Industrial Engineering,
University of Roma Tor Vergata, I-00173 Roma, Italy, A. Fucci, A. Salamon, F.
Sargeni 202020Also at Department of Electronic Engineering, University of Roma
Tor Vergata, I-00173 Roma, Italy

Dipartimento di Fisica dell’Università e INFN, Sezione di Torino, Torino,
Italy
R. Arcidiacono 212121Also at Università degli Studi del Piemonte Orientale,
I-13100 Vercelli, Italy, B. Bloch-Devaux, M. Boretto 131313Present address:
CERN, European Organization for Nuclear Research, CH-1211 Geneva 23,
Switzerland, E. Menichetti, E. Migliore, D. Soldi

INFN, Sezione di Torino, Torino, Italy
C. Biino, A. Filippi, F. Marchetto

Instituto de Física, Universidad Autónoma de San Luis Potosí, San Luis Potosí,
Mexico
J. Engelfried, N. Estrada-Tristan 222222Also at Universidad de Guanajuato,
Guanajuato, Mexico

Horia Hulubei National Institute of Physics for R&D in Physics and Nuclear
Engineering, Bucharest-Magurele, Romania
A. M. Bragadireanu, S. A. Ghinescu, O. E. Hutanu

Joint Institute for Nuclear Research, Dubna, Russia
A. Baeva, D. Baigarashev, D. Emelyanov, T. Enik, V. Falaleev, V. Kekelidze, A.
Korotkova, L. Litov 151515Also at Faculty of Physics, University of Sofia,
BG-1164 Sofia, Bulgaria, D. Madigozhin, M. Misheva 232323Present address:
Institute of Nuclear Research and Nuclear Energy of Bulgarian Academy of
Science (INRNE-BAS), BG-1784 Sofia, Bulgaria, N. Molokanova, S. Movchan, I.
Polenkevich, Yu. Potrebenikov, S. Shkarovskiy, A. Zinchenko 222Deceased

Institute for Nuclear Research of the Russian Academy of Sciences, Moscow,
Russia
S. Fedotov, E. Gushchin, A. Khotyantsev, Y. Kudenko 242424Also at National
Research Nuclear University (MEPhI), 115409 Moscow and Moscow Institute of
Physics and Technology, 141701 Moscow region, Moscow, Russia, V. Kurochka, M.
Medvedeva, A. Mefodev

Institute for High Energy Physics - State Research Center of Russian
Federation, Protvino, Russia
S. Kholodenko, V. Kurshetsov, V. Obraztsov, A. Ostankov 222Deceased, V.
Semenov 222Deceased, V. Sugonyaev, O. Yushchenko

Faculty of Mathematics, Physics and Informatics, Comenius University,
Bratislava, Slovakia
L. Bician 131313Present address: CERN, European Organization for Nuclear
Research, CH-1211 Geneva 23, Switzerland, T. Blazek, V. Cerny, Z. Kucerova

CERN, European Organization for Nuclear Research, Geneva, Switzerland
J. Bernhard, A. Ceccucci, H. Danielsson, N. De Simone 252525Present address:
DESY, D-15738 Zeuthen, Germany, F. Duval, B. Döbrich, L. Federici, E.
Gamberini, L. Gatignon, R. Guida, F. Hahn 222Deceased, E. B. Holzer, B.
Jenninger, M. Koval 262626Present address: Charles University, 116 36 Prague
1, Czech Republic, P. Laycock 333Present address: Brookhaven National
Laboratory, Upton, NY 11973, USA, G. Lehmann Miotto, P. Lichard, A. Mapelli,
K. Massri, M. Noy, V. Palladino 272727Present address: Physics Department,
Imperial College London, London, SW7 2BW, UK, M. Perrin-Terrin 282828Present
address: Aix Marseille University, CNRS/IN2P3, CPPM, F-13288, Marseille,
France${}^{,}\,$292929Also at Université Catholique de Louvain, B-1348
Louvain-La-Neuve, Belgium, J. Pinzino 303030Present address: Department of
Physics, University of Toronto, Toronto, Ontario, M5S 1A7,
Canada${}^{,}\,$313131Also at INFN, Sezione di Pisa, I-56100 Pisa, Italy, V.
Ryjov, S. Schuchmann 323232Present address: Institut für Physik and PRISMA
Cluster of excellence, Universität Mainz, D-55099 Mainz, Germany, S. Venditti

University of Birmingham, Birmingham, United Kingdom
T. Bache, M. B. Brunetti 333333Present address: Department of Physics,
University of Warwick, Coventry, CV4 7AL, UK, V. Duk 343434Present address:
INFN, Sezione di Perugia, I-06100 Perugia, Italy, V. Fascianelli 353535Present
address: Dipartimento di Psicologia, Università di Roma La Sapienza, I-00185
Roma, Italy, J. R. Fry, F. Gonnella, E. Goudzovski, L. Iacobuzio, C.
Lazzeroni, N. Lurkin 999Present address: Université Catholique de Louvain,
B-1348 Louvain-La-Neuve, Belgium, F. Newson, C. Parkinson 999Present address:
Université Catholique de Louvain, B-1348 Louvain-La-Neuve, Belgium, A. Romano,
A. Sergi, A. Sturgess, J. Swallow 111Corresponding authors: J. Swallow, R.
Volpe, email<EMAIL_ADDRESS><EMAIL_ADDRESS>

University of Bristol, Bristol, United Kingdom
H. Heath, R. Page, S. Trilov

University of Glasgow, Glasgow, United Kingdom
B. Angelucci, D. Britton, C. Graham, D. Protopopescu

University of Lancaster, Lancaster, United Kingdom
J. Carmignani, J. B. Dainton, R. W. L. Jones, G. Ruggiero

University of Liverpool, Liverpool, United Kingdom
L. Fulton, D. Hutchcroft, E. Maurice 363636Present address: Laboratoire
Leprince Ringuet, F-91120 Palaiseau, France, B. Wrona

George Mason University, Fairfax, Virginia, USA
A. Conovaloff, P. Cooper, D. Coward 373737Also at SLAC National Accelerator
Laboratory, Stanford University, Menlo Park, CA 94025, USA, P. Rubin

11footnotetext: Present address: Laboratori Nazionali di Frascati, I-00044
Frascati, Italy22footnotetext: Also at CERN, European Organization for Nuclear
Research, CH-1211 Geneva 23, Switzerland

##  1 Introduction

Some scenarios Beyond the Standard Model of particle physics (BSM) include a
new light feebly interacting particle $X$, which can be produced in
$K^{+}\rightarrow\pi^{+}X$ decays. In a hidden sector portal framework the new
$X$ particle mediates interactions between standard model (SM) and hidden
sector fields [1]. In the Higgs portal scenario, $X$ is a scalar that mixes
with the SM Higgs boson; this is realised in inflationary [2], scale invariant
[3], and relaxion [4] models, which additionally have cosmological
implications. A massless $X$ particle would have the properties of a neutral
boson arising from the spontaneous breaking of a global $U(1)$ symmetry [5]:
$X$ may then acquire mass through explicit symmetry breaking. One example,
arising from the breaking of a Peccei-Quinn (PQ) symmetry, is an axion [6, 7],
which would be a signature of the PQ mechanism and credibly solve the strong
CP problem [8, 9]. Such an axion could be flavor non-diagonal [10].
Alternatives, from breaking of the lepton number and flavour symmetries
respectively, are majorons [11] or familons [7, 12]. A QCD axion with mass
$\mathcal{O}(10^{-4}\,\text{eV})$ could be a dark matter candidate, and
specific axion models can also solve the SM flavor problem [13]. In a broader
class of models, $X$ is considered as an axion-like particle (ALP) that acts
as a pseudoscalar mediator [14]. Alternatively the introduction of a light,
feebly-coupled, spin-1 boson can effectively generate through its axial
couplings the phenomenology related to an invisible spin-0 ALP [15].

Searches for $X$ production in the $K^{+}\rightarrow\pi^{+}X$ decay have the
potential to constrain many BSM models. The $K^{+}\rightarrow\pi^{+}X$ decay
is characterised by an incoming $K^{+}$, an outgoing $\pi^{+}$ and missing
energy-momentum, as is the rare $K^{+}\rightarrow\pi^{+}\nu\bar{\nu}$ decay.
An interpretation of the NA62 $K^{+}\rightarrow\pi^{+}\nu\bar{\nu}$ studies
using 2017 data [16] in terms of a search for the $K^{+}\rightarrow\pi^{+}X$
decay is presented here. Upper limits are established on
$\text{BR}(K^{+}\rightarrow\pi^{+}X)$ and interpreted in terms of two BSM
scenarios.

##  2 Beamline, detector and dataset

The NA62 experiment, beamline and detector are described in detail in [17] and
a schematic of the detector is shown in Fig. 1. A right-handed coordinate
system, $(x,y,z)$, is defined with the target at the origin and the beam
travelling towards positive $z$, the $y$ axis is vertical (positive up) and
the $x$-axis is horizontal (positive left). A $400\,\text{GeV/}c$ proton beam
extracted from the CERN Super Proton Synchrotron (SPS) impinges on a beryllium
target creating a $75\,\text{GeV/}c$ secondary hadron beam with a $1\%$ rms
momentum spread and a composition of $70\%$ pions, $23\%$ protons and $6\%$
kaons. Kaons ($K^{+}$) are positively tagged with $70\,\text{ps}$ timing
precision by the KTAG detector, a differential Cherenkov counter filled with
nitrogen gas. The momentum and position of the $K^{+}$ are measured by the
GigaTracker (GTK), a spectrometer formed of three silicon pixel tracker
stations and a set of four dipole magnets. GTK measurements have momentum,
direction and time resolutions of $0.15\,\text{GeV}/c$, $16\,\mu\text{rad}$
and $100\,\text{ps}$, respectively. After traversing the GTK magnets, a
magnetized scraper used to sweep away muons, and a bending magnet (B), the
beam at the FV entrance has a rectangular profile of $52\times
24\,\text{mm}^{2}$ and a divergence of $0.11\,\text{mrad}$.

Fig. 1: Schematic top view of the NA62 beamline and detector. The “CHOD” label
indicates both the CHOD and NA48-CHOD hodoscopes described in the text. Also
shown is the trajectory of a beam particle in vacuum which crosses all the
detector apertures, thus avoiding interactions with material. A dipole magnet
between MUV3 and SAC deflects the beam particles out of the SAC acceptance.

The experiment is designed to study $K^{+}$ decays occurring in the
$60\,\text{m}$ fiducial volume (FV) starting $2.6\,\text{m}$ downstream of
GTK3 and housed inside a $117\,\text{m}$ long vacuum tank, containing a
magnetic spectrometer, and ending at the ring imaging Cherenkov counter
(RICH). Momentum and position measurements for charged particles produced in
$K^{+}$ decays in the FV are provided by the magnetic spectrometer composed of
four STRAW tracking stations, two on either side of a dipole magnet (M). This
spectrometer provides a momentum measurement with resolution $\sigma_{p}/p$ of
$0.3$–$0.4\%$. The RICH is filled with neon gas at atmospheric pressure and
provides particle identification for charged particles, and a time measurement
with a precision better than $100\,\text{ps}$. Two adjacent scintillator
hodoscopes (CHOD and NA48-CHOD), provide time measurements for charged
particles with a $200\,\text{ps}$ resolution.

A system of veto detectors is key to the experiment. Interactions of beam
particles in GTK3 are detected by the charged particle anti-counter (CHANTI),
formed of six stations of scintillator bar counters. Downstream, a photon veto
system is used to reject the $K^{+}\rightarrow\pi^{+}\pi^{0}$ background. This
analysis selects $\pi^{+}$ particles with momenta in the range
$15$–$35\,\text{GeV/}c$. This means that a $\pi^{0}$ from the
$K^{+}\rightarrow\pi^{+}\pi^{0}$ background has momentum of at least
$40\,\text{GeV/}c$ and the subsequent $\pi^{0}\rightarrow\gamma\gamma$ decay,
$\text{BR}=98.8\%$, produces two energetic photons which can be detected with
high efficiency. There are twelve large angle veto (LAV) stations positioned
to ensure hermetic coverage for photon emission angles of
$8.5$–$50\,\text{mrad}$. The liquid krypton calorimeter (LKr) provides
coverage for $1$–$8.5\,\text{mrad}$. The small angle photon veto (SAV) covers
angles below $1\,\text{mrad}$ using two sampling calorimeters of shashlyk
design (IRC and SAC).

Downstream of the LKr are two hadronic sampling calorimeters (MUV1 and MUV2).
Together with the LKr, these provide particle identification information
through the pattern of energy deposition. Electrons/positrons produce
electromagnetic showers that are well-contained in the LKr, which has a depth
of 27 radiation lengths. Pions may pass through the LKr without losing all of
their energy and can produce a hadronic shower in MUV1 and MUV2. In contrast,
muons are minimum ionising particles in the calorimetric system. The MUV3
detector is positioned downstream of a $0.8\,\text{m}$ iron absorber and
consists of a plane of scintillator tiles. It provides measurements of muons
with $400\,\text{ps}$ time resolution.

A two-level trigger system is employed with a hardware level 0 (L0) selection
followed by a level 1 (L1) decision made by software algorithms. The primary
trigger stream of the experiment is dedicated to collection of
$K^{+}\rightarrow\pi^{+}\nu\bar{\nu}$ events and uses information from the
CHOD, RICH, LKr, MUV3 at L0 [18] and KTAG, LAV, STRAW at L1 [16]. The
NA48-CHOD also provides a $99\%$ efficient minimum bias trigger, used for
collection of $K^{+}\rightarrow\pi^{+}\pi^{0}$ events that are used for
normalisation. The data sample collected in 2017 for the study of the
$K^{+}\rightarrow\pi^{+}\nu\bar{\nu}$ decay is used for this analysis.

##  3 Signal selection

The observable for the $K^{+}\rightarrow\pi^{+}X$ search is the reconstructed
squared missing mass

$m^{2}_{\mathrm{miss}}=(P_{K}-P_{\pi})^{2}\,,$

where $P_{K}$ and $P_{\pi}$ are the $K^{+}$ and $\pi^{+}$ 4-momenta, derived
from the measured 3-momenta of the GTK and STRAW tracks under the $K^{+}$ and
$\pi^{+}$ mass hypotheses, respectively. The event selection is identical to
that used for the $K^{+}\rightarrow\pi^{+}\nu\bar{\nu}$ measurement [16] and
is summarised below.

Candidate events must have fewer than three reconstructed STRAW tracks with no
negatively charged tracks. Only one track can fulfil additional criteria to
become a $\pi^{+}$ candidate but, for example, an additional out-of-time halo
muon track may exist. The time assigned to the $\pi^{+}$ candidate is
calculated using the mean times measured in the STRAW, NA48-CHOD and RICH
weighted by their respective measured resolutions. A $\pi^{+}$ candidate track
must have momentum in the range $15$–$35\,\text{GeV}/c$ and be within the
sensitive regions of the downstream detectors (RICH, CHODs, LKr and MUV1,2,3)
with geometrically and time-coincident associated signals recorded in the
CHODs, LKr and RICH.

The candidate track must be consistent with the $\pi^{+}$ hypothesis for the
RICH reconstructed mass and likelihood. The candidate must also satisfy a
multivariate classifier based on calorimetric information. On average, for
$15$–$35\,\text{GeV}/c$ tracks, the two methods achieve $\pi^{+}$
identification efficiencies of $82\%$ and $78\%$, with probabilities of
misidentification of $\mu^{+}$ as $\pi^{+}$ of $2.3\times 10^{-3}$ and
$6.3\times 10^{-6}$, respectively. A MUV3 veto condition rejects events with
signals geometrically associated with the track within a time window of
$7\,\text{ns}$. No signals are allowed in any LAV station (or SAV) within
$3\,(7)\,\text{ns}$ of the $\pi^{+}$ time. No LKr clusters are allowed beyond
a distance of $100\,\text{mm}$ from the $\pi^{+}$ impact point within cluster-
energy dependent time windows of $10$ to $100\,\text{ns}$. The STRAW, CHODs
and LKr are used to veto events with additional activity, including tracks
produced by photon interactions upstream of the calorimeters and partially
reconstructed multi-track decays. Overall rejection of
$\pi^{0}\rightarrow\gamma\gamma$ decays is achieved with an inefficiency of
$1.3\times 10^{-8}$.

A $K^{+}$ is tagged upstream by the KTAG if Cherenkov photons are detected
within $2\,\text{ns}$ of the $\pi^{+}$ track time in at least five out of its
total of eight sectors. A GTK track is associated with the $K^{+}$ if its time
is within $0.6\,\text{ns}$ of the KTAG time and the closest distance of
approach (CDA) to the $\pi^{+}$ track is less than $4\,\text{mm}$. The
$K^{+}$/$\pi^{+}$ matching is based on time coincidence and spatial
information and has an efficiency of $75\%$. The average probability for wrong
(accidental) association with pileup GTK tracks is $1.3\%$ ($3.5\%$) when the
$K^{+}$ track is (is not) correctly reconstructed.

Upstream backgrounds arise from a combination of early $K^{+}$ decays
(upstream of the FV), beam particle interactions in the GTK stations,
additional GTK tracks, and large-angle $\pi^{+}$ scattering in the first STRAW
station. To minimise such backgrounds, the vertex formed between the selected
$K^{+}$ and $\pi^{+}$ tracks must be inside the FV with no additional activity
in the CHANTI within $3\,\text{ns}$ of the $\pi^{+}$ candidate time.
Additionally, a ‘box cut’ is applied requiring that the projection of the
$\pi^{+}$ candidate track back to the final collimator (COL) is outside the
area defined by $|x|<100\,\text{mm}$ and $|y|<500\,\text{mm}$.

The $m^{2}_{\mathrm{miss}}$ observable is used to discriminate between a
peaking two-body $K^{+}\rightarrow\pi^{+}X$ signal and backgrounds. Two signal
regions are defined, called region 1 and region 2, to minimise large
backgrounds from $K^{+}\rightarrow\pi^{+}\pi^{0}$,
$K^{+}\rightarrow\mu^{+}\nu_{\mu}$ and $K^{+}\rightarrow\pi^{+}\pi^{+}\pi^{-}$
decays. The reconstructed $m^{2}_{\mathrm{miss}}$ for region 1 must be between
$0$ and $0.01\,\text{GeV}^{2}/c^{4}$ and that for region 2 between $0.026$ and
$0.068\,\text{GeV}^{2}/c^{4}$. Additional momentum-dependent constraints
supplement the definition of the signal regions using alternative squared
missing mass variables, constructed either by replacing the GTK measurement of
the beam 3-momentum with the average beam momentum and direction, or the STRAW
3-momentum measurement with one measured by the RICH under the $\pi^{+}$ mass
hypothesis. These requirements reject events with incorrect reconstruction of
$m^{2}_{\mathrm{miss}}$ due to momenta mismeasurements and improve background
rejection, but decrease acceptance at the boundaries of the signal regions.

##  4 Signal and background models

Geant4-based [19] Monte Carlo simulations of $K^{+}\rightarrow\pi^{+}X$ decays
are performed with the assumption that $X$ is stable, for $X$ masses covering
the search range at $1.4\,\,\text{MeV}/c^{2}$ intervals. This value
corresponds to intervals of the squared missing mass that are always smaller
than its resolution. These simulations include decay kinematics, interactions
in material, and the responses of the detectors. In this study, a scan is
performed searching for $K^{+}\rightarrow\pi^{+}X$ signals with $X$ mass,
$m_{X}$, in the ranges $0$–$110\,\text{MeV/}c^{2}$ and
$154$–$260\,\text{MeV/}c^{2}$. These $m_{X}$ ranges extend beyond the
$K^{+}\rightarrow\pi^{+}\nu\bar{\nu}$ signal regions because of the resolution
of the reconstructed $m^{2}_{\mathrm{miss}}$ observable. The resolution of
$m^{2}_{\mathrm{miss}}$, $\sigma_{m^{2}_{\text{miss}}}$, as a function of
simulated $m_{X}$ is shown in Fig. 2 (left). The reconstructed
$m^{2}_{\mathrm{miss}}$ resolution for a control sample of selected
$K^{+}\rightarrow\pi^{+}\pi^{0}$ events is found to be $4\%$ better in
simulations than in data. The resolution derived from simulations is therefore
corrected by increasing it by $4\%$ and a systematic uncertainty of $10\%$ is
assigned to the $m^{2}_{\mathrm{miss}}$ resolution. The acceptance for the
selection described in section 3, as obtained using simulations, is displayed
in Fig. 2 (centre). The single event sensitivity, $\text{BR}_{SES}$, defined
as the branching ratio corresponding to the observation of one signal event,
is calculated by following the procedure adopted for the
$K^{+}\rightarrow\pi^{+}\nu\bar{\nu}$ analysis using the
$K^{+}\rightarrow\pi^{+}\pi^{0}$ decay for normalisation [16]; the resulting
values are shown in Fig. 2 (right). The uncertainty of $BR_{SES}$ is $10\%$
and is mainly systematic. The largest contributions to this uncertainty are
associated with the trigger efficiency, signal and normalisation
reconstruction and selection efficiencies [16], and differences between
$K^{+}\rightarrow\pi^{+}\nu\bar{\nu}$ and $K^{+}\rightarrow\pi^{+}X$
kinematics.

Fig. 2:  Resolution of the $m^{2}_{\mathrm{miss}}$ observable (left),
acceptance (centre) and single event sensitivity, $\text{BR}_{SES}$, (right)
for $K^{+}\rightarrow\pi^{+}X$, as functions of mass hypothesis $m_{X}$.

The sensitivity for low $X$ masses is limited by the
$K^{+}\rightarrow\pi^{+}\nu\bar{\nu}$ signal region definition
$m^{2}_{\mathrm{miss}}>0$, which is necessary to suppress the background from
$K^{+}\rightarrow\mu^{+}\nu_{\mu}$ decays. This effect reduces the acceptance
by half for $m_{X}=0$, and equivalently at each signal region boundary (Fig. 2
centre).

The acceptance for $X$ with finite lifetime, $\tau_{X}$ and $m_{X}\neq 0$, is
computed under the following assumptions: $X$ decays only to visible SM
particles; decays upstream of MUV3 are detected with $100\%$ efficiency. The
efficiency is $99.9\%$, and the uncertainty in this quantity is included in
the systematic uncertainty. The acceptance for a set of $\tau_{X}$ values is
calculated by weighting simulated events by the probability that $X$ does not
decay upstream of MUV3. The acceptance increases as a function of lifetime
reaching a plateau for $\tau_{X}>10\,\text{ns}$. For
$m_{X}<20\,\text{MeV/}c^{2}$, losses of acceptance at lower lifetimes are
compensated by the increase in the Lorentz factor.

The background contributions for the $K^{+}\rightarrow\pi^{+}X$ search are the
same as for the $K^{+}\rightarrow\pi^{+}\nu\bar{\nu}$ analysis with the
addition of the $K^{+}\rightarrow\pi^{+}\nu\bar{\nu}$ decay itself, which
becomes the dominant background. The SM description of the
$K^{+}\rightarrow\pi^{+}\nu\bar{\nu}$ decay is assumed. The total expected
background and the reconstructed $m^{2}_{\mathrm{miss}}$ distributions for
each component are obtained from auxiliary measurements, as described in [16].
The resulting numbers of background events in the signal regions are
summarised in Table 1. The contributions from kaon decays other than
$K^{+}\rightarrow\pi^{+}\nu\bar{\nu}$ are grouped in the row other $K^{+}$
decays, and their distribution in $m^{2}_{\mathrm{miss}}$ is known with good
accuracy. For the upstream background, an additional systematic uncertainty of
$30\%$ is included, to account for the uncertainty in the estimation of its
distribution in $m^{2}_{\mathrm{miss}}$ resulting from the limited size of the
control sample used for the auxiliary measurements. The total background is
described, as a function of the reconstructed $m^{2}_{\mathrm{miss}}$, by
fitting polynomial functions to the expectations in signal regions 1 and 2, as
shown in Fig. 3.

Fig. 3:  Distributions of the expected reconstructed squared missing mass,
$m^{2}_{\mathrm{miss}}$, for background processes, obtained from simulations
and data-driven procedures, displayed as stacked histograms with bin width
$0.00067\,\text{GeV}^{2}/c^{4}$. In each signal region, the polynomial
function used to describe the total background is shown. Table 1: Summary of
the predicted numbers of background events in the signal regions and the
observed events. The statistical uncertainty for SM
$K^{+}\rightarrow\pi^{+}\nu\bar{\nu}$ is negligible and the external
uncertainty arises from the uncertainty of the SM
$K^{+}\rightarrow\pi^{+}\nu\bar{\nu}$ branching ratio.

Region 1 Region 2 $K^{+}\rightarrow\pi^{+}\nu\bar{\nu}$ (SM) $0.55$
$\pm\,0.07_{syst}\pm 0.13_{ext}$ $1.61$ $\pm\,0.11_{syst}\pm 0.22_{ext}$
Upstream background $0.21\pm 0.12_{stat}\pm 0.12_{syst}$ $0.68\pm
0.21_{stat}\pm 0.26_{syst}$ Other $K^{+}$ decays $0.26\pm 0.04_{stat}\pm
0.05_{syst}$ $0.31\pm 0.04_{stat}\pm 0.06_{syst}$ Total background $1.02\pm
0.13_{stat}\pm 0.15_{syst}\pm 0.13_{ext}$ $2.60\pm 0.21_{stat}\pm
0.28_{syst}\pm 0.22_{ext}$ Observed events $0$ $2$

##  5 Statistical analysis

The search procedure involves a fully frequentist hypothesis-test using a
shape analysis with observable $m^{2}_{\mathrm{miss}}$ and an unbinned profile
likelihood ratio test statistic. Each $X$ mass hypothesis is treated
independently. The parameter of interest,
$\text{BR}(K^{+}\rightarrow\pi^{+}X)$, is related to the expected number of
signal events, $n_{S}$, by
$\text{BR}(K^{+}\rightarrow\pi^{+}X)=n_{S}\times\text{BR}_{SES}$.

The likelihood function has the form:

$\mathcal{L}=\frac{(n_{tot})^{n}e^{-n_{tot}}}{n!}\times\prod_{j}^{n}\Bigg{[}\frac{n_{B}}{n_{tot}}f_{B}\Big{(}m^{2}_{\mathrm{miss}},_{j}\Big{)}+\frac{n_{S}}{n_{tot}}f_{S}\Big{(}m^{2}_{\mathrm{miss}},_{j}|\mu_{X},\sigma_{X}\Big{)}\Bigg{]}\times\prod^{N_{\text{nuis}}}_{i}\mathcal{C}_{i}(p_{\text{meas}}^{i}|p_{\text{nuis}}^{i})\,$

where $n$ is the observed number of events, $n_{tot}=n_{B}+n_{S}$ and $n_{B}$
is the expected number of background events; $f_{B}(m^{2}_{\mathrm{miss}})$ is
a polynomial function of $m^{2}_{\mathrm{miss}}$ normalised to unity which
describes the total background in the signal region relevant for a certain
mass hypothesis $m_{X}$; and $f_{S}(m^{2}_{\mathrm{miss}}|\mu_{X},\sigma_{X})$
is the Gaussian function, normalised to unity, with parameters $\mu_{X}$ and
$\sigma_{X}$ obtained from a fit to the distribution of the reconstructed
simulated events. Index $j$ runs over the $n$ observed events and their
reconstructed $m^{2}_{\mathrm{miss}}$ are denoted
$m^{2}_{\mathrm{miss}},_{j}$. The $N_{\text{nuis}}$ nuisance parameters
considered, $p_{\text{nuis}}^{i}$, are
$n_{B},\,\text{BR}_{SES},\,\mu_{X},\,\sigma_{X}$, and are estimated by
auxiliary measurements. These estimations,
$p_{\text{meas}}^{i}=\hat{n}_{B},\,\hat{\text{BR}}_{SES},\,\hat{\mu}_{X},\,\hat{\sigma}_{X}$,
are treated as global observables [20]. The constraint terms,
$\mathcal{C}_{i}(p_{\text{meas}}^{i}|p_{\text{nuis}}^{i})$, are the
probability density functions describing the distribution of each nuisance
parameter. The constraint term for $n_{B}$ is a Poisson distribution with mean
value $(\hat{n}_{B}/\sigma_{B})^{2}$ where $\hat{n}_{B}$ and $\sigma_{B}$ are
the central value and uncertainty of the background expectation [21]. The
constraint term for $\text{BR}_{SES}$ is a log-normal function with parameters
corresponding to a relative uncertainty of $10\%$. A Gaussian constraint term
is used for $\mu_{X}$, with relative uncertainty depending on the mass
hypothesis $m_{X}$. A log-normal constraint term is used for $\sigma_{X}$,
with the mean corresponding to the estimated value after the $4\%$ correction
(described in Section 4), and relative uncertainty of 10%. The normalised
polynomial functions, describing the background distribution in
$m^{2}_{\mathrm{miss}}$, are considered to be known exactly.

For each mass hypothesis the fully frequentist test is performed according to
the CLs method [22] to exclude the presence of a signal with $90\%$ confidence
level (CL) for the observed data. A cross-check was performed, using single
bin counting experiments in windows of width equal to four times
$\sigma_{m^{2}_{\text{m}iss}}$ around each mass hypothesis, with a hybrid
frequentist treatment using a log-likelihood ratio test statistic. A
comparable expected sensitivity was obtained.

##  6 Results and discussion

Two candidate $K^{+}\rightarrow\pi^{+}X$ events are observed [16] at
reconstructed $m_{\text{miss}}$ values of $196$ and $252\,\text{MeV}/c^{2}$.
Upper limits are established on $\text{BR}(K^{+}\rightarrow\pi^{+}X)$ at
$90\%$ CL for each $X$ mass hypothesis: expected and observed upper limits,
assuming stable or invisibly decaying $X$, are displayed in Fig. 4 (left). The
observed upper limits are compared to the previous results from the E787/E949
experiments [23] in Fig. 4 (right), as a function of $m_{X}$ and for different
values of $\tau_{X}$, assuming $X$ decays to visible SM particles. The
strongest limits of $5\times 10^{-11}$ are obtained at large $X$ masses
($160$–$250\,\text{MeV}/c^{2}$) and long $X$ lifetimes ($>5\,\text{ns}$).
Under the assumption of stable or invisibly decaying $X$ these upper limits
improve by a factor of $\mathcal{O}(10)$ in signal region 2, and are
competitive in region 1. For unstable $X$, assuming decays only to visible SM
particles, the same pattern holds in general. However, in region 1 the limits
obtained improve across an increasingly large range of mass hypotheses as the
assumed lifetime becomes shorter. Despite differences in experimental set-up
between E787/E949 (stopped $K^{+}$ decay-at-rest) and NA62 (highly boosted
$K^{+}$ decay-in-flight), the two results exhibit similar dependence on
$\tau_{X}$. This is because the ratios of the Lorentz factor for the $X$
particle to the decay length are similar in the two experiments.

Fig. 4:  Left: upper limits on $\text{BR}(K^{+}\rightarrow\pi^{+}X)$ for each
mass hypothesis, $m_{X}$, tested. Right: model-independent observed upper
limits as functions of the mass and lifetime assumed for $X$ for NA62 (this
work, solid lines) and E787/E949 [23] (dashed lines).

In a Higgs portal model with a dark sector scalar mixing with the Higgs boson,
$X$ production and decay are driven by the mixing parameter $\sin^{2}\theta$
(model BC4 [1, 24]). This gives rise to $K^{+}\rightarrow\pi^{+}X$ decays with
branching ratio proportional to $\sin^{2}\theta$. The constraints derived on
$\sin^{2}\theta$ from this search, alongside results from other studies, are
shown in Fig. 5.

In a scenario where $X$ is an ALP with couplings proportional to SM Yukawa
couplings (model BC10 [1, 14]) the $K^{+}\rightarrow\pi^{+}X$ decay occurs
with a branching ratio proportional to the square of the coupling constant
$g_{Y}$. The constraints on $g_{Y}$ derived from this and other searches are
shown in Fig. 6.

If $X$ decays only to invisible particles, such as dark matter, bounds on the
coupling parameter ($\sin^{2}\theta$ or $g_{Y}$ for the scalar and ALP models,
respectively) are directly derived from its relationship with the branching
ratio, with results shown in the right-hand panels of Figs. 5 and 6. If $X$
decays only to visible SM particles, $\tau_{X}$ is inversely proportional to
the coupling parameters [14, 24], limiting the reach of this analysis for
large coupling because of lower acceptance for shorter lifetimes. The
$X\rightarrow e^{+}e^{-}$ decays dominate the visible decay width up to the
di-muon threshold beyond which an additional channel opens and $\tau_{X}$
decreases, limiting the sensitivity of this search. The model-dependent
relationship between the lifetime and coupling therefore determines the shape
of the exclusion regions shown in the left-hand panels of Figs. 5 and 6.

Fig. 5:  Excluded regions of the parameter space $(m_{S},\sin^{2}\theta)$ for
a dark scalar, $S$, of the BC4 model [1] decaying only (left) to visible SM
particles as in the BC4 model and (right) invisibly. The exclusion bound from
the present search for the decay $K^{+}\rightarrow\pi^{+}S$ is labelled as
“$K^{+}\rightarrow\pi^{+}+\text{inv.}$” and is shaded in red. In the $\pi^{0}$
mass region the independent NA62 search for
$\pi^{0}\rightarrow\text{invisible}$ decays [25] provides constraints, shown
in purple. Other bounds, shown in grey, are derived from the experiments E949
[23], CHARM [24], NA48/2 [26], LHCb [27, 28] and Belle [29].

Fig. 6:  Excluded regions of the parameter space $(m_{a},g_{Y})$ for an ALP,
$a$, of the BC10 model [1] decaying only (left) to visible particles and
(right) invisibly. The exclusion bound from the present search for the decay
$K^{+}\rightarrow\pi^{+}a$ is labelled as
“$K^{+}\rightarrow\pi^{+}+\text{inv.}$” and is shaded in red. In the $\pi^{0}$
mass region the independent NA62 search for
$\pi^{0}\rightarrow\text{invisible}$ decays [25] provides constraints, shown
in purple. Other bounds, shown in grey, are derived from the experiments E949
[23], $K_{\mu 2}$ [30], CLEO [31], CHARM [32], KTeV [33], LHCb [27, 28] and
from Big Bang nucleosynthesis (BBN) [1].

##  7 Conclusions

A search for the $K^{+}\rightarrow\pi^{+}X$ decay, where $X$ is a long-lived
feebly interacting particle, is performed through an interpretation of the
$K^{+}\rightarrow\pi^{+}\nu\bar{\nu}$ analysis of data collected in 2017 by
the NA62 experiment at CERN. Two candidate $K^{+}\rightarrow\pi^{+}X$ events
are observed, in agreement with the expected background. Upper limits on
$\text{BR}(K^{+}\rightarrow\pi^{+}X)$ are established at $90\%$ CL, with the
strongest limits of $5\times 10^{-11}$ at large $X$ masses
($160$–$250\,\text{MeV}/c^{2}$) and long $X$ lifetimes ($>5\,\text{ns}$),
improving on current results by up to a factor of $\mathcal{O}(10)$. An
interpretation of these results to constrain BSM models is presented in
scenarios where $X$ is a dark scalar mixing with the Higgs boson or is an ALP
with couplings to fermions.

## Acknowledgements

It is a pleasure to express our appreciation to the staff of the CERN
laboratory and the technical staff of the participating laboratories and
universities for their efforts in the operation of the experiment and data
processing.

The cost of the experiment and its auxiliary systems was supported by the
funding agencies of the Collaboration Institutes. We are particularly indebted
to: F.R.S.-FNRS (Fonds de la Recherche Scientifique - FNRS), Belgium; BMES
(Ministry of Education, Youth and Science), Bulgaria; NSERC (Natural Sciences
and Engineering Research Council), funding SAPPJ-2018-0017 Canada; NRC
(National Research Council) contribution to TRIUMF, Canada; MEYS (Ministry of
Education, Youth and Sports), Czech Republic; BMBF (Bundesministerium für
Bildung und Forschung) contracts 05H12UM5, 05H15UMCNA and 05H18UMCNA, Germany;
INFN (Istituto Nazionale di Fisica Nucleare), Italy; MIUR (Ministero
dell’Istruzione, dell’Università e della Ricerca), Italy; CONACyT (Consejo
Nacional de Ciencia y Tecnología), Mexico; IFA (Institute of Atomic Physics)
Romanian CERN-RO No.1/16.03.2016 and Nucleus Programme PN 19 06 01 04,
Romania; INR-RAS (Institute for Nuclear Research of the Russian Academy of
Sciences), Moscow, Russia; JINR (Joint Institute for Nuclear Research), Dubna,
Russia; NRC (National Research Center) “Kurchatov Institute” and MESRF
(Ministry of Education and Science of the Russian Federation), Russia; MESRS
(Ministry of Education, Science, Research and Sport), Slovakia; CERN (European
Organization for Nuclear Research), Switzerland; STFC (Science and Technology
Facilities Council), United Kingdom; NSF (National Science Foundation) Award
Numbers 1506088 and 1806430, U.S.A.; ERC (European Research Council)
“UniversaLepto” advanced grant 268062, “KaonLepton” starting grant 336581,
Europe.

Individuals have received support from: Charles University Research Center
(UNCE/SCI/ 013), Czech Republic; Ministry of Education, Universities and
Research (MIUR “Futuro in ricerca 2012” grant RBFR12JF2Z, Project GAP), Italy;
Russian Foundation for Basic Research (RFBR grants 18-32-00072, 18-32-00245),
Russia; Russian Science Foundation (RSF 19-72-10096), Russia; the Royal
Society (grants UF100308, UF0758946), United Kingdom; STFC (Rutherford
fellowships ST/J00412X/1, ST/M005798/1), United Kingdom; ERC (grants 268062,
336581 and starting grant 802836 “AxScale”); EU Horizon 2020 (Marie
Skłodowska-Curie grants 701386, 842407, 893101).

## References

  * [1] J. Beacham et al., J. Phys. G 47 (2020) 010501.
  * [2] F. Bezrukov, D. Gorbunov, JHEP 05 (2010) 010.
  * [3] J.D. Clarke, R. Foot, R.R. Volkas, JHEP 02 (2014) 123.
  * [4] A. Banerjee et al., JHEP 07 (2020) 153.
  * [5] M. Tanabashi et al., Phys. Rev. D 98 (2018) 030001.
  * [6] S. Weinberg, Phys. Rev. Lett. 40 (1978) 223.
  * [7] F. Wilczek, Phys. Rev. Lett. 49 (1982) 1549.
  * [8] R.D. Peccei, H.R. Quinn, Phys. Rev. Lett. 38 (1977) 1440.
  * [9] R.D. Peccei, H.R. Quinn, Phys. Rev. D 16 (1977) 1791.
  * [10] M. Hindmarsh, P. Moulatsiotis, Phys. Rev. D 59 (1999) 055015.
  * [11] G.B. Gelmini, M. Roncadelli, Phys. Lett. B 99 (1981) 411.
  * [12] A. Davidson, K.C. Wali, Phys. Rev. Lett. 48 (1982) 11.
  * [13] L. Calibbi et al., Phys. Rev. D 95 (2017) 095009.
  * [14] J.M. Dolan et al., JHEP 03 (2015) 171 [Erratum: JHEP 07 (2015) 103].
  * [15] P. Fayet, Phys. Rev. D 74 (2006) 054034.
  * [16] E. Cortina Gil et al., JHEP 11 (2020) 42.
  * [17] E. Cortina Gil et al., JINST 12 (2017) P05025.
  * [18] R. Ammendola et al., Instrum. Meth. A 929 (2019) 1.
  * [19] J. Allison et al., Nucl. Instrum. Meth. A 835 (2016) 186.
  * [20] K. Cranmer, arXiv:1503.07622.
  * [21] R. Cousins et al., Nucl. Instrum. Meth. A 595 (2008) 480.
  * [22] A.L. Read, J. Phys. G 28 (2002) 2693.
  * [23] A.V. Artamonov et al., Phys. Rev. D 79 (2009) 092004.
  * [24] M.W. Winkler, Phys. Rev. D 99 (2019) 015018.
  * [25] E. Cortina Gil et al., arXiv:2010.07644, submitted to JHEP (2020).
  * [26] J.R. Batley et al., Phys. Lett. B 769 (2017) 67.
  * [27] R. Aaij et al., Phys. Rev. D 95 (2017) 071101.
  * [28] R. Aaij et al. Phys. Rev. Lett. 115 (2015) 161802.
  * [29] J.T. Wei et al., Phys. Rev. Lett. 103 (2009) 171801.
  * [30] T. Yamazaki et al., Phys. Rev. Lett. 52 (1984) 1089.
  * [31] R. Ammar et al., Phys. Rev. Lett. 87 (2001) 271801.
  * [32] F. Bergsma et al., Phys. Lett. B 157 (1985) 458.
  * [33] A. Alavi-Harati et al., Phys. Rev. Lett. 84 (2000) 5279.


#  siMLr: Machine Learning inside the SIR model for COVID-19 Forecasting

Roberto Vega
University of Alberta
<EMAIL_ADDRESS>
&Leonardo Flores
Independent researcher
<EMAIL_ADDRESS>
&Russell Greiner
University of Alberta
<EMAIL_ADDRESS>

###### Abstract

Accurate forecasts of the number of newly infected people during an epidemic
are critical for making effective timely decisions. This paper addresses this
challenge using the siMLr model, which incorporates machine learning (ML) into
the epidemiological SIR model. For each region, siMLr tracks the changes in
the policies implemented at the government level, which it uses to estimate
the time-varying parameters of an SIR model for forecasting the number of new
infections 1- to 4-weeks in advance. It also forecasts the probability of
changes in those government policies at each of these future times, which is
essential for the longer-range forecasts. We applied siMLr to data from
regions in Canada and in the United States, and show that its MAPE (mean
average percentage error) performance is as good as SOTA forecasting models,
with the added advantage of being an interpretable model. We expect that this
approach will be useful not only for forecasting COVID-19 infections, but also
in predicting the evolution of other infectious diseases.

##  1 Introduction

Since its identification in December 2019, COVID-19 has posed critical
challenges for the public health and economies of essentially every country in
the world [6; 7; 15]. Government officials have taken a wide range of measures
in an effort to contain this pandemic, including closing schools and
workplaces, setting restrictions on air travel, and establishing stay at home
requirements [8]. Accurately forecasting the number of new infected people in
the short and medium term is critical for the timely decisions about policies
and for the proper allocation of medical resources [2; 14]. For the general
population (e.g., business owners, parents with children in school age),
additionally forecasting when the government will take action might help them
to plan accordingly.

There are three basic approaches for predicting the dynamics of an epidemic:
compartmental models, statistical methods, and ML-based methods [2; 22].
Compartmental models subdivide a population into mutually exclusive categories
(compartments), with a set of dynamical equations that explain the transitions
among categories [3]. The _Susceptible-Infected-Removed_ (SIR) model [12] is a
common choice for the modelling of infectious diseases. Statistical methods,
as the name implies, extract general statistics from the data to fit
mathematical models that explain the evolution of the epidemic [14]. Finally,
ML-based methods use machine learning algorithms to analyze historical data
and find patterns that lead to accurate predictions of the number of new
infected people [18; 22].

Arguably, when any approach is used to make high-stake decisions, it is
important that it be not just accurate, but also interpretable: It should give
the decision-maker enough information to justify the recommendation [20].
Here, we propose siMLr, which is an interpretable probabilistic graphical
model (PGM) that combines compartmental models and ML-based methods. As its
name suggests, it incorporates machine learning (ML) within an SIR model.

siMLr uses a _mixture of experts_ approach [9], where the contribution of each
_expert_ to the final forecast depends on the changes in the government
policies implemented at various earlier time points. When there is no recent
change in policies (2 to 4 weeks before the week to be predicted), siMLr
relies on an SIR model with time-varying parameters that are fitted using
machine learning methods. When a change in policy occurs, siMLr instead relies
on a simpler model that predicts that the new number of infections will remain
constant. Note that forecasting the number of new infections 1 and 2 weeks in
advance ($\Delta I_{1}$ and $\Delta I_{2}$) is relatively easy as siMLr knows,
at the time of the prediction, whether the policy has changed recently.
However, for 3- or 4-week forecasts ($\Delta I_{3}$ and $\Delta I_{4}$), our
model needs to estimate the likelihood of a future change of policy. siMLr
incorporates prior domain knowledge to estimate such policy-change
probabilities.

This work makes three important contributions. (1) It empirically shows that
an SIR model with time-varying parameters can describe the complex dynamics of
COVID-19. (2) It describes an interpretable model that predicts the new number
of infections 1 to 4 weeks in advance, achieving state-of-the-art results, in
terms of mean average percentage error (MAPE), on data from regions from
Canada and the United States. (3) It effectively estimates the probability of
observing a change in the government policies, 1 to 4 weeks in advance.

The rest of Section 1 describes the related work and basics of the SIR
compartmental model. Section 2 then describes in detail our proposed siMLr
approach. Finally, Section 3 shows the results of the predicting the number of
new infections in the United States and provinces of Canada.

###  1.1 Basic SIR model

The _Susceptible-Infected-Removed_ (SIR) compartmental model [12] is a
mathematical model of infectious disease dynamics that divide the population
into 3 disjoint groups [3]: _Susceptible (S)_ refers to the set of people who
have never been infected, but can acquire the disease if they interact with
someone with the infection. _Infected (I)_ refers to the set of people who
have and can transmit the infection at a given point in time. _Removed (R)_
refers to the people who have either recovered or died from the infection, and
so cannot transmit the disease anymore. The transitions between the groups is
given by the following differential equations:

$\frac{dS}{dt}\ =\ -\beta\ \frac{S(t)I(t)}{N}\qquad\quad\frac{dI}{dt}\ =\
\beta\ \frac{S(t)I(t)}{N}-\gamma I(t)\qquad\quad\frac{dR}{dt}\ =\ \gamma\
I(t)$ (1)

This simple model, which assumes homogeneous and constant population, is fully
defined by the parameters $\beta$ (transmission rate) and $\gamma$ (recovery
rate). The intuition behind this model is that every infected patient gets in
contact with $\beta$ people. Since only the susceptible people can become
infected, the chance of interacting with a susceptible person is simply the
proportion of susceptible people in the entire population, $N=S+I+R$.
Likewise, at every time point, $\gamma$ proportion of the infected people is
removed from the system. Figure 1(a) depicts the general behaviour of an SIR
model.

Figure 1:  (a) General behaviour of the SIR model. (b) The number of
infections predicted by the SIR model with fixed parameters, fitted to the US
data for 1 week prediction. (c) Similar to (b), but with time-varying
parameters.

###  1.2 Related work

The main idea behind combining compartmental models with machine learning is
to replace the fixed parameters of the former with time-varying parameters
that can be learned from data [1; 2; 4; 14; 16; 22]. However, it is still not
clear if these approaches are able to explain (and forecast) the different
waves of COVID-19 infections.

Particularly relevant to our approach is the work by Arik et al. [2], who used
latent variables and autoencoders to model extra compartments in an extended
_Susceptible-Exposed-Infected-Removed_ (SEIR) model [3]. Those additional
compartments bring further insight into how the disease impacts the
population; however, we argue that they are not needed for an accurate
prediction of the number of new infections; see results in Section 3. One
limitation of their model is a decrease in performance when the trend in the
number of new infections changes [2]. We hypothesize that those changes in
trend are related to the government policies that are in place at a specific
point in time. siMLr is able to capture those changes by tracking the policies
implemented at the government level, allowing it to also predict when those
policies are likely to change.

A different line of work uses machine learning methods to directly predict the
number of new infections, without incorporating epidemiological models [10;
11; 17; 19]. Like our approach, Yeung et al. added the non-pharmaceutical
interventions (policies) as features in their models; then they used
traditional machine learning models to make their predictions under the
assumption that the data is not a time-series, but a collection of points that
are independent and identically distributed [23]. Their approach is limited to
make predictions up to 2 weeks in advance, since information about the
policies that will be implemented in the future is not available at inference
time. Our siMLr approach differs by being interpretable and also forecasting
policy changes, which allows it to extend the horizon of the $\Delta I$
predictions.

##  2 Our siMLr model

We view siMLr as a probabilistic graphical model that uses a mixture of
experts approach to forecast the number of new COVID-19 infections, 1 to 4
weeks in advance; Figure 2 depicts this graphical model as a plate model [13].
The blue nodes are estimated at every time point, while the values of the
green nodes are either known as part of the historical data, or inferred in a
previous time point.

The random variables in Figure 2 are assumed to have the following
distributions:

$\begin{array}[]{r@{\quad|\quad}c@{\quad\sim\quad}l}\mbox{CT}_{t+1}&\mbox{CP}_{t-3},\
\mbox{CP}_{t-2},\
\mbox{CP}_{t-1}&Categorical_{K\in\\{-1,0,1\\}}(\theta_{CT})\\\
\beta_{t+1}&\beta_{t-2},\ \beta_{t-1},\ \beta_{t},\
\mbox{CT}_{t+1}&\mathcal{N}(\mu_{\beta},\Sigma_{\beta})\\\
\gamma_{t+1}&\gamma_{t-2},\ \gamma_{t-1},\ \gamma_{t},\
\mbox{CT}_{t+1}&\mathcal{N}(\mu_{\gamma},\Sigma_{\gamma})\\\
\mbox{SIR}_{t+1}&\beta_{t+1},\
\gamma_{t+1}&\mathcal{N}(\mu_{SIR},\Sigma_{SIR})\\\
U_{t}&\mbox{SIR}_{t-2},\mbox{SIR}_{t-1},\
\mbox{SIR}_{t}&Categorical_{K\in\\{-1,0,1\\}}(\theta_{U})\\\
O_{t}&W_{t}&Categorical_{K\in\\{0,1\\}}(\theta_{O})\\\ \mbox{CP}_{t+1}&O_{t},\
U_{t}&Categorical_{K\in\\{-1,0,1\\}}(\theta_{CP})\end{array}$ (2)

where $t$ indexes the current week, $SIR_{t}=[S_{t},I_{t},R_{t}]$,
$\mu_{SIR}\in\mathbb{R}^{3}$ is given below by Equation
3$,\mu_{\beta}=(\alpha_{0,CT_{t+1}})+(\alpha_{1,CT_{t+1}})\beta_{t-1}+(\alpha_{2,CT_{t+1}})\beta_{t-2}+(\alpha_{3,CT_{t+1}})\beta_{t-3}$
and
$\mu_{\gamma}=(\omega_{0,CT_{t+1}})+(\omega_{1,CT_{t+1}})\gamma_{t-1}+(\omega_{2,CT_{t+1}})\gamma_{t-2}+(\omega_{3,CT_{t+1}})\gamma_{t-3}$
are linear combinations of the three previous values of $\beta$ and $\gamma$
(respectively). The coefficients of those linear combinations depend on the
value of the random variable $\mbox{CT}_{t+1}$. We did not specify a
distribution for the node $\mbox{\tt New infections}_{t+1}$ because its value
is deterministically computed as $S_{t}-S_{t+1}$.

Informally, the assignment $\mbox{CT}_{t}=-1$ means that we expect a change in
trend from an increasing number of infections to a decreasing one. The
opposite happens when $\mbox{CT}_{t}=1$, while $\mbox{CT}_{t}=0$ means that we
expect the population to follow the current trend (either increasing or
decreasing). These changes in trend are assumed to depend on changes in the
government policies 2 to 4 weeks prior to the week of our forecast – e.g.,  we
use $\\{\mbox{CT}_{t-3},\,\mbox{CT}_{t-2},\,\mbox{CT}_{t-1}\\}$ when
predicting the number of new infections at $t+1$, $\Delta I_{t+1}$; and we
need $\\{\mbox{CT}_{t},\,\mbox{CT}_{t+1},\,\mbox{CT}_{t+2}\\}$ when predicting
$\Delta I_{t+4}$, Note that, at time $t$, we will not know $\mbox{CT}_{t+1}$
nor $\mbox{CT}_{t+2}$.

The status of $\mbox{CT}_{t+1}$ defines the coefficients that relate
$\beta_{t+1}$ and $\gamma_{t+1}$ with their three previous values
$\beta_{t},\beta_{t-1},\beta_{t-2}$ and
$\gamma_{t},\gamma_{t-1},\gamma_{t-2}$, respectively. Since $\beta_{t+1}$ and
$\gamma_{t+1}$ fully parameterize the SIR model in Equation 1, we can estimate
the new number of infected people, $\Delta I_{t+1}$, from these parameters (as
well as the SIR values at time $t$).

The random variables $U_{t}\in\\{-1,0,1\\}$ and $O_{t}\in\\{0,1\\}$ are
auxiliary variables designed to predict the probability of observing a change
in policy at time $t+1$. Intuitively, $U_{t}$ represents the “urgency” of
modifying a policy. As the number of cases per 100K inhabitants and the rate
of change between the number of cases in two consecutive time points
increases, the urgency to set stricter government policies increases. As the
number (and rate of change) of cases decreases, the urgency to relax the
policies increases. Finally, $O_{t}$ models the “willingness” to execute a
change in government policies. As the number of time points without a change
increases, so does this “willingness”.

Figure 2:  Modeling siMLr as a probabilistic graphical model for forecasting
new cases of COVID-19. The blue nodes are estimated at each time point, while
the green ones are either based on past information, or where estimated in a
previous iteration.

###  2.1 SIR with time-varying parameters

We can approximate an SIR model by transforming the differential Equation 1
into an equation of differences and then solving for $S_{t},I_{t}$ and
$R_{t}$:

$S_{t}\ =\ -\beta\frac{S_{t-1}I_{t-1}}{N}+S_{t-1}\qquad I_{t}\ =\
\beta\frac{S_{t-1}I_{t-1}}{N}-\gamma I_{t-1}+I_{t-1}\qquad R_{t}\ =\ \gamma
I_{t-1}+R_{t-1}$ (3)

While the SIR model is non-linear with respect to the states (S, I, R), it is
linear with respect to the parameters $\beta$ and $\gamma$. Therefore, under
the assumption of constant and known population size (i.e.,
$N=S_{t}+I_{t}+R_{t}$) we can re-write the set of Equations 3 as:

$\displaystyle\begin{bmatrix}S_{t}\\\ I_{t}\end{bmatrix}$
$\displaystyle=\quad\begin{bmatrix}-\frac{S_{t-1}I_{t-1}}{N}&0\\\
\frac{S_{t-1}I_{t-1}}{N}&-I_{t-1}\\\ \end{bmatrix}\begin{bmatrix}\beta\\\
\gamma\end{bmatrix}\quad+\quad\begin{bmatrix}S_{t-1}\\\ I_{t-1}\end{bmatrix}$
(4) $\displaystyle R_{t}$ $\displaystyle=\quad N-S_{t}-I_{t}$

Given a sequence of states $x_{1},\dots,x_{n}$, where $x_{t}\ =\ [S_{t}\
I_{t}]^{T}$, it is possible to estimate the optimal parameters of the SIR
model as:

$(\beta^{*},\
\gamma^{*})\quad=\quad\operatorname*{arg\,min}_{\beta,\gamma}\sum_{i=1}^{n}||x_{i}-\hat{x}_{i}||^{2}+\lambda_{1}(\beta-\beta_{0})^{2}+\lambda_{2}(\gamma-\gamma_{0})^{2}$
(5)

where $\hat{x}_{i}$ is computed using Equation 4, and $\lambda_{1}$ and
$\lambda_{2}$ are optional regularization parameters that allow the
incorporation of the priors $\beta_{0}$ and $\gamma_{0}$. For the case of
Gaussian priors – i.e., $\beta\sim\mathcal{N}(\beta_{0},\sigma_{\beta}^{2})$
and $\gamma\sim\mathcal{N}(\gamma_{0},\sigma_{\gamma}^{2})$ – we use
$\lambda_{1}=\frac{1}{2\sigma_{\beta}^{2}}$ and
$\lambda_{2}=\frac{1}{2\sigma_{\gamma}^{2}}$ [18].

In the traditional SIR model, we set $\lambda_{1}=\lambda_{2}=0$ and fit a
single $\beta$ and $\gamma$ to the entire time series. However, as shown in
Figure 1(a), an SIR model with fixed parameters is unable to accurately model
several waves of infections. As illustration, Figure 1(b) shows the
predictions produced by fitting an SIR with fixed parameters (Equation 5) to
the US data from 29/March/2020 to 3/May/2021, and then using those parameters
to make predictions one week in advance, over this same interval. That is,
using this learned $(\beta,\gamma)$, and the number of people in the $S$, $I$,
and $R$ compartments on 28/March/2020, we predicted the number of observed
cases during the week of 29/March/2020 to 4/April/2020. Then, using the
reported $[S,I,R]$ values on 4/April/2020, we predicted the number of new
infections in the week of 5/April/2020 to 11/April/2020 and so on. Note that
even though the parameters $\beta$ and $\gamma$ were found using the entire
time series – i.e., using information that was not available at the time of
prediction – the resulting model still does a poor job fitting the reported
data.

Figure 1(c), on the other hand, was created by allowing $\beta$ and $\gamma$
to change every week. Here, we first found the parameters that fit the data
from 29/March/2020 to 4/April/2020 – call them $\beta_{1}$ and $\gamma_{1}$ –
then used those parameters along with the SIR state on 28/March/2020 to
predict the number of new infections one week ahead. By repeating this
procedure during the entire time series we obtained an almost perfect fit to
the data. Of course, these are also not “legal” predictions since they too use
information that is not available at prediction time – i.e., they used the
number of reported infections during this first week to find the parameters,
which were then used to estimate the number of cases over this time. However,
this “cheating” example shows that an SIR model, with the best time varying
parameters, can model the complex dynamics of COVID-19, including the
different waves of infections. Recall from Figure 1(b) that this is not the
case in the SIR model with fixed parameters: even in the best scenario
(knowing the optimal parameters based on the true values over months of the
same data you want to fit!), it still cannot properly fit the data.

###  2.2 Estimating $\beta_{t+1}$ and $\gamma_{t+1}$

Naturally, the challenge is “legally” computing the appropriate values of
$\beta_{t+1}$ and $\gamma_{t+1}$, for each week, using only the data that is
known at time $t$. Figure 2 shows that computing $\beta_{t+1}$ and
$\gamma_{t+1}$ depends on the status of the random variable $\mbox{CT}_{t+1}$.
When $\mbox{CT}_{t+1}=0$ – i.e., there is no change in the current trend – we
assume that:

$\begin{gathered}\beta_{t+1}\quad\sim\quad\mathcal{N}(\alpha_{0}+\alpha_{1}\beta_{t}+\alpha_{2}\beta_{t-1}+\alpha_{3}\beta_{t-2},\quad\sigma_{\beta}^{2})\\\
\gamma_{t+1}\quad\sim\quad\mathcal{N}(\omega_{0}+\omega_{1}\gamma_{t}+\omega_{2}\gamma_{t-1}+\omega_{3}\gamma_{t-2},\quad\sigma_{\gamma}^{2})\end{gathered}$
(6)

At time _t_ , we can use the historical _daily_ data $x_{1},x_{2},\dots,x_{t}$
to find the _weekly_ parameters $\beta_{1},\beta_{2},\dots,\beta_{t/7}$ and
$\gamma_{1},\gamma_{2},\dots,\gamma_{t/7}$. The first weekly pair
$(\beta_{1},\gamma_{1})$ is found by fitting Equation 5 to
$x_{1},\dots,x_{7}$; $(\beta_{2},\gamma_{2})$ to $x_{8},\dots,x_{14}$; and so
on. Finally, we find the parameters $\boldsymbol{\alpha}$ and
$\boldsymbol{\omega}$ in Equation 6 by maximizing the likelihood of the
computed pairs. After finding those parameters, it is straightforward to infer
$(\beta_{t+1},\gamma_{t+1})$ [13].

Note that we decided to estimate $\beta_{t+1}$ and $\gamma_{t+1}$ as a
function of the 3 previous values of those parameters; this allows them to
incorporate the velocity and acceleration at which the parameters change. That
is, we can estimate the velocity at which $\beta$ change as
$v_{\beta,t}=\beta_{t}-\beta_{t-1}$ and its acceleration as
$a_{\beta,t}=v_{\beta,t}-v_{\beta,t-1}$. Then, estimating
$\beta_{t}=\theta_{0}+\theta_{1}\beta_{t-1}+\theta_{2}v_{\beta,t-1}+\theta_{3}a_{\beta,t-1}$
is equivalent to the model in Equation 6. The same reasoning applies to the
computation of $\gamma_{t}$. We call this approach the “trend-following
varying-time parameters SIR”, tf-v-SIR.

For the case of $\mbox{CT}_{t}=-1$ and $\mbox{CT}_{t}=1$ (which represents a
change in trend from increasing number of infections to decreasing number of
infections or vice-versa), we set $\beta_{t+1}$ and $\gamma_{t+1}$ to values
such that the predicted number of new cases at week $t+1$ is identical to the
one at week $t$. We call this the “Same as the Last Observed Week” (SLOW)
model.

###  2.3 Estimating $\mbox{CT}_{t+1},\ \mbox{CP}_{t+1},\ \mbox{O}_{t}$

The random variables
$\mbox{CT}_{t+1}\in\\{-1,0,1\\},\mbox{CP}_{t+1}\in\\{-1,0,1\\}$ and
$\mbox{O}_{t}\in\\{0,1\\}$ in Figure 2 are all discrete nodes with discrete
parents, meaning their probability mass functions are fully defined by
conditional probability tables (CPTs). Learning the parameters of such CPTs
from data is challenging due to the scarcity of historical information. The
random variable $\mbox{CT}_{t+1}$ depends on the random variable changes in
policy (CP) at times $t-1,t-2,t-3$; however, there are very few changes in
policy in a given region, meaning it is difficult to accurately estimate those
probabilities from data. For the random variable O, which represents the
“willingness” of the government to implement a change in policy, there is no
observable data at all. We therefore relied on prior expert knowledge to set
the parameters of the conditional probability tables for these random
variables. Figure 3 shows the actual CPTs used in our experiments.

Figure 3:  Conditional probability tables used by siMLr. The names of the
variables refer to the nodes that appear on Figure 2

###  2.4 Estimating $U_{t}$

For modelling the random variable $U_{t}$, which represents the “Urgency to
change the trend”, we use an NN-CPD (neural-network conditional probability
distribution), which is a modified version of the multinomial logistic
conditional probability distribution [13].

######  Definition 2.1 (NN-CPD).

Let $Y\in\\{1,\dots,m\\}$ be an ${m}$-valued random variable with _k_ parents
$X_{1},\dots,X_{k}$ that each take on numerical values. The conditional
probability distribution $P(Y\ |\ X_{1},\dots,X_{k})$ is an NN-CPD if there is
an function $z\ =\ f_{\theta}(X_{1},\dots,X_{k})\ \in\mathbb{R}^{m}$,
represented as a neural network with parameters $\theta$, such that $p(Y=i\ |\
x_{1},\dots,x_{k})=\exp(z_{i})/\sum_{j}\exp(z_{j})$, where $z_{i}$ represents
the _i-th_ entry of $z$.

Similar to the case of discrete nodes with discrete parents, the main problem
in estimating the parameters $\theta$ of the NN-CPD is the scarcity of data.
To compensate for the lack of data, we again rely on domain knowledge and use
the probabilistic labels approach proposed by Vega et al. [21]. Using this
approach, we can use a dataset with few training instances along with their
probabilities to learn the parameters of a neural network more efficiently.

To compute $P(U_{t}\ |\ \mbox{SIR}_{t-2},\mbox{SIR}_{t-1},\ \mbox{SIR}_{t})$,
we extract two features: $c_{t}=\mbox{10E5}(S_{t-1}-S_{t})/N$, which
represents the number of new infections per 100K inhabitants; and
$v_{t}=c_{t}-c_{t-1}$, which estimates the rate of change of $c_{t}$, then
define $P(U_{t}\ |\ \mbox{SIR}_{t-2},\mbox{SIR}_{t-1},\
\mbox{SIR}_{t})=P(U_{t}\ |\ c_{t},v_{t})$.

To learn the parameters $\theta$, we created the dataset shown in Figure 4 and
a neural network with a single hidden layer with 64 units with sigmoid
activation functions, and 3 output units with softmax activation. As the
number and rate of change of new infections increases, the “urgency” to impose
stricter policies increases. The opposite effect (relaxing the policies)
occurs with a decrease in the number and rate of change of new infections.
When the new infections and its rate of change are relatively stable, then
there is a preference to keep the current policies in place.

Figure 4: Dataset used to create the NN-CPD for the variable $U_{t}$ and its
visualization. Values closer to 1 (yellow) increase $p(U_{t}=1\ |\
C_{t},V_{t})$. Values closer to 0 (green) increase $p(U_{t}=0\ |\
C_{t},V_{t})$. Values closer to $-1$ (blue) increase $p(U_{t}=-1\ |\
C_{t},V_{t})$

##  3 Experiments and results

We evaluated the performance of siMLr, in terms of MAPE, for forecasting the
number of new infections 1 to 4 weeks in advance, in data from United States
(as a country) and the 6 biggest provinces of Canada: Alberta (AB), British
Columbia (BC), Manitoba (MN), Ontario (ON), Quebec (QB), and Saskatchewan
(SK). For each of the regions, the predictions are done on a weekly basis,
over the 39 weeks from 26/Jul/2020 to 1/May/2021. This time span captures
different waves of infections.

At the end of every week, we fitted the siMLr parameters using the data that
was available until that week. For example, on 25/Jul/2020 we used all the
data available from 1/Jan/2020 to 25/Jul/2020 to fit the parameters of siMLr.
Then, we made the predictions for the new number of new infections during the
weeks: 26/Jul/2020 – 1/Aug/2020 (1 week in advance), 2/Aug/2020 – 8/Aug/2020
(2 weeks in advance), 9/Aug/2020 – 15/Aug/2020 (3 weeks in advance), and
16/Aug/2020 – 22/Aug/2020 (4 weeks in advance). After this, we then fitted the
parameters with data up to 1/Aug/2020 and repeated the same process, for 38
more iterations, until we covered the entire range of predictions.

We compared the performance of siMLr with the simple SIR compartmental model
with time-varying parameters learned using Equation 6 but no other random
variable (tf-v-SIR), and with the simple model that forecasts that the number
of cases 1 to 4 weeks in advance is the “Same as the Last Observed Week”
(SLOW). For the United States data, we also compared the performance of siMLr
against the publicly available predictions at the COVID-19 Forecast Hub
(https://github.com/reichlab/covid19-forecast-hub) [5].

For training, we used the publicly available (Creative Commons Attribution CC
BY standard) dataset OxCGRT [8], which contains the policies implemented by
different countries, as well as the time period over which they were
implemented. We limited our analysis to 3 policy decision: Workplace closing,
Stay at home requirements, and Cancellation of public events. For information
about the new number of reported cases and deaths, we used the publicly
available (Creative Commons Attribution 4.0 International) COVID-19 Data
Repository by the Center for Systems Science and Engineering at Johns Hopkins
University (https://github.com/CSSEGISandData/COVID-19) [6].

###  3.1 Data preprocessing

Before inputting the time-series data to siMLr, we performed some basic
preprocessing: (1) The original data contains the cumulative number of
reported infections/deaths on a daily basis. We trivially transformed this
time-series into the number of new daily infections/deaths. (2) We considered
negative values from the new daily infections/deaths time-series as missing,
assuming these negative values arose due to inconsistencies during the data
reporting procedure. (3) We “filled-in” the missing values. When the number of
new infections/deaths was missing at day $d$, we assumed that the entry at
$d+1$ contained the cases for both $d$ and $d+1$, and divided the number of
new infections/deaths evenly between both days. (4) We eliminated outliers.
For each day $d$, with number of reported new infections, $\Delta I_{d}$, we
computed the mean $(\mu_{d})$ and standard deviation $(\sigma_{d}$ of the set
$\Delta I_{d-10},\dots,\Delta I_{d-1}$; we then set $\Delta
I_{d}:=\min\\{\Delta I_{d},\mu_{d}+4\sigma_{d}\\}$. (5) We used the number of
new infections and new deaths to produce the SIR vector
$SIR_{t}=[S_{t},I_{t},R_{t}]$. Here, we assume everyone in that region was
susceptible at the start time – i.e., $S_{0}=$ total population. At each new
time, we transfer the number of new infections from $S$ to $I$, and the number
of new deaths and recovered from $I$ to $R$. If the number of new recovered
people is not reported, we assume that each new infected person that did not
die within the next 15 days, has recovered.

###  3.2 Results and discussion

Figure 5 shows the MAPE of the 1- to 4-week forecasts for the United States
(US) and the 6 biggest provinces of Canada. Note that siMLr has a consistently
lower MAPE than tf-v-SIR and SLOW.

Figure 6 illustrates the actual predictions of siMLr one week in advance (and
tf-v-SIR and SLOW) for the province of Alberta, Canada; and two weeks in
advance for the United States as a country. These two cases exemplify the
behaviour of siMLr.

As noted above, there is a 2- to 4-week lag after a policy changes, before we
see the effects. This means the task of making 1-week forecasts is relatively
simple, as the relevant policy (at times $t-3$ to $t-1$) is fully observable.
This allows siMLr to directly compute $\mbox{CT}_{t+1}$, which can then choose
whether to continue using the SIR with time-varying parameters if no policy
changed at time $t-1$, $t-2$, or $t-3$, versus use the SLOW predictor if the
policy changed. (See its CPT in Figure 3). Figure 6(a) shows a change in the
trend of reported new cases at week 22. However, just by looking at the
evolution of number of new infections before week 22, there is no way to
predict this change, which is why tf-v-SIR predicts that the number of new
infections will continue growing. However, since siMLr observed a change in
the government policies at week 20, it realized it could no longer rely on its
estimation of parameters and so switched to the SLOW model, which is why it
was more accurate here. A similar behaviour occurs in week 34, when the 3rd
wave of cases in Alberta started. Due to a relaxation in the policies on week
31, siMLr (at week 31) correctly predicted a change of trend around weeks
33–35.

Predictions at the country level are more complicated, since most of the time
policies are implemented at the state (or province) level instead of
nationally. For making predictions for an entire country, as well as
predictions 3 or 4 weeks in advance, siMLr first predicts, then uses, the
likelihood of observing a change in trend, at every week. In these cases, the
random variable $CT_{t+1}$ no longer acts like a “switch”, but instead it
mixes the predictions of the tf-v-SIR and SLOW models, according to the
probability of observing a change in the trend.

Figure 6(b) shows that whenever there is a stable trend in the number of new
reported infections – which suggests there has been no recent policy changes –
siMLr relies on the predictions of the tf-v-SIR model; however, as the number
(and rate of change) of new infections increases, so does the probability of
observing a change in the policy. Therefore, siMLr starts giving more weight
to the predictions of the SLOW model. Note this behaviour in the same figure
during weeks 13 – 20.

Figure 6(c) shows how our proposed siMLr approach compares with the 18 models
that submitted predictions to the CDC during the same span of time. Note that
siMLr and the model _LNQ-ens1_ are the best performing models, with no
statistically significant difference (p > 0.05 on a paired t-test) with
respect to MAPE. A striking result is how hard it is to beat the simple SLOW
model (presented as _COVIDhub-baseline_ in the CDC dataset). Out of the 19
models considered here, only 5 (including siMLr) do better than this simple
baseline when predicting 3 to 4 weeks ahead! This brings some insight into the
challenge of making accurate prediction in the medium term – probably due to
the need to predict, then use, policy change information.

##  4 Conclusions

Figure 5:  Comparison of siMLr, SIR model with time-varying parameters, and
SLOW. Figure 6:  (a) 1-week forecasts siMLr, tf-v-SIR, and SLOW, for Alberta,
Canada. (b) 2-week forecasts, of the same models, for US data. (c) Comparison
of siMLr versus models submitted to the CDC (on US data).

###  4.1 Limitations and societal impact

We do not foresee any potential negative societal impacts of this work. There
are, however, two main limitations to siMLr. siMLr’s underlying PGM makes it
easy to interpret. However, this structure does not imply causality: Although
changes in the observed policy influence changes in the trend of new reported
cases, the opposite is also true in reality. Also, while siMLr is able to
estimate the likelihood of a change in policy, it does not identify which
policy might be implemented next.

The second limitation is that siMLr relies on conditional probabilities tables
that are hard to learn due to lack of data, which forced us to build them
based on domain knowledge. If these estimations are disproportionately wrong,
then the predictions might be also misleading. Also, different regions might
have different “thresholds” for taking action. Although for our experiments we
set those CPTs a priori for US and Canada, they might differ for other
countries.

Despite these two limitations, siMLr produced state-of-the-art results in both
forecasting in the US as a country, and at the provincial level in Canada.

###  4.2 Contributions

This paper demonstrated that a model that explicitly models and incorporates
government policy decisions can accurately produce 1- to 4-week forecasts of
the number of COVID-19 infections. This involved showing (1) that an SIR model
with time-varying parameters is enough to describe the complex dynamics of
this pandemic, including the different waves of infections. Modelling siMLr as
a probabilistic graphical model not only makes it interpretable, but also
makes it easy to incorporate domain knowledge that compensates for the
relatively scarce data. siMLr’s excellent performance – comparable to state-
of-the-art systems in this competitive task – show that it is possible to
design interpretable machine learning models without sacrificing performance.

We expect that this approach will be useful not only for modelling COVID-19,
but other infectious diseases as well. We also hope that its interpretability
will leads to its adoption by researchers, and users, in epidemiology and
other non-ML fields.

## References

  * Anastassopoulou et al. [2020] C. Anastassopoulou, L. Russo, A. Tsakris, and C. Siettos.  Data-based analysis, modelling and forecasting of the covid-19 outbreak.  _PloS one_ , 15(3):e0230405, 2020.
  * Arik et al. [2020] S. Arik, C.-L. Li, J. Yoon, R. Sinha, A. Epshteyn, L. Le, V. Menon, S. Singh, L. Zhang, and M. Nikoltchev.  Interpretable sequence learning for covid-19 forecasting.  _NeurIPS_ , 33, 2020.
  * Blackwood and Childs [2018] J. C. Blackwood and L. M. Childs.  An introduction to compartmental modeling for the budding infectious disease modeler.  _Letters in Biomathematics_ , 5(1):195–221, 2018\.
  * Chen et al. [2020] Y.-C. Chen, P.-E. Lu, C.-S. Chang, and T.-H. Liu.  A time-dependent sir model for covid-19 with undetectable infected persons.  _IEEE Transactions on Network Science and Engineering_ , 7(4):3279–3294, 2020.
  * Cramer et al. [2021] E. Y. Cramer, V. K. Lopez, J. Niemi, G. E. George, J. C. Cegan, I. D. Dettwiller, W. P. England, M. W. Farthing, R. H. Hunter, B. Lafferty, et al.  Evaluation of individual and ensemble probabilistic forecasts of covid-19 mortality in the us.  _medRxiv_ , 2021.
  * Dong et al. [2020] E. Dong, H. Du, and L. Gardner.  An interactive web-based dashboard to track covid-19 in real time.  _The Lancet infectious diseases_ , 20(5):533–534, 2020.
  * Fauci et al. [2020] A. S. Fauci, H. C. Lane, and R. R. Redfield.  Covid-19—navigating the uncharted, 2020.
  * Hale et al. [2021] T. Hale, N. Angrist, R. Goldszmidt, B. Kira, A. Petherick, T. Phillips, S. Webster, E. Cameron-Blake, L. Hallas, S. Majumdar, et al.  A global panel database of pandemic policies (oxford covid-19 government response tracker).  _Nature Human Behaviour_ , 5(4):529–538, 2021\.
  * Jacobs et al. [1991] R. A. Jacobs, M. I. Jordan, S. J. Nowlan, and G. E. Hinton.  Adaptive mixtures of local experts.  _Neural computation_ , 3(1):79–87, 1991.
  * Jin et al. [2021] X. Jin, Y.-X. Wang, and X. Yan.  Inter-series attention model for covid-19 forecasting.  In _SIAM International Conference on Data Mining (SDM)_ , pages 495–503. SIAM, 2021.
  * Kafieh et al. [2021] R. Kafieh, R. Arian, N. Saeedizadeh, Z. Amini, N. D. Serej, S. Minaee, S. K. Yadav, A. Vaezi, N. Rezaei, and S. Haghjooy Javanmard.  Covid-19 in iran: forecasting pandemic using deep learning.  _Computational and Mathematical Methods in Medicine_ , 2021, 2021\.
  * Kermack and McKendrick [1927] W. O. Kermack and A. G. McKendrick.  A contribution to the mathematical theory of epidemics.  _Proceedings of the royal society of london. Series A, Containing papers of a mathematical and physical character_ , 115(772):700–721, 1927.
  * Koller and Friedman [2009] D. Koller and N. Friedman.  _Probabilistic graphical models: principles and techniques_.  MIT press, 2009.
  * Liao et al. [2020] Z. Liao, P. Lan, Z. Liao, Y. Zhang, and S. Liu.  Tw-sir: time-window based sir for covid-19 forecasts.  _Scientific reports_ , 10(1):1–15, 2020.
  * Liu and Manzoor [2020] H. Liu and A. e. a. Manzoor.  The covid-19 outbreak and affected countries stock markets response.  _Int’l J Environmental Research and Public Health_ , 17(8):2800, 2020.
  * Liu and Stechlinski [2012] X. Liu and P. Stechlinski.  Infectious disease models with time-varying parameters and general nonlinear incidence rate.  _Applied Mathematical Modelling_ , 36(5):1974–1994, 2012.
  * Mojjada et al. [2020] R. K. Mojjada, A. Yadav, A. Prabhu, and Y. Natarajan.  Machine learning models for covid-19 future forecasting.  _Materials Today: Proceedings_ , 2020.
  * Murphy [2012] K. P. Murphy.  _Machine learning: a probabilistic perspective_.  MIT press, 2012.
  * Omran et al. [2021] N. F. Omran, S. F. Abd-el Ghany, H. Saleh, A. A. Ali, A. Gumaei, and M. Al-Rakhami.  Applying deep learning methods on time-series data for forecasting covid-19 in egypt, kuwait, and saudi arabia.  _Complexity_ , 2021, 2021.
  * Rudin [2019] C. Rudin.  Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead.  _Nature Machine Intelligence_ , 1(5):206–215, 2019.
  * Vega et al. [2021] R. Vega, P. Gorji, Z. Zhang, X. Qin, A. Rakkunedeth, J. Kapur, J. Jaremko, and R. Greiner.  Sample efficient learning of image-based diagnostic classifiers via probabilistic labels.  In _International Conference on Artificial Intelligence and Statistics_ , pages 739–747. PMLR, 2021.
  * Watson et al. [2021] G. L. Watson, D. Xiong, L. Zhang, J. A. Zoller, J. Shamshoian, P. Sundin, T. Bufford, A. W. Rimoin, M. A. Suchard, and C. M. Ramirez.  Pandemic velocity: Forecasting covid-19 in the us with a machine learning & bayesian time series compartmental model.  _PLoS computational biology_ , 17(3):e1008837, 2021.
  * Yeung et al. [2021] A. Y. Yeung, F. Roewer-Despres, L. Rosella, and F. Rudzicz.  Machine learning–based prediction of growth in confirmed covid-19 infection cases in 114 countries using metrics of nonpharmaceutical interventions and cultural dimensions: Model development and validation.  _Journal of Medical Internet Research_ , 23(4):e26628, 2021.