#  _You better watch out_ : US COVID-19 wave dynamics versus vaccination
strategy

Giacomo Cacciapaglia<EMAIL_ADDRESS>Institut de Physique des
deux Infinis de Lyon (IP2I), UMR5822, CNRS/IN2P3, F-69622, Villeurbanne,
France University of Lyon, Université Claude Bernard Lyon 1, F-69001, Lyon,
France Corentin Cot<EMAIL_ADDRESS>Institut de Physique des deux Infinis
de Lyon (IP2I), UMR5822, CNRS/IN2P3, F-69622, Villeurbanne, France University
of Lyon, Université Claude Bernard Lyon 1, F-69001, Lyon, France Anna Sigridur
Islind<EMAIL_ADDRESS>Department of Computer Science, Reykjavík University,
Menntavegur 1, 102 Reykjavík, Iceland María Óskarsdóttir<EMAIL_ADDRESS>Department of Computer Science, Reykjavík University, Menntavegur 1, 102
Reykjavík, Iceland Francesco Sannino<EMAIL_ADDRESS>CP3-Origins & the
Danish Institute for Advanced Study, University of Southern Denmark, Campusvej
55, DK-5230 Odense, Denmark  Dipartimento di Fisica E. Pancini, Università di
Napoli Federico II & INFN sezione di Napoli, Complesso Universitario di Monte
S. Angelo Edificio 6, via Cintia, 80126 Napoli, Italy.

###### Abstract

Abstract: We employ the epidemic Renormalization Group (eRG) framework to
understand, reproduce and predict the COVID-19 pandemic diffusion across the
US. The human mobility across different geographical US divisions is modelled
via open source flight data alongside the impact of social distancing for each
such division. We analyse the impact of the vaccination strategy on the
current pandemic wave dynamics in the US. We observe that the ongoing
vaccination campaign will not impact the current pandemic wave and therefore
strict social distancing measures must still be enacted. To curb the current
and the next waves our results indisputably show that vaccinations alone are
not enough and strict social distancing measures are required until sufficient
immunity is achieved. Our results are essential for a successful vaccination
strategy in the US.

##  I Introduction

The United States (US), raged by the SARS-CoV-2 virus, are paying an immense
toll in terms of the loss of human lives and jobs, with a dreadful impact on
society and economy. Understanding and predicting the time evolution of the
pandemic plays a key role in defining prevention and control strategies.
Short-term forecasts have been obtained, since the early days, via effective
methods Perc  _et al._ (2020); Hâncean  _et al._ (2020); Zhou  _et al._
(2020). Furthermore, time-honored mathematical models can be used, like
compartmental models Prem  _et al._ (2020); Scala  _et al._ (2020); Friston
_et al._ (2020); Sonnino  _et al._ (2020); Abou-Ismail (2020) of the SIR type
Kermack  _et al._ (1927) or complex networks Zhan  _et al._ (2018); Perc  _et
al._ (2017); Wang  _et al._ (2015). Nevertheless, it remains very hard to
understand and forecast the wave pattern of pandemics like COVID-19 Scudellari
(2020).

In this work, we employ the _epidemic Renormalization Group_ (eRG) framework,
recently developed in Della Morte  _et al._ (2020); Cacciapaglia and Sannino
(2020a). It can be mapped Cacciapaglia and Sannino (2020a); Della Morte and
Sannino (2020) into a time-dependent compartmental model of the SIR type
Kermack  _et al._ (1927). The eRG framework provides a single first order
differential equation, apt to describing the time-evolution of the cumulative
number of infected cases in an isolated region Della Morte  _et al._ (2020).
It has been extended in Cacciapaglia and Sannino (2020a) to include
interactions among multiple regions of the world. The main advantage over SIR
models is its simplicity, and the fact that it relies on symmetries of the
system instead of a detailed description. As a result, no computer simulation
is needed in order to understand the time-evolution of the epidemic even at
large scales Cacciapaglia and Sannino (2020a). Recently, the framework has
been extended to include the multi-wave pattern Cacciapaglia and Sannino
(2020b); Cacciapaglia  _et al._ (2020a) observed in the COVID-19 and other
pandemics Taubenberger and Morens (2006).

The Renormalization Group approach Wilson (1971a, b) has a long history in
physics with impact from particle to condensed matter physics and beyond. Its
application to epidemic dynamics is complementary to other approaches Li  _et
al._ (2019); Zhan  _et al._ (2018); Perc  _et al._ (2017); Wang  _et al._
(2015, 2016); Danby (1985); Brauer (2019); Miller (2012); Murray (2002);
Fishman  _et al._ (2014); Pell  _et al._ (2018). Here we demonstrate that the
framework is able to reproduce and predict the pandemic diffusion in the US
taking into account the human mobility across different geographical US
divisions, as well as the impact of social distancing within each one. To gain
an insight and to better monitor the human exchange we make use of open source
flight data among the states. We calibrate the model on the first wave
pandemic, raging from March to August, 2020. With these insights, we then
analyse and understand the current second wave, raging in all the divisions.
The eRG framework can also be easily adapted to take into account vaccinations
Wang  _et al._ (2016). We propose a new framework and use it to quantify the
impact of the vaccination campaign, started on December 14th, on the current
and future wave dynamics. Our results are in agreement with previous work
based on compartmental models Paltiel  _et al._ (2020), and confirm that the
current campaign will have limited impact on the ongoing wave.

##  II Methodology

In this section we briefly review our methods that include the open source
flight data description, their interplay with the eRG mathematical model
framework and, last but not least, the interplay with vaccine deployment and
implementation.

###  II.1 Data description

The flight data comes from the OpenSky Network, which is a non-profit
association that provides open access to real-world air traffic control
dataset for research purposes Schäfer  _et al._ (2014). The OpenSky COVID-19
Flight Dataset 111opensky-network.org was made available in April 2020 and is
currently updated on a monthly basis, with the purpose of supporting research
on the spread of the pandemic and the associated economic impact. This dataset
has been used to investigate mobility in the early months of the pandemic
Islind  _et al._ (2020) as well as the pandemic’s effect on economic
indicators Bank of England (2020).

The data provides information about the origin and destination airports as
well as the date and time of all flights worldwide. For our analyses we
considered domestic flights in the US only. We aggregated the data, to obtain
the number of flights between all pairs of airports per day, from the
beginning of April until the end of October, 2020. Subsequently, the airports
in each state and the number of flights associated with them were combined, to
give the number of within and between state flights, on a day to day basis for
the whole period.

The number of daily infected cases, which is also used for analysis in this
paper, is provided by the open source online repository Opendatasoft
222public.opendatasoft.com/explore/dataset/testing-data-covid19-usa/.

Figure 1: Illustration of the geographical divisions of the US used in this
study.

###  II.2 Mathematical modeling

The states within the US have different population and demographic
distribution. A state-by-state mathematical modeling, therefore, is challenged
by statistical artifacts. For these reasons we group the states following the
census divisions 333US Census Bureau, as summarized in Table 1 and illustrated
in Fig.1. Note, that contrary to the official definitions, we include Maryland
and Delaware in Mid-Atlantic instead of South Atlantic. The main reason is
that the population of these two states is more connected to states in Mid-
Atlantic, as proven by the diffusion timing of the virus.

Division composition
---
Division names |  Division code |  States within the division
New England |  NE |  Massachusetts, Connecticut, New Hampshire, Maine, Rhode Island and Vermont
Mid-Atlantic |  MA |  New York, Pennsylvania, New Jersey, Maryland and Delaware
South Atlantic |  SA |  Florida, Georgia, North Carolina, Virginia, South Carolina and West Virginia
East South Central |  ESC |  Tennessee, Alabama, Kentucky and Mississippi
West South Central |  WSC |  Texas, Louisiana, Oklahoma and Arkansas
East North Central |  ENC |  Illinois, Ohio, Michigan, Indiana and Wisconsin
West North Central |  WNC |  Missouri, Minnesota, Iowa, Kansas, Nebraska, South Dakota and North Dakota
Mountains |  M |  Arizona, Colorado, Utah, Nevada, New Mexico, Idaho, Montana and Wyoming
Pacific |  P |  California, Washington, Oregon, Hawaii and Alaska
Table 1: States of the US integrated into 9 divisions. Maryland and Delaware
are moved from South Atlantic to Mid-Atlantic.

Building upon our successful understanding of the COVID-19 temporal evolution
Cacciapaglia  _et al._ (2020b) we apply our framework to the US case. Building
on that framework we employ the following eRG set of first order differential
equations Cacciapaglia and Sannino (2020a) to describe the time-evolution of
the cumulative number of infected cases within the US divisions:

$\displaystyle\frac{d\alpha_{i}}{dt}=\gamma_{i}\alpha_{i}\left(1-\frac{\alpha_{i}}{a_{i}}\right)+\sum_{j\neq
i}\frac{k_{ij}}{n_{mi}}(e^{\alpha_{j}-\alpha_{i}}-1)\,,$ (1)

where

$\alpha_{i}(t)=\rm ln\ \mathcal{I}_{i}(t)\ ,$ (2)

with $\mathcal{I}_{i}(t)$ being the cumulative number of infected cases per
million inhabitants for the division $i$ and $\ln$ indicating its natural
logarithm. These equations embody, within a small number of parameters, the
pandemic spreading dynamics across coupled regions of the world via the
temporal evolution of $\alpha_{i}(t)$. The parameters $\gamma_{i}$ and $a_{i}$
can be extracted by the data within each single wave. The fit methodology is
described in Della Morte  _et al._ (2020); Cacciapaglia and Sannino (2020a).

Figure 2: Weekly new number of cases for all the 9 divisions.

In the US, it is well known that the COVID-19 pandemic started in NE and MA
(mainly in New York City) and then spread to the other divisions. Thus, we
define the _US first wave_ period from March to the end of August as shown in
Fig. 2. In particular, one observes a peak of new infected in NE and MA around
April, while for the other divisions the main peak occurs around July. We also
observe an initial feature in the latter divisions that we did not attempt to
model except for ENC (mostly located in Chicago) and WNC. For the two latter
divisions, we considered these as two independent first wave components. The
_US second wave_ is thus associated with the episode starting in October,
2020.

As a first method, working under the assumption that the US pandemic indeed
originated in New York (MA), we first determine the $k_{ij}$ matrix entries
between the division MA and the others. The values are chosen to reasonably
reproduce the delay between the main peaks of the first wave in pairs of
divisions (C.f. the top section in Table 2). Interestingly, with the exception
of NE, the entries of the $k$ matrix are comparable to the ones we used for
Europe Cacciapaglia  _et al._ (2020b). For NE, a large coupling is needed due
to the tight connections between the two regions, in particular New York City
with the neighbouring states and Massachusetts.

As a second method, we used the flight data to estimate the number of
travellers between different divisions, under the assumption that the $k_{ij}$
matrix entries are proportional to this set of data. To have a realistic
matrix for $k_{ij}$, we first take the mean number of flights from division
$i$ to division $j$ during the period from April 1st to May 31st for the first
wave, and from September 1st to October 31st for the second wave. Then, we
multiply the number of flights by an effective average number of passengers,
and normalize it by $10^{6}$, following the definition of $k_{ij}$
Cacciapaglia and Sannino (2020a). For the first wave, the optimal average
number of passenger is found to be $10$, while for the second wave we find an
optimal value of $5$. Note that these values do not correspond to the actual
number of passengers in the flights: in fact, the values of the couplings
$k_{ij}$ also take into account the probability of the passengers to carry the
infection as compared to the average in the division of origin. When the value
is low it might suggest that the sample of passengers in a flight is less
infectious than average, as people with symptoms tend not to travel. Controls
at airports may also contribute to this. The key information we extract from
the flight data is the relative flux of infections among different divisions.

The results are listed in the middle and bottom sections of Table 2. We keep
the same value from the previous fit only for MA-NE. The reason behind this
choice is the tight connection between the two divisions, where most of the
human mobility is imputable to road transport.

$k_{ij}$ values (1st wave fits)
---
Division code |  NE |  MA |  SA |  ESC |  WSC |  ENC |  WNC |  M |  P
MA |  $0.72$ |  $0$ |  $0.0014$ |  $0.00075$ |  $0.0017$ |  $0.0023$ |  $0.0005$ |  $0.002$ |  $0.0053$
First wave $k_{ij}$ values (Flight data, from April 1st to May 31st)
Division code |  NE |  MA |  SA |  ESC |  WSC |  ENC |  WNC |  M |  P
NE |  $0$ |  $0.72$ |  $0.0045$ |  $0.00088$ |  $0.00087$ |  $0.0024$ |  $0.00052$ |  $0.00067$ |  $0.00091$
MA |  $0.72$ |  $0$ |  $0.019$ |  $0.0056$ |  $0.0041$ |  $0.012$ |  $0.0025$ |  $0.0031$ |  $0.0059$
SA |  $0.0043$ |  $0.018$ |  $0$ |  $0.0085$ |  $0.013$ |  $0.019$ |  $0.0057$ |  $0.0050$ |  $0.0067$
ESC |  $0.00092$ |  $0.0053$ |  $0.0093$ |  $0$ |  $0.0051$ |  $0.0068$ |  $0.0023$ |  $0.0035$ |  $0.0065$
WSC |  $0.00095$ |  $0.0038$ |  $0.014$ |  $0.0055$ |  $0$ |  $0.0092$ |  $0.0054$ |  $0.011$ |  $0.010$
ENC |  $0.0025$ |  $0.012$ |  $0.018$ |  $0.0063$ |  $0.0086$ |  $0$ |  $0.0082$ |  $0.0079$ |  $0.0099$
WNC |  $0.00038$ |  $0.0022$ |  $0.0056$ |  $0.0019$ |  $0.0046$ |  $0.0070$ |  $0$ |  $0.0055$ |  $0.0027$
M |  $0.00050$ |  $0.0020$ |  $0.0042$ |  $0.0026$ |  $0.011$ |  $0.0072$ |  $0.0043$ |  $0$ |  $0.028$
P |  $0.00084$ |  $0.0055$ |  $0.0063$ |  $0.0050$ |  $0.010$ |  $0.0092$ |  $0.0033$ |  $0.030$ |  $0$
Second wave $k_{ij}$ values (Flight data, from September 1st to October 31st)
Division code |  NE |  MA |  SA |  ESC |  WSC |  ENC |  WNC |  M |  P
Region-X |  $0.0066$ |  $0.028$ |  $0.029$ |  $0.013$ |  $0.019$ |  $0.027$ |  $0.014$ |  $0.03$ |  $0.03$
NE |  $0$ |  $0.72$ |  $0.0028$ |  $0.00046$ |  $0.00031$ |  $0.0015$ |  $0.00026$ |  $0.00041$ |  $0.00082$
MA |  $0.72$ |  $0.$ |  $0.011$ |  $0.002$ |  $0.0017$ |  $0.0064$ |  $0.0013$ |  $0.0021$ |  $0.0029$
SA |  $0.0026$ |  $0.011$ |  $0$ |  $0.005$ |  $0.005$ |  $0.0096$ |  $0.003$ |  $0.0033$ |  $0.0035$
ESC |  $0.00041$ |  $0.0019$ |  $0.0051$ |  $0$ |  $0.0019$ |  $0.0028$ |  $0.00087$ |  $0.0012$ |  $0.0015$
WSC |  $0.00028$ |  $0.0015$ |  $0.0049$ |  $0.0018$ |  $0$ |  $0.0028$ |  $0.0016$ |  $0.004$ |  $0.0034$
ENC |  $0.0014$ |  $0.0062$ |  $0.0089$ |  $0.0028$ |  $0.003$ |  $0$ |  $0.0039$ |  $0.0043$ |  $0.0045$
WNC |  $0.00024$ |  $0.0013$ |  $0.0028$ |  $0.0009$ |  $0.0017$ |  $0.0038$ |  $0$ |  $0.0028$ |  $0.0016$
M |  $0.00032$ |  $0.0017$ |  $0.0029$ |  $0.0011$ |  $0.0054$ |  $0.004$ |  $0.0026$ |  $0$ |  $0.014$
P |  $0.00074$ |  $0.0028$ |  $0.0032$ |  $0.0014$ |  $0.0046$ |  $0.0041$ |  $0.0018$ |  $0.015$ |  $0$
Table 2: Values of the $k_{ij}$ entries among US divisions. In the top
section, the values between Mid-Atlantic (MA) and the other divisions are
obtained from fits of the first wave timing. In the central and bottom
sections, the complete matrix (except the entries between MA and NE) is
obtained using flight data for the first wave (from April 1st to May 31st) and
the second (from September 1st to October 31st), respectively.

By the end of November, we clearly observe a new rise in the number of
infections, signalling the onset of a second wave pandemic in the US (see Fig.
2). Using our framework, we model and then simulate the second wave across the
different US divisions.

Finally, to check the geographical diffusion of the virus during the various
phases of the pandemic in the US, we define an indicator of the uniformity of
the new case incidence Cacciapaglia  _et al._ (2020a). This indicator can be
defined week by week via a $\chi^{2}$-like variable, given by:

$\chi^{2}(t)=\frac{1}{9}\sum^{9}_{i=1}\left(\frac{\mathcal{I}_{i}^{\prime}(t)}{\left<\mathcal{I}^{\prime}(t)\right>}-1\right)^{2}\,,$
(3)

where $\mathcal{I}_{i}^{\prime}(t)$ is the number of new cases per week in
division $i$ at time $t$ and $\left<\mathcal{I}^{\prime}(t)\right>$ the mean
of the same quantity in the 9 divisions. The parameter $\chi^{2}$ quantifies
the geographical diffusion of the SARS-CoV-2 virus in the US: the smaller its
value, the more uniform the pandemic spread within the whole country. The
result is shown in Fig. 3: during the first peak in April (light gray shade),
the value of $\chi^{2}$ is large, signalling that the epidemic diffusion is
localized in a few divisions; during the second peak of the first wave (gray
shade), the value has dropped, signalling that the epidemic has been spreading
to all divisions. Finally, the data for the ongoing second wave (dark gray
shade) shows that $\chi^{2}$ is dropping towards zero, as expected for a more
diffuse incidence of infections.

Figure 3: Evolution of the uniformity indicator $\chi^{2}$ over time (weekly
basis). The shaded bands indicate the period when epidemic peaks are recorded.

###  II.3 Vaccine deployment and implementation

Various vaccines have been developed for the COVID-19 pandemic, and their
deployment in the US has already started on December 14th
444https://www.washingtonpost.com. The effect of the immunization due to the
vaccine has been studied in the context of compartmental models, like SEIR
Paltiel  _et al._ (2020). In our mathematical model, the simplest and most
intuitive effect is a reduction of both the total number of infections during
a single wave, $a_{i}$, and/or the effective diffusion rate of the virus
$\gamma_{i}$, in each division.

To validate this working hypothesis, and understand how the vaccination of a
portion of the population affects the values of $a$ and $\gamma$ in the eRG
framework, we studied the effect of immunization in a simple percolation
model, which has been shown to be in the same class of universality as simple
compartmental models Cardy and Grassberger (1985). To do so, we set up a
Monte-Carlo simulation, consisting of a square grid whose nodes are associated
to a susceptible individual. Each node can be in four exclusive states:
Susceptible (S), Infected (I), Recovered (R) or Vaccinated (V). At each step
in time in the simulation, for each node we generate a random number $r$
between 0 and 1: if the node is in state S in proximity to a node in state I
and $r<\gamma_{\ast}$, we switch its state to I, else it remains S; if the
node is in state I and $r<\epsilon_{\ast}$, we switch its state to R, else it
remains I; if the node is in state R or V, it will not change. This model
reproduces the diffusion of the infection, where $\gamma_{\ast}$ is the
infection probability on the lattice and $\epsilon_{\ast}$ is the recovery
rate. Finally, we fit the data from the simulation to the solution of a simple
eRG equation to extract $\gamma$ and $a$. The vaccination is implemented by
setting a random fraction $R_{v}$ of nodes to the state V before the
simulation starts. The values of $a$ and $\gamma$ as a function of the
fraction of vaccinations are shown in Fig.4: we observe that both parameters
are reduced by the same percentage as the vaccination up to $R_{v}\lesssim
25$%. Above this value of vaccinated nodes, the simulation is unstable and the
result cannot be trusted. This result, nevertheless, demonstrates that the
vaccination reduces both parameters $a$ and $\gamma$ proportionally
reinforcing our expectation.

Figure 4: $a$ and $\gamma$ fit parameters versus initial percentage of
vaccinated nodes for $\gamma_{\ast}=0.6$ and $\epsilon_{\ast}=0.4$.

In a realistic scenario, the vaccination of the population can only be
implemented in a gradual way, so that the total vaccination campaign has a
duration in time. We can thereby assume that a fraction $R_{v}$ of the
population is vaccinated in a time interval $\Delta_{t}$. The rate of
vaccinations is therefore $c=R_{v}/\Delta_{t}$. This implies that the
variation in $\gamma$, during the time interval from $t_{v}$ to
$t_{v}+\Delta_{t}$, is given by:

$\frac{d\gamma(t)}{dt}=-c\ \gamma(t_{v})\,,$ (4)

where $\gamma(t_{v})$ is the effective infection rate before the start of the
vaccination campaign. The solution for the time-dependent effective infection
rate is

$\displaystyle\gamma(t)$ $\displaystyle=$
$\displaystyle\gamma(t_{v})[1-c(t-t_{v})]\,,$ (5)

until $t=t_{v}+\Delta_{t}$, after which $\gamma$ remains constant again at a
reduced value $\gamma(t_{v})\ (1-R_{v})$.

To find the variation of $a(t)$ within the vaccination interval $t_{v}$ to
$t_{v}+\Delta_{t}$, we assume that the not-yet-infected individuals are
vaccinated at the same rate $c$ as the total population. Thus, at any given
time, the variation in the number of individuals that will be exposed to the
infection, $e^{a(t)}$, is proportional to the difference
$e^{a(t)}-e^{\alpha(t)}$. This leads to the following differential equation:

$\frac{da(t)}{dt}=-c\ (1-e^{\alpha(t)-a(t)})\,.$ (6)

This equation needs to be solved in a coupled system with the eRG one. Note
that the derivative is zero outside of the time interval
$[t_{v},t_{v}+\Delta_{t}]$. In the numerical solutions for the effect of the
vaccine, we will add one equation for each $a_{i}(t)$, assuming that the
vaccination rate $c$ is the same in all divisions.

| First wave parameters (fitted) | Second wave parameters
---|---|---
Division | Code | $a$ | $\gamma$ | $a$ | $\gamma$
New England |  NE |  $9.397(7)$ |  $0.416(7)$ |  $11.006$ |  $0.214$
Mid-Atlantic |  MA |  $9.496(7)$ |  $0.516(9)$ |  $10.882$ |  $0.206$
South Atlantic |  SA |  $9.691(5)$ |  $0.370(3)$ |  $10.885$ |  $0.185$
East South Central |  ESC |  $9.63(2)$ |  $0.331(6)$ |  $11.201$ |  $0.207$
West South Central |  WSC |  $9.720(6)$ |  $0.340(3)$ |  $10.713$ |  $0.213$
East North Central |  ENC |  $8.88(3)$ |  $0.300(8)$ |  $11.074$ |  $0.250$
West North Central |  WNC |  $8.66(2)$ |  $0.342(6)$ |  $11.060$ |  $0.263$
Mountain |  M |  $9.478(9)$ |  $0.330(4)$ |  $11.089$ |  $0.213$
Pacific |  P |  $9.33(1)$ |  $0.291(5)$ |  $11.535$ |  $0.171$
Table 3: Parameters of the eRG model for the first and second wave in the 9
divisions. For the first wave, we report the values from the fit, including
the $1\sigma$ error. For the second wave, the values are chosen to reproduce
the current data, adjourned to December 16th.

##  III Results

###  III.1 Validating the eRG on the first wave data

The epidemic data (C.f. Fig.2) shows that the MA division (New York City) was
first hit hard by the COVID-19 pandemic, and was followed closely by NE. The
other divisions witnessed a comparable peak of new infections 3-4 months
later. Note that we are using cases normalized per million to facilitate the
comparison between divisions with varying population. As a first study, we
want to test the $eRG$ equations (1) against the hypothesis that the epidemic
has been diffusing from MA to the other divisions. The parameters $a_{i}$ and
$\gamma_{i}$ are fixed by fitting the data, as shown in Table 3. Thus, the
timing of the peaks in the divisions is determined by the entries of the
$k_{ij}$ matrix. Determining all 81 entries from the data is not possible, as
we only have 9 epidemiological curves. Thus, we assume that only the couplings
between the source MA and any other division are responsible. The results of
the fits are shown in the top block of Table 2, and will be used as a control
benchmark.

Except for $k_{21}$ that links NE and MA, all the other $k_{2j}$ are of order
$10^{-3}$, thus confirming the range we found for the European second wave
Cacciapaglia  _et al._ (2020b). The value of $k_{21}$ is of order unity, which
implies that there is a stronger connection between the two divisions. This
may be explained by the fact that there exist a significant flow of people
between New York City and the neighbouring states (including Massachusetts) in
New England. Work commutes and weekend travelling by car explains the required
high number of travellers per week. Another interesting feature is the
presence of a small peak of infections for ENC and WNC, around March. This
feature cannot originate from the MA division, as that would imply a $k$-value
of order $10$, which is clearly unrealistic Cacciapaglia and Sannino (2020a).
The only viable solution is that the epidemic hit these two divisions from
abroad. On the other hand, the second peak observed around August can be
explained by the interaction with MA.

The values of $k_{ij}$ are, in principle, determined by the flow of people
between different divisions. Thus, we could use any set of mobility data Yang
_et al._ (2020) to estimate the relative numbers of the entries, while the
normalization also depends on the effective infection power of the traveling
individuals and it can be determined from the data. With the help of mobility
data, we can reduce the 81 parameters to a single one. Due to the large
distances across divisions, we decided to focus on the flight data, as
described in the methodology section. The values of the entries are reported
in the middle section of Table 2. Note that for MA-NE we used the same value
obtained from the previous fit, as the people’s flow is mostly dominated by
land movements.

xxxx First wave |  xxx Second wave
---|---
|
|
|
Figure 5: Simulation of the spread of the first wave (left plots) and the
second wave (right plots) using flight-data-derived kappa matrix. For the
first wave, MA is used as a seed region, while for the second wave a
combination of the first waves among divisions acts as the seed region
(Region-X). The vertical dashed lines in the right column plots mark the date
when the simulation was done. The data points in the grayed region where not
used to tune the eRG solutions.

Using this matrix of $k_{ij}$ to simulate the spread of the first wave across
the country, as originating from MA, we obtain the curves in the left panels
of Fig.5. For nearly all divisions, we obtain the correct timing for the peak,
with the exception of SA and ESC (for ENC and WNC, the anomaly may be linked
to the presence of a mild early peak and the absence of a prominent second
peak). The results are more accurate for divisions far from MA, thus
validating the method as the diffusion of the virus seems to depend on the
people travelling (by air) among divisions. For SA, the predicted curve is
substantially anticipated compared to the data: this discrepancy may be
explained by the presence of an air hub in Atlanta, GE, so that many of the
passengers of flights landing there do not stop in the division but instead
take an immediate connecting flight.

###  III.2 Understanding the second wave

The US states are currently witnessing a second wave, which is ravaging in all
the 9 divisions with comparable intensity. Previous studies in the eRG
framework have uncovered two possible origins for an epidemic wave to start:
one is the coupling with an external region with a raging epidemic
Cacciapaglia and Sannino (2020a), the second is the instability represented by
a strolling phase in between waves Cacciapaglia and Sannino (2020b);
Cacciapaglia  _et al._ (2020a). We have shown that the former mechanism can
account for the peak structure during the first wave.

As a first step, we will try to use the same method to understand the second
wave. Since travelling to the US from abroad has been strongly reduced and
regulated, we will consider the divisions that witnessed a peak in July-August
as source for the second wave. To this purpose, we define a Region-X
Cacciapaglia and Sannino (2020a) as an average sum of all the divisions with a
pandemic peak occurring in the July-August period. The parameters are chosen
to reproduce the number of cases in the totality of the relevant 7 divisions
(SA, ESC, NSC, ENC, WNC, M and P) normalized by the total population. For each
division, we optimized $a_{i}$ and $\gamma_{i}$ to reproduce the current data
adjourned at December 16th (C.f. Table 3). For the couplings $k_{ij}$ we use
the flight data, except for the usual MA-NE couplings (C.f. bottom section of
Table 2). Finally, the $k_{0j}$ connecting the 9 divisions to the source
Region-X are computed by summing the $k$ entries between the division $j$ and
the 7 divisions used to model Region-X (also derived from flight data).

The results of the eRG equations are shown in the right panels of Fig. 5,
showing a good agreement. The values of the $k_{0j}$ of the Region-X are one
order of magnitude larger than the others. This fact can be interpreted by the
presence of hotspots in each division which also contribute significantly to
the new wave. In other words, traveling among divisions cannot be the only
responsible factor for the onset of the second wave in the US. This hypothesis
can also be validated by studying the uniformity of the distribution of the
new infections in various states during the three peaks, as shown in Fig.3.
Comparing the three peak regions, we see that the uniformity indicator is
systematically decreasing, thus indicating a more geographically uniform
presence of the virus.

It is also interesting to notice that the value of $\gamma_{i}$ for the second
wave is systematically smaller than the infection rate during the first wave.
This is in agreement with the results we found in Cacciapaglia and Sannino
(2020b); Cacciapaglia  _et al._ (2020a), where we modelled the multi-wave
structure of the pandemic via an instability inside each region. The result of
this simple analysis supports the hypothesis that the virus is now endemic for
all states in the US, thus a multi-wave pattern will continue to emerge.
Traveling among states (or divisions) is less relevant at this stage.

The result of our eRG analysis shows that the current wave will end in March-
April 2021. Note, however, that we have not taken into account the potential
disastrous effect of the Christmas and New Year holidays, which could lead to
an increase in the infection rates. In some divisions there is a increase at
the end of November, which can be attributed to the Thanksgiving holiday.

###  III.3 Effect of the current vaccination strategy

xxxx South Atlantic |  xxxxWest North Central
---|---
|
Figure 6: Evolution of the number of infections without vaccination ($c=0$)
and with a vaccination rate of $0.64\%$/week, $1\%$/week and $2\%$/week
starting on December 14th and stopping at $20\%$ of the population vaccinated.
We show the results for two sample divisions: South Atlantic and West North
Central.

Following the development of multiple vaccines for the SARS-COV-2 virus
555https://www.nature.com/articles/d41586-020-03626-1, vaccination campaigns
have started in many countries, including the US. This will influence the
development of the current wave, and help in curbing the future ones. The
vaccination campaign started on December 14th in the US
666https://www.washingtonpost.com. We also know that the US has purchased 100
million doses from Pfizer (plus an additional 100 million from Moderna)
777https://www.forbes.com, so that at least $20\%$ of the population may be
vaccinated in this first campaign. As of December 28, $0.64\%$ of the total
population has been vaccinated 888https://ourworldindata.org in a little over
one week, thus in our study we will use this as a benchmark weekly rate. The
data listed above defines our starting benchmark for the current vaccination
campaign.

To study the effect of the vaccinations, we have solved the eRG equations for
the second wave, with the addition of the reduction of $a_{i}$ and
$\gamma_{i}$, as detailed in the methodology section. We show the result for
two sample divisions in Fig.6 (dashed curves) as compared to the same
solutions without vaccines (solid curves). A vaccination at a $0.64\%$ rate
per week does not affect the peak of new infections. As a reference, we also
increased the vaccination rates to $1\%$ and $2\%$: in these cases, an
important flattening of the epidemic curve can be observed for SA, where the
vaccination started early compared to the peak of infections. This situation
may be realized, as the vaccine is being administered to the population that
is more at risk of being infected by the virus. In the other extreme case,
represented by WNC, the vaccine is ineffective in changing the current wave
because the peak has already been attained before the vaccination campaign
started.

Our results confirm that the current vaccination strategy, which is performed
during a peak episode, is not effective to substantially slow down the spread
of the virus. On the other hand, the effectiveness for future waves is not a
question. It would be, in fact, very efficient to be able to administer the
vaccine to a larger portion of the population before the start of the next
wave.

###  III.4 Update of the vaccination to the first quarter 2021

As shown in the right column in Fig. 5, our simulation of the second wave,
done in mid December 2020, reproduces very well the epidemiological data up to
March 17, 2021. The only exception is Pacific, which has seen a sharper drop
in the number of new infections. Furthermore, one can clearly see a rebounce
in January that can be accredited to the Christmas holidays. Nevertheless,
this small effect does not have a significant impact on the agreement of our
prediction with the data.

In the first quarter of 2021, the vaccination campaign has also taken off
steadily, with nearly a quarter of the US population having received at least
one shot of vaccine. Furthermore, since February 27 the FDA has authorised the
use of the Janssen mono-dose vaccine 999https://www.fda.gov/emergency-
preparedness-and-response/coronavirus-disease-2019-covid-19/janssen-
covid-19-vaccine, which is now being administered together with the two-dose
Pfizer and Moderna vaccines. The data show that the rate of vaccinations has
been increasing approximately linearly with time, thus we updated the
prediction to take into account a vaccination fraction $c(t)$ growing linearly
with time:

$c(t)=u\;t\,,$ (7)

where the numerical value for each division are shown in Table 4.

The new results are shown in Fig. 7 for the 9 divisions. We consider both
people that received at least one shot (partial vaccination, in dash-dotted
lines) and fully vaccinated ones (dashed lines), with an offset of 4 weeks
between the beginning of the two. We consider them as two extreme cases,
defining a systematic error in our account of vaccinations. In all divisions,
the effect is minor, as the vaccination campaign has started too close to the
peak of the second wave. The only exception is Pacific, where taking into
account vaccinations substantially improves the agreement with the data.

The updated results confirm that a vaccination campaign operating during a
wave will not significantly affect the timing and height of the peak. Social
distancing and containment measures remain necessary. Conversely, vaccinating
a large portion of the population will certainly curb the eventual next wave.

| Percentage of the population vaccinated | Vaccination rate slope $u$
---|---|---
Division | Code | Partial vaccination | Full vaccination | Partial vaccination | Full vaccination
New England |  NE |  $30.1$ |  $16.4$ |  $0.00450$ |  $0.00641$
Mid-Atlantic |  MA |  $27.8$ |  $14.1$ |  $0.00415$ |  $0.00554$
South Atlantic |  SA |  $24.3$ |  $13.5$ |  $0.00363$ |  $0.00530$
East South Central |  ESC |  $23.1$ |  $12.6$ |  $0.00345$ |  $0.00492$
West South Central |  WSC |  $23.0$ |  $11.9$ |  $0.00344$ |  $0.00467$
East North Central |  ENC |  $25.6$ |  $14.7$ |  $0.00382$ |  $0.00575$
West North Central |  WNC |  $26.3$ |  $14.9$ |  $0.00393$ |  $0.00584$
Mountain |  M |  $25.7$ |  $14.8$ |  $0.00386$ |  $0.00579$
Pacific |  P |  $26.6$ |  $14.0$ |  $0.00397$ |  $0.00548$
Table 4: Percentage of the population vaccinated with at least one dose and
with two doses in each division as of the date of 24th of march 2021. Figure
7: Results of the eRG solutions for the second wave with a vaccination
campaign based on the data. Here, we consider a vaccination rate linearly
increasing in time, with slopes given in Table 4. The eRG parameters are the
same used for Fig. 5, based on data until December 28, 2020.

##  IV Discussion

In this paper we employ the epidemic Renormalization Group (eRG) framework in
order to understand, reproduce and predict the diffusion of the COVID-19
pandemic across the US as well as the effect of vaccination strategies. By
using flight data, we are able to see the changes in mobility across the
divisions, and observe how these changes affect the spread of the virus.
Furthermore, we show that the impact of the vaccination campaign on the
current wave of the pandemic in the US is marginal. Based on that, the
importance of social distancing is still relevant. Furthermore, we demonstrate
that the current wave is due to the endemic diffusion of the virus. Therefore,
building upon our previous results Cacciapaglia  _et al._ (2020a), in order to
control the next pandemic wave the number of daily new cases per million must
be around or less than 10-20 during the next inter-wave period. This
conclusion is further corroborated in Priesemann  _et al._ (2020) for Europe.

We learnt that the number of infected individuals in the current wave are not
affected measurably by the vaccination campaign. However, it is foreseeable
that it will impact specific compartments such as the overall number of
deceased individuals. Our study included an immunization rate between 0.64% to
2% of the total population each week. We also updated the results with the
actual rates of vaccination in the different divisions, as of March 24, 2020.
The results of our eRG model agree remarkably well with the new data from
December 28, 2020, to March 17, 2021. To curb the current and the next waves,
our results indisputably show that vaccinations alone are not enough and
strict social distancing measures are required until sufficient immunity is
achieved.

Our results should be employed by governments and decision makers to implement
local and global measures and, most importantly, the results of this paper can
be used as a foundation for vaccination campaign strategies for governments.

Given that pandemics are recurrent events, our results go beyond COVID-19 and
are universally applicable. What we have seen in the data for the US is that
it started in New York and, from there, it diffused to the rest of the
country. It is, therefore, important to contain future pandemics at an early
stage.

## References

  * Perc  _et al._ (2020) M. Perc, N. Gorišek Miksić, M. Slavinec, and A. Stožer, Frontiers in Physics 8, 127 (2020).
  * Hâncean  _et al._ (2020) M.-G. Hâncean, M. Perc, and L. Juergen, R. Soc. open sci. 7, 200780 (2020).
  * Zhou  _et al._ (2020) T. Zhou, Q. Liu, Z. Yang, J. Liao, K. Yang, W. Bai, and W. Zhang, J Evid. Based Med. 13, 3 (2020).
  * Prem  _et al._ (2020) K. Prem, Y. Liu, T. W. Russel, A. J. Kucharski, R. Eggo, and N. Davies, The Lancet Public Health 5, issue 5, E261 (2020).
  * Scala  _et al._ (2020) A. Scala, A. Flori, A. Spelta, E. Brugnoli, M. Cinelli, W. Quattrociocchi, and F. Pammolli, Sci Rep 10, 13764 (2020).
  * Friston  _et al._ (2020) K. J. Friston, T. Parr, P. Zeidman, A. Razi, G. Flandin, J. Daunizeau, O. J. Hulme, A. J. Billig, V. Litvak, C. J. Price, R. J. Moran, and C. Lambert, “Second waves, social distancing, and the spread of covid-19 across america,” (2020), arXiv:2004.13017 [q-bio.PE] .
  * Sonnino  _et al._ (2020) G. Sonnino, F. Mora, and P. Nardone, “A stochastic compartmental model for covid-19,” (2020), arXiv:2012.01869 [physics.med-ph] .
  * Abou-Ismail (2020) A. Abou-Ismail, SN Compr Clin Med. , 1 (2020).
  * Kermack  _et al._ (1927) W. O. Kermack, A. McKendrick, and G. T. Walker, Proceedings of the Royal Society A 115, 700 (1927).
  * Zhan  _et al._ (2018) X.-X. Zhan, C. Liu, G. Zhou, Z.-K. Zhang, G.-Q. Sun, J. J. Zhu, and Z. Jin, Applied Mathematics and Computation 332, 437 (2018).
  * Perc  _et al._ (2017) M. Perc, J. J. Jordan, D. G. Rand, Z. Wang, S. Boccaletti, and A. Szolnoki, Physics Reports 687, 1 (2017).
  * Wang  _et al._ (2015) Z. Wang, M. A. Andrews, Z.-X. Wu, L. Wang, and C. T. Bauch, Physics of Life Reviews 15, 1 (2015).
  * Scudellari (2020) M. Scudellari, Nature 584, 22 (2020).
  * Della Morte  _et al._ (2020) M. Della Morte, D. Orlando, and F. Sannino, Front. in Phys. 8, 144 (2020).
  * Cacciapaglia and Sannino (2020a) G. Cacciapaglia and F. Sannino, Sci Rep 10, 15828 (2020a), arXiv:2005.04956 [physics.soc-ph] .
  * Della Morte and Sannino (2020) M. Della Morte and F. Sannino, “Renormalisation Group approach to pandemics as a time-dependent SIR model,” (2020), arXiv:2007.11296 [physics.soc-ph] .
  * Cacciapaglia and Sannino (2020b) G. Cacciapaglia and F. Sannino, “Evidence for complex fixed points in pandemic data,”  (2020b), arXiv:2009.08861 [physics.soc-ph] .
  * Cacciapaglia  _et al._ (2020a) G. Cacciapaglia, C. Cot, and F. Sannino, “Multiwave pandemic dynamics explained: How to tame the next wave of infectious diseases,” (2020a), arXiv:2011.12846 [physics.soc-ph] .
  * Taubenberger and Morens (2006) J. K. Taubenberger and D. M. Morens, Rev Biomed 17(1), 69 (2006).
  * Wilson (1971a) K. G. Wilson, Phys. Rev. B 4, 3174 (1971a).
  * Wilson (1971b) K. G. Wilson, Phys. Rev. B 4, 3184 (1971b).
  * Li  _et al._ (2019) L. Li, J. Zhang, C. Liu, H.-T. Zhang, Y. Wang, and Z. Wang, Applied Mathematics and Computation 347, 566 (2019).
  * Wang  _et al._ (2016) Z. Wang, C. T. Bauch, S. Bhattacharyya, A. d’Onofrio, P. Manfredi, M. Perc, N. Perra, M. Salathé, and D. Zhao, Physics Reports 664, 1 (2016).
  * Danby (1985) J. M. A. Danby,  _Computing applications to differential equations modelling in the physical and social sciences_ (Reston Publishing Company, Reston VA (USA), 1985).
  * Brauer (2019) F. Brauer, Journal of Biological Dynamics 13, 23 (2019).
  * Miller (2012) J. C. Miller, Bulletin of mathematical biology 74, 2125 (2012).
  * Murray (2002) J. D. Murray,  _Mathematical biology_ , 3rd ed., Interdisciplinary applied mathematics (Springer, New York (USA), 2002).
  * Fishman  _et al._ (2014) D. Fishman, E. Khoo, and A. Tuite, PLOS Currents Outbreaks 6 (2014).
  * Pell  _et al._ (2018) B. Pell, Y. Kuang, C. Viboud, and G. Chowell, Epidemics 22, 62 (2018), the RAPIDD Ebola Forecasting Challenge.
  * Paltiel  _et al._ (2020) A. D. Paltiel, J. L. Schwartz, A. Zheng, and R. P. Walensky, Health Affairs (2020).
  * Schäfer  _et al._ (2014) M. Schäfer, M. Strohmeier, V. Lenders, I. Martinovic, and M. Wilhelm, in  _IPSN-14 Proceedings of the 13th International Symposium on Information Processing in Sensor Networks_ (IEEE, 2014) pp. 83–94.
  * Islind  _et al._ (2020) A. S. Islind, M. Óskarsdóttir, and H. Steingrímsdóttir, “Changes in mobility patterns in europe during the covid-19 pandemic: Novel insights using open source data,” (2020), arXiv:2008.10505 [cs.CY] .
  * Bank of England (2020) Bank of England, “Monetary policy report and interim financial stability report - may 2020,” https://www.bankofengland.co.uk/report/2020/monetary-policy-report-financial-stability-report-may-2020 (2020).
  * Cacciapaglia  _et al._ (2020b) G. Cacciapaglia, C. Cot, and F. Sannino, Sci Rep 10, 15514 (2020b), arXiv:2007.13100 [physics.soc-ph] .
  * Cardy and Grassberger (1985) J. L. Cardy and P. Grassberger, Journal of Physics A: Mathematical and General 18, L267 (1985).
  * Yang  _et al._ (2020) C. Yang, D. Sha, Q. Liu, Y. Li, H. Lan, W. W. Guan, T. Hu, Z. Li, Z. Zhang, J. H. Thompson, and et al., International Journal of Digital Earth 13, 1186–1211 (2020).
  * Priesemann  _et al._ (2020) V. Priesemann, M. M. Brinkmann, S. Ciesek, S. Cuschieri, T. Czypionka, G. Giordano, D. Gurdasani, C. Hanson, N. Hens, E. Iftekhar, _et al._ , The Lancet (2020).

## Author contribution

This work has been designed and performed conjointly and equally by the
authors, who have equally contributed to the writing of the article. ASI and
MO have extracted and processed the data from flights; CC has worked on the
numerical results from the eRG equations and analysed the epidemiological
data.

## Competing interests

The authors declare no competing interests.


We investigate algebraic and compositional properties of abstract multiway rewriting systems, which are archetypical structures underlying the formalism of the Wolfram model. We demonstrate the existence of higher homotopies in this class of rewriting systems, where homotopical maps are induced by the inclusion of appropriate rewriting rules taken from an abstract rulial space of all possible such rules. Furthermore, we show that a multiway rewriting system with homotopies up to order $n$ may naturally be formalized as an $n$-fold category, such that (upon inclusion of appropriate inverse morphisms via invertible rewriting relations) the infinite limit of this structure yields an ${\infty}$-groupoid. Via Grothendieck's homotopy hypothesis, this ${\infty}$-groupoid thus inherits the structure of a formal homotopy space. We conclude with some comments on how this computational framework of homotopical multiway systems may potentially be used for making formal connections to homotopy spaces upon which models relevant to  physics may be instantiated.

§ INTRODUCTION

The Wolfram model<cit.> <cit.>  (see also <cit.> <cit.> <cit.>  <cit.>) is an explicitly computational framework based upon abstract rewriting that attempts to formalize the fundamental mathematical structures underlying space, time and matter. The model seeks to capture the ways in which simple abstract rewriting rules may be composed, so as to yield more complex structures (suggestively titled “universes”) that admit certain emergent “laws of physics”. This general approach of seeking a constructivist formalization of certain fundamental structures has strong parallels in recent work in the foundations of both physics and mathematics, for instance in relation to synthetic geometry and cohesive homotopy type theory<cit.><cit.><cit.><cit.>. The general framework of the Wolfram model eschews the rigid continuum description of spacetime in terms of Lorentzian manifolds in favor of a more rudimentary description of the combinatorics and topology of the intrinsic causal relationships between events, which, in turn, makes manifest some of the computational architecture that underlies contemporary physics. The archetypical structures that appear within this framework are so called multiway systems - non-deterministic abstract rewriting systems equipped with a notion of causal structure. The “multiway” moniker here designates the fact that all permissible applications of the rewriting rules are instantiated in all possible orderings, leading to multiple chains of rewrite terms that are partially-ordered by causality. Depending upon the precise interpretation of these terms (or states), Wolfram model multiway systems may be realized as rewriting systems over hypergraphs, character strings, ZX-diagrams<cit.><cit.> (as formalized by Coecke and Duncan<cit.>), or  string diagrams<cit.>.

In this paper, we consider abstract multiway rewriting systems, i.e. systems in which the terms/states are purely abstract entities, since doing so will be more convenient for studying the algebraic and compositional properties of such systems at the highest possible level of generality; these general results may then be easily specialized or extrapolated to particular classes of rewriting systems. Specifically, our analysis reveals a rich homotopical structure admissible on these abstract multiway rewriting systems, where these homotopies are induced by the inclusion of certain “higher-order” rules taken from an abstract rulial space, designating the space of all possible rewriting rules of a given signature (formally, at least for the case of hypergraphs, a monoidal category of cospans). Although a generic abstract rewriting system may be thought of as simply being an $F$-coalgebra for the power set functor  ${\mathcal{P}}$  (which, at least in the case of hypergraph rewriting, is further equipped with a natural symmetric monoidal structure), we show here that a multiway rewriting system equipped with homotopies up to order $n$ may elegantly be formalized as an $n$-fold category. Furthermore, upon including inverse morphisms (via the inclusion of invertible rewriting relations), the infinite limit of this structure naturally yields an ${\infty}$-groupoid, which therefore inherits the structure of a formal homotopy space via Grothendieck's homotopy hypothesis <cit.>.  For specific models of ${\infty}$-categories, these structures correspond precisely to topological spaces, hence allowing one to make formal connections between the abstract computational framework of Wolfram model multiway systems and the topological spaces (such as Lorentzian manifolds, or projective Hilbert spaces) upon which laws of physics are traditionally instantiated.

§ PRELIMINARIES OF WOLFRAM MODEL MULTIWAY SYSTEMS

This section provides an overview of the basic abstract rewriting constructions that constitute the Wolfram model, with a particular emphasis on their non-deterministic aspects, as captured via multiway rewriting systems. We present both combinatorial and compositional descriptions of such systems, as well as an abstract description of their rulial space (i.e. the space of all possible rewriting rules) for the case of (hyper)graph rewriting, by means of a double-pushout (DPO) rewriting formalism over (selective) adhesive categories. The contents of this section serve merely as a preliminary introduction to the known constructions of Wolfram model multiway rewriting systems; the informed reader is encouraged to proceed directly to Section <ref>, in which the main arguments of this paper are presented.

§.§ Combinatorial & Compositional Descriptions of Multiway  Systems

Combinatorially, a multiway rewriting system is simply a directed, acyclic graph of states, determined by abstract rewriting rules that inductively generate a (potentially infinite) multiway evolution graph, together with a partial order on the rewrite rule applications, determined by their causal structure. This is straightforward to formalize in terms of abstract rewriting systems <cit.> in which the underlying rewrite relation is not (necessarily) confluent<cit.>.

An “abstract rewriting system” (or “ARS”) is a set, denoted $A$ (whose elements are known as “objects”), equipped with a binary relation, denoted ${\to}$ (known as the “rewrite relation”).

A “multiway evolution graph”, denoted ${G_{multiway}}$, is a directed, acyclic graph corresponding to the evolution of a (generically non-confluent) abstract rewriting system ${\left( A, \to \right)}$, in which the set of vertices corresponds to the set of objects ${V \left( G_{multiway} \right)}$, and in which the directed edge ${A \to B}$ exists in ${E \left( G_{multiway} \right)}$ if and only if there exists an application of the rewrite relation that transforms object $A$ to object $B$.

Hence, directed edges will connect vertices $A$ and $B$ in ${G_{multiway}}$ if and only if ${A \to B}$ in the underlying rewriting system, and a directed path will connect vertices $A$ and $B$ if and only if ${A \to^{*} B}$, where ${\to^{*}}$ denotes the reflexive transitive closure of ${\to}$, i.e if and only if there exists a finite rewrite sequence of the form:

\begin{equation}
A \to A^{\prime} \to A^{\prime \prime} \to \cdots \to B^{\prime} \to B.
\end{equation}

An example of a multiway evolution graph for a simple Wolfram model evolution (i.e. a non-deterministic hypergraph rewriting system) is shown in Figure <ref>.

The multiway evolution graph corresponding to the first 3 steps in the non-deterministic evolution history of the hypergraph substitution rule ${\left\lbrace \left\lbrace x, y \right\rbrace, \left\lbrace y, z \right\rbrace \right\rbrace \to \left\lbrace \left\lbrace w, y \right\rbrace, \left\lbrace y, z \right\rbrace, \left\lbrace z, w \right\rbrace, \left\lbrace x, w \right\rbrace \right\rbrace}$, starting from a simple double self-loop initial condition ${\left\lbrace \left\lbrace 0, 0 \right\rbrace, \left\lbrace 0, 0 \right\rbrace \right\rbrace}$.

In the most general case, in which the rewrite relation ${\to}$ is treated as an indexed union of sub-relations, i.e. ${\to = \to_1 \cup \to_2 \cup \dots}$, with label set ${\Lambda}$ for the indices, and in which the set of objects $A$ is arbitrary (i.e. its elements could represent graphs, hypergraphs, string diagrams, terms, character strings, etc.), the resulting labeled abstract rewriting system ${\left( A, \Lambda, \to \right)}$ permits an elegant compositional description in terms of $F$-coalgebras. Specifically, note that the system ${\left( A, \Lambda, \to \right)}$ is now simply a bijective function from $A$ to a subset of the power set of $A$ indexed by ${\Lambda}$, i.e. ${\mathcal{P} \left( \Lambda \times A \right)}$:

\begin{equation}
p \mapsto \left\lbrace \left( \alpha, q \right) \in \Lambda \times A : p \to^{\alpha} q \right\rbrace.
\end{equation}

Recall now that an $F$-coalgebra for an endofunctor ${F : \mathbf{C} \to \mathbf{C}}$ consists of an object $A$ in ${\mathrm{ob} \left( \mathbf{C} \right)}$ equipped with a morphism ${\alpha : A \to F A}$ in ${\mathrm{hom} \left( \mathbf{C} \right)}$, hence denoted ${\left( A, \alpha \right)}$. Thus, since the power set construction on ${\mathbf{Set}}$ is a covariant endofunctor ${\mathcal{P} : \mathbf{Set} \to \mathbf{Set}}$, we see that the abstract rewriting system ${\left( A, \to \right)}$ consists of an object $A$ equipped with an additional morphism of ${\mathbf{Set}}$, namely the rewrite relation ${\to}$:

\begin{equation}
\to : A \to \mathcal{P} A.
\end{equation}

The causal structure of the abstract rewriting system also endows the multiway system with an additional partial order relation on the rewrite applications themselves<cit.>, which can in turn be described compositionally in terms of a partial monoidal structure<cit.><cit.><cit.>. However, for the purposes of the present paper, the causal structure can (for the most part) be safely ignored.

§.§ Rulial Space and Double-Pushout (DPO) Rewriting

In a large class of cases, such as (hyper)graph rewriting, open graph rewriting, string diagram rewriting (including the ZX-calculus), etc., the rewrite relation ${\to}$ may be further specified as a double-pushout (DPO) rewriting system; for instance, this was formalized by Dixon and Kissinger<cit.> for the case of open graphs using a selective adhesive category of open spans, and the same formalism was subsequently applied to both ZX-diagrams and arbitrary hypergraph rewriting systems as special cases<cit.><cit.>. Specifically, the rewriting rules are defined as spans of monomorphisms ${\rho}$ of the form:

\begin{equation}
\rho = \left( l : K \to L, r : K \to R \right),
\end{equation}

with the left- and right-hand-sides of the rule specified by objects $L$ and $R$, respectively, and with object $K$ designating the interface graph. A match for the rule ${\rho}$ within an object $G$ is simply a morphism ${m : L \to G}$, with the rule being applicable if and only if there exists a pair of pushout diagrams of the form <cit.>:

\begin{equation}
\begin{tikzcd}
L \arrow[d, "m"] & K \arrow[l, "l"] \arrow[d, "k"] \arrow[r, "r"] & R \arrow[d, "n"]\\
G & D \arrow[l, "f"] \arrow[r, "g"] & H
\end{tikzcd}.
\end{equation}

These rules are defined and applied in the context of an adhesive category<cit.>, in which every pushout along a monomorphism satisfies the van-Kampen square condition, such that, for every commutative diagram of the form:

\begin{equation}
\begin{tikzcd}
B^{\prime} \arrow[ddd, "f_{h}^{\prime}"] \arrow[dr, "h_B"] & & & A^{\prime} \arrow[lll, "g_h"] \arrow[dl, "h_A"] \arrow[ddd, "f_h"]\\
& B \arrow[d, "f^{\prime}"] & A \arrow[l, "g"] \arrow[d, "f"] &\\
& D & C \arrow[l, "g^{\prime}"] &\\
D^{\prime} \arrow[ur, "h_D"] & & & C^{\prime} \arrow[lll, "g_{h}^{\prime}"] \arrow[ul, "h_C"]
\end{tikzcd},
\end{equation}

such the following subdiagrams are both pullbacks:

\begin{equation}
\begin{tikzcd}
B^{\prime} \arrow[d, "h_B"] & A^{\prime} \arrow[l, "g_h"] \arrow[d, "h_A"]\\
B & A \arrow[l ,"g"]
\end{tikzcd}, \qquad \text{ and } \qquad
\begin{tikzcd}
A \arrow[d, "f"] & A^{\prime} \arrow[l, "h_A"] \arrow[d, "f_h"]\\
C & C^{\prime} \arrow[l, "h_C"]
\end{tikzcd},
\end{equation}

the pushouts and pullbacks satisfy a compatibility condition. In fact, selective adhesive categories<cit.> allow DPO rewriting to be defined in an even broader class of cases, in which one is working within a full subcategory ${\mathbf{C}^{\prime}}$ of an adhesive category ${\mathbf{C}}$, with embedding functor ${S : \mathbf{C}^{\prime} \to \mathbf{C}}$, with the only condition being that $S$ preserves monomorphisms.

This construction allows us to formalize the notion of rulial space, i.e. the space of all possible rewriting rules of a given signature, at least for the case of Wolfram model hypergraph rewriting systems (and systems such as string diagram rewriting systems that are known to be special cases):

The “rulial space” of Wolfram model systems is the category of cospans of a DPO rewriting system, defined over the selective adhesive category of Wolfram model hypergraphs.

The rulial space of Wolfram model systems functorially acquires the structure of a selective adhesive category. Indeed, as a consequence of the concurrency and parallelism theorems of algebraic graph transformation theory<cit.>, the rulial space also inherits a natural monoidal structure<cit.> (which is, in turn, inherited by all Wolfram model multiway systems, since they are obtained by the composition of certain rules in rulial space).

Even though the rewriting formulation above borrows from DPO rewriting as per  <cit.>, it is worth noting other recent developments in rewriting theory, that  will no doubt be useful for our model in future work. We point the reader to the following relevant works: (i) mathematical foundations of  string diagram rewriting, reported in  <cit.>,  <cit.>; (ii) higher dimensional rewriting theory <cit.>,  <cit.>; (iii) algebraic graph transformations and compositional systems  <cit.>,  <cit.>.

§ HOMOTOPICAL MULTIWAY SYSTEMS AS $N$-FOLD CATEGORIES

Given the compositional setup outlined above, we now address explicit constructions  of higher homotopies upon abstract multiway systems, and introduce a systematic algorithm for identifying rewriting rules resulting in these homotopies. Subsequently, we prove that a multiway system equipped with homotopies up to order $n$ may be formalized as an $n$-fold category, such that the infinite limit of such a higher-order multiway system yields an ${\infty}$-groupoid (upon the admission of invertible rewriting rules).

§.§ A Computational Framework for Constructing Homotopical Multiway  Systems

To avoid us needing to deal with the added complexities of full hypergraph rewriting systems, let us consider instead a very simple rewriting system on character strings, as specified by the rule ${A \to AB}$, generating the multiway evolution graph shown in Figure <ref>. We have deliberately selected a very minimal example, although the basic algorithm described here applies in principle to all classes of multiway systems, as discussed below. Every path in this multiway evolution graph corresponds to a proof of a proposition, as shown on the left of Figure <ref> for the particular case of a proof of the proposition that ${AA \to ABBBABBB}$, subject to the axiom ${A \to AB}$. A key feature of multiway systems is that one can consider multiple paths connecting the same pair of vertices, which may be interpreted type-theoretically as corresponding to the existence of multiple proofs of the same proposition, as shown on the right of Figure <ref>.

The multiway evolution graph corresponding to the first 8 steps in the non-deterministic evolution history of the string substitution rule ${A \to AB}$, starting from the two-character initial condition $AA$.

On the left, a multiway evolution graph with a red highlighted path between vertices ${AA}$ and ${ABBBABBB}$, corresponding to a proof of the proposition ${AA \to ABBBABBB}$, subject to the string substitution axiom ${A \to AB}$. On the right, the same multiway evolution graph showing multiple highlighted paths (red and yellow) between vertices ${AA}$ and ${ABBBABBB}$, illustrating the existence of multiple proofs of the proposition ${AA \to ABBBABBB}$.

A key feature of homotopy type theory is the ability to describe proofs of equivalences between proofs in terms of homotopies between paths<cit.>. In the discrete setup of a multiway evolution graph, a concrete realization of the homotopy map between paths can be given by simply introducing bi-directional edges mapping vertices on one path to corresponding vertices on the other path, and back again, as shown in Figure <ref>. However, this approach to introducing homotopy maps is somewhat ad hoc; it is far more natural to introduce these bi-directional relationships between multiway vertices as a new set of rewrite relations that can be appended to the original multiway system rules, as shown in Figure <ref>. For the particular system at hand, these additional rules are:

\begin{multline}
\hspace{2cm}  \left\lbrace AAB \to ABA, \,  AABB \to ABBA, \, AABBB \to ABBBA, \right.\\
\left. ABABBB \to ABBBAB, \,  ABBABBB \to ABBBABB \right\rbrace  \hspace{3cm}
\end{multline}

and their respective inverses (which have significance for the forthcoming constructions of groupoids):

\begin{multline}
\hspace{2cm} \left\lbrace ABA \to AAB, \, ABBA \to AABB, \, ABBBA \to AABBB, \right.\\
\left. ABBBAB \to ABABBB, \, ABBBABB \to ABBABBB \right\rbrace \hspace{3cm}
\end{multline}

A multiway evolution graph with purple highlighted paths between vertices on one (red) path and corresponding vertices on the other (yellow) path, and back again, interpreted as homotopic maps between the two associated proofs of the proposition ${AA \to ABBBABBB}$.

A multiway evolution graph with explicit rules added to map vertices on one (red) path to corresponding vertices on the other (yellow) path, and back again, interpreted as homotopic maps between the two associated proofs of the proposition ${AA \to ABBBABBB}$, shown in purple.

This general algorithmic procedure for introducing homotopy maps between paths in an arbitrary multiway system may now be iterated, so as to introduce higher homotopy maps between homotopies; for instance, order-3 homotopies (3-cells) between order-2 homotopies (2-cells) and so on (more of this will be reported in forthcoming work  <cit.>).

§.§ Higher Categorical Formulation of Homotopical Multiway Systems

We now attempt to formalize the above constructions in the language of higher category theory. Starting with the order-2 homotopies (2-cells or 2-morphisms) between paths (1-morphisms) we make the following proposition:

The multiway rewriting system shown in Figure <ref> is an order-2 homotopy rewriting system with 2-morphisms between paths, thus yielding a double category.

Let us first recall the standard definition of a double category.

A double category $D$, denoted    D_1 [r, shift left]   [r, shift right]      D_0     is defined in terms of the following conditions:

* The objects of $D$ are the objects of $D_0$

* $D$ has vertical morphisms, which are the morphisms of ${D_0}$

* Additionally, $D$ also has horizontal morphisms, which are the objects of ${D_1}$

* Finally, the 2-morphisms of $D$ (also referred to as squares or 2-cells) are the morphisms of ${D_1}$

Here, a 2-cell in $D$ can be represented by the following commutative square:

\begin{equation}
\begin{tikzcd}[column sep={2cm,between origins},row sep={2cm,between origins}]
a \dar{l} \rar["f",""{name=U, below}""] & c \dar{m} \\
b \rar["g" below,""{name=D}""] & d \\
\arrow[Rightarrow, from=U, to=D,"\phi"]
\end{tikzcd}
\end{equation}

where $a$, $b$, $c$, $d$ are objects; $l$, $m$ are vertical arrows; $f$, $g$ are horizontal arrows; and $\phi$ denotes the 2-cell.

Note that a double category is an internal category in Cat. Vertical composition in a double category is given by composition in the categories $D_0$ and $D_1$, while horizontal composition is given by composition on   D_1 [r, shift left]   [r, shift right]      D_0    itself, by virtue of it being a category internal to Cat. This set-up will be useful for proving the desired proposition.

Proof of Proposition <ref>:
Following the functorial description provided in Section <ref>, any abstract  multiway system may be thought of as a category ${\cal M}_0$ whose  objects are given by states (rewrite terms) of the multiway system, and whose paths of causally-ordered rewrite chains refer to morphisms of this category. Composition is given by concatenation of directed paths. Transitive closure along directed paths ensures associativity and self-loops on nodes (trivially obtained via identity rewrite rules) ensure identity laws of a category.

From the definition of a double category given above, ${\cal M}_0$ is precisely the category $D_0$, with objects being the vertices of a multiway evolution graph (a prototypical example being given in Figure <ref>),  and paths formed by multiway evolution (also shown in Figure <ref>) corresponding to vertical morphisms in $D_0$.
The horizontal arrows in Figure <ref>, colored in purple, are induced by the introduction of additional rules to the rewriting system and serve the purpose of connecting (vertical) paths in the multiway system. Hence, these can be identified as objects of another category, which we denote as ${\cal M}_1$ (${\cal M}_1$ will be identified with $D_1$ once we specify its morphisms).
Imposing commutativity over squares formed by the vertical and horizontal arrows of the multiway system gives us the morphisms of ${\cal M}_1$; these squares, exemplified in Figure <ref>, are 2-morphisms of the double category     M_1 [r, shift left]   [r, shift right]      M_0   .

Given this double category, the 2-morphism in the multiway evolution graph between the highlighted red and yellow paths in Figure <ref> is precisely given by vertical compositions of squares. These are the 2-cells of the new multiway rewriting system, enhanced with additional rules, thus making it a 2-homotopy rewriting system. (Note that the triangles at the start and end points of the 1-paths are also squares with an identity morphism at those points.)

This leads us directly to the following corollary:

A 2-homotopy multiway system extends to a double groupoid upon admitting invertible rules to those used in the multiway construction (including the homotopy rules).

Proof of Corollary <ref>:
Adding invertible rules to the multiway system equips every rewrite arrow with an inverse arrow, such that their composition is the identity map on either rewrite term. Given an order-2 homotopy  multiway system, the presence of invertible arrows implies that all 1-morphisms and 2-morphisms are isomorphisms. This extends the multiway system in Proposition <ref> to a double groupoid.

Furthermore, as shown in the previous section, with appropriate additional rules, one can continue the above construction inductively to obtain rewriting systems with progressively higher homotopies between paths. This yields the following proposition:

So long as additional rewrite rules between cells, up to order $n - 1$, are admissible, then the multiway system shown in Figure <ref> can be enhanced to yield an $n$-homotopy rewriting system with $n$-morphisms between paths, thus yielding an $n$-fold category.

Proof of Proposition <ref>:
The proof of this proposition follows by a simple inductive construction.

A double category itself is a  special case of an $n$-fold category. In Proposition <ref>, we proved this proposition for the $n = 2$ case.

For the $n = 3$ case, we need to construct an order-3 homotopy rewriting system equipped with 3-morphisms, which are cubes bounded by squares of the given double category. This can be done by adding specific rewrite rules that give arrows between paths (1-cells) in such a way that one obtains multiple 2-cells, which eventually permit additional arrows between them, giving 3-morphisms as cubes bounded by squares. These cubes then define a 3-fold category as follows:

A 3-fold category, denoted     D_2 [r, shift left]   [r, shift right]      D_1 [r, shift left]   [r, shift right]      D_0      is defined in terms of the following conditions:

* The objects are the objects of $D_0$

* The vertical arrows are the morphisms of $D_0$

* The horizontal arrows are the objects of $D_1$

* The vertical squares are the morphisms of $D_1$

* The horizontal squares are the objects of $D_2$

* The cube bounded by vertical and horizontal squares is a morphism of $D_2$

The above data has to satisfy the commutative diagram:

\begin{equation}
\begin{tikzcd}[row sep=1.5em]
A \arrow[rr,"f",""{name=0, below}] \arrow[dr,swap,"a",""{name=1}""] \arrow[dd,swap,"h",""{name=2}""] &&
B \arrow[dd,swap,"h'" near end,""{name=3}""] \arrow[dr,"b",""{name=4, below}""] \\
& A'\arrow[rr, crossing over, "f'" left,""{name=U, below}"", yshift=-0.05ex]
&& B' \arrow[dd,"k'",""{name=5}""] \\
C \arrow[rr,"g" near start,""{name=N, below}""] \arrow[dr,swap,"c",""{name=6}"", below] && D \arrow[dr,"d",""{name=7, below}""] \\
& C' \arrow[rr,"g'" near start,""{name=M}""] \arrow[uu,<-,crossing over,"k" near end]&& D'
\arrow[from=1,to=4, dash, shorten=10mm, phantom,""{name=K}""]
\arrow[r,from=6,to=7, dash, shorten=10mm, phantom,""{name=L}""]
\arrow[Rightarrow, line width=0.3pt, from=K, to=L,"\psi"] \\
\arrow[r, dash, line width=0.3pt, from=K, to=L]
\arrow[Rightarrow, from=0, to=N, xshift=2ex, crossing over, shorten=1mm, "\eta", near start, swap]
\arrow[Rightarrow, from=U, to=M, xshift=2ex, crossing over, shorten=1mm, "\omega", near end, swap]
\arrow[Rightarrow, from=1, to=6, shorten=2mm, "\alpha",near start]
\arrow[Rightarrow, from=4, to=7, shorten=1mm, "\beta"]
\arrow[r, from=1, to=4, shorten =11mm, xshift=10ex, yshift=-3.03ex, crossing over]
\arrow[r, dash, from=1, to=4, shorten=11mm, xshift=-0.3ex, yshift=-3.03ex, crossing over]
\end{tikzcd}
\end{equation}

where arrows are 1-morphisms, squares are 2-morphisms and cubes are 3-morphisms (3-cells) of the 3-fold category.   Compositions, identity and associative laws on each of these entities follow from their usual definitions in $D_0$, $D_1$ and $D_2$ where appropriately applicable.

The 3-cells in a multiway system can then be obtained by vertical, horizontal and sideways compositions of cubes of a 3-fold category.

Likewise, one can continue this process of constructing higher morphisms indefinitely, so long as additional rewrite rules between cells of up to order $n - 1$ are admissible such that $n$-morphisms are now $n$-hypercubes defined within an $n$-fold category, and $n$-cells of such multiway systems are simply compositions of $n$-hypercubes in $n$-directions. The ensuing  $n$-fold category can then be expressed as:

\[  \begin{tikzcd} {\cal M}_{n - 1} \arrow[r, shift left] \arrow[r, shift right]  & {\cal M}_{n - 2} \,\,  \cdots\cdots \,\, {\cal M}_2 \arrow[r, shift left] \arrow[r, shift right] & {\cal M}_1 \arrow[r, shift left] \arrow[r, shift right]  & {\cal M}_0 \end{tikzcd} \]

where the objects of ${\cal M}_i$ (for $0 < i < {n-2}$) are $i$-dimensional hypercubes in ${\mathbb R}^n$ whose normal vectors are oriented along the $z$-axis in ${\mathbb R}^n$. The morphisms of ${\cal M}_i$ are $(i + 1)$-dimensional hypercubes in ${\mathbb R}^n$ whose normal vectors are oriented orthogonal to the $z$-axis in ${\mathbb R}^n$. Furthermore, the objects of ${\cal M}_{n - 1}$  are $(n - 1)$-dimensional hypercubes with normal vectors oriented along the $z$-axis in ${\mathbb R}^n$ and the morphism in ${\cal M}_{n - 1}$ is the commutative $n$-hypercube composed of the above.

This iterative definition of an $n$-fold category corroborates the previously described construction, in which homotopies of up to order $n$ in multiway systems are realized by including additional rewrite rules that introduce new arrows to obtain squares (and so on), until one eventually obtains a commutative $n$-hypercube.

The benefit gained by realizing $n$-fold categories in our constructions is the ease of expressing multiple compositions by gluing hypercubes up to order $n$. Composable arrays of hypercubes can then be used to construct all ($i \leq n$)-cells of the category.

For example, in a 2-fold category, we can define a composable array of 2-dimensional elements (squares) to be such that any array is composable with its immediate neighbors (the associative and commutative laws imply that the composition is well defined). This process easily extends in $n$-fold categories  using $n$-dimensional elements ($n$-hypercubes).

The following remark is in order:

As in Corollary <ref>, a $n$-homotopy multiway system extends to an $n$-fold groupoid upon admitting invertible rules that ensure invertibility of all higher morphisms.

Given the iterative construction above, where one induces higher homotopies in multiway systems via the inclusion of supplementary rules, one can now ask what the $n \to \infty$ limit of this construction yields?

The $n \to \infty$ limit of a  multiway rewriting system with invertible homotopy rules, corresponding to higher morphisms (so long as they are admissible), is an $\infty$-groupoid.

Proof of Proposition <ref>:
Following Proposition <ref> and Remark <ref>, multiway rewriting systems equipped with homotopies up to order $n$, along with invertible rules, yield an $n$-fold groupoid.
The $n \to \infty$ limit of an $n$-fold groupoid is precisely an $\infty$-groupoid. This limit exists so long as the infinite hierarchy of rules required to iteratively construct higher morphisms is admissible on that multiway rewriting system.

§ IMPLICATIONS AND OUTLOOK

This paper has illustrated how higher categorical structures arise naturally within the abstract non-deterministic rewriting systems of the Wolfram model, at least when formalized as multiway systems.  We have shown how higher homotopies induced on multiway systems via specific rewriting rules correspond to morphisms of an $n$-fold category, and we have both demonstrated an explicit computational procedure for constructing these homotopies on multiway rewriting systems, as well as provided the formal correspondence of these systems to $n$-fold categories and $n$-fold groupoids. Interestingly, the $n \to \infty$ limit of this structure yields an $\infty$-groupoid, with the latter being relevant from the point of view of Grothendieck's homotopy hypothesis <cit.>; we speculate that this correspondence might potentially offer a new way to understand how the kinds of spatial structures relevant to the foundations of physics might potentially arise from purely combinatorial rewriting systems <cit.><cit.><cit.>. Our results here provide a preliminary framework for investigating such connections to topological and geometric spaces from an underlying computational perspective, which we outline in forthcoming work <cit.>.

On the other hand, complementary to what we have shown here, there is an ongoing research program involving the application of methods from homotopy type theory and infinity-category theory, that seeks to formalize, among other things, the mathematical foundations of synthetic spatial structures (in particular, cohesive topological and geometric spaces) <cit.><cit.><cit.>. Constructions of synthetic geometry are also of great relevance for the foundations of physics. In particular, there are active efforts in the formalization of quantum field theories on a cohesive $\infty$-topos   <cit.>. In a sense, the Wolfram model can be thought of as a “constructivist” realization of some of the above ideas, expressed via a concrete model of computation, namely abstract rewriting systems. The precise connections of Wolfram model multiway rewriting  systems to homotopy types and infinity-categories is also a topic for future work.

§.§.§ Acknowledgments
The authors would like to thank Stephen Wolfram for his encouragement, and for the many useful conversations and suggestions that helped shape this work.