#  LYCEN 2021-01
XXX

Giacomo<EMAIL_ADDRESS>Corentin
<EMAIL_ADDRESS>Adele de<EMAIL_ADDRESS>
Stefan<EMAIL_ADDRESS>Francesco
<EMAIL_ADDRESS>and Shahram
<EMAIL_ADDRESS>
#  LYCEN 2021-01
Epidemiological theory of virus variants

Giacomo<EMAIL_ADDRESS>Corentin
<EMAIL_ADDRESS>Adele de<EMAIL_ADDRESS>
Stefan<EMAIL_ADDRESS>Francesco
<EMAIL_ADDRESS>and Shahram
<EMAIL_ADDRESS>
1 Institut de Physique des 2 Infinis (IP2I) de Lyon, CNRS/IN2P3, UMR5822,

69622 Villeurbanne, France

2 Université de Lyon, Université Claude Bernard Lyon 1, 69001 Lyon, France

3 Politecnico di Torino, Torino, Italy

4 Scuola Superiore Meridionale, Largo S. Marcellino, 10, 80138 Napoli NA,
Italy

5 Dipartimento di Fisica, E. Pancini, Univ. di Napoli, Federico II and INFN
sezione di Napoli,

Complesso Universitario di Monte S. Angelo Edificio 6, via Cintia, 80126
Napoli, Italy

6 CP3-Origins and D-IAS, Univ. of Southern Denmark, Campusvej 55, DK-5230
Odense, Denmark

Abstract:

We propose a physical theory underlying the temporal evolution of competing
virus variants that relies on the existence of (quasi) fixed points capturing
the large time scale invariance of the dynamics. To motivate our result we
first modify the time-honoured compartmental models of the SIR type to account
for the existence of competing variants and then show how their evolution can
be naturally re-phrased in terms of flow equations ending at quasi fixed
points. As the natural next step we employ (near) scale invariance to organise
the time evolution of the competing variants within the effective description
of the _epidemic Renormalization Group_ framework. We test the resulting
theory against the time evolution of COVID-19 virus variants that validate the
theory empirically.

###### Contents

  1. 1 Introduction
  2. 2 Compartmental Models as Mutation Templates
    1. 2.1 Simple SIIR Model
    2. 2.2 Impact on the Cumulative Number of Infected
    3. 2.3 Towards RG Flows
      1. 2.3.1 Flow for $\sigma_{1,2}<1$
      2. 2.3.2 Flow for $\sigma_{1,2}>1$
      3. 2.3.3 Flow for $\sigma_{1}<1$ and $\sigma_{2}>1$
  3. 3 Mutation eRG (MeRG) Approach
    1. 3.1 The Model
    2. 3.2 Structure of the $\beta$-Functions in the $\mathbb{P}$-plane
    3. 3.3 Critical Surfaces and Mutations
    4. 3.4 Multi Wave Dynamics and CeRG
  4. 4 Explicit Examples from the COVID-19 pandemic
    1. 4.1 Second wave in California
    2. 4.2 Second and third waves in the UK
    3. 4.3 First and second waves in South Africa
  5. 5 Conclusions

##  1 Introduction

In the wake of the ongoing pandemic of SARS-CoV-2, epidemiologists are
currently witnessing a surge of data showing not only the spread and time
evolution of the virus but also its genetic evolution, notably the emergence
of mutations. While understanding their biological properties and assessing
the danger they pose for humans is of tantamount importance, the very
occurrence of new variants is also a crucial component in the time evolution
of the pandemic itself. Indeed, the companion paper [1] analyses the emergence
of new variants as one of the driving forces behind the dynamics of the SARS-
CoV-2 pandemic (notably the multi-wave structure) through the analysis of the
above mentioned data by means of machine learning and numerical techniques.
The current paper is aimed at providing a theoretical framework to model the
interplay and competition of different virus variants within a given
population in an efficient manner.

The application of mathematical modelling to describe and predict the spread
of infectious diseases has a history going back more than a century [2, 3, 4,
5, 6, 7, 8, 9, 10, 11, 12]. The pioneering SIR model [13] is an example of
compartmental models, which are based on dividing the population into
different classes and model the spread of the disease via a set of first order
differential equations in time that describe the flow of individuals between
these different compartments. Such models are deterministic in the sense that
the solutions are uniquely determined through the initial conditions
supplemented to the differential equations (apart from certain parameters
related to the infectivity of the disease and the rate at which individuals
may recover from an infection). These models can be refined by further
subdividing the population into more compartments depending on biological,
geographical and social particularities of the situation. We refer the reader
to some of the excellent reviews that are available in the literature for more
details, _e.g._ [14, 15, 16, 17, 18].

Complementary to this, there exist other models of a stochastic nature: in
these models, the microscopic processes leading to the spread of the disease
are understood in a probabilistic sense and time is typically a discretised
variable. Models of this type include lattice and percolation models, many of
which are inspired by chemical or diffusion processes. We refer the reader to
the excellent reviews [19, 20] for more details. These models are related to
the compartmental models mentioned above through processes that effectively
reduce the number of degrees of freedom, such as mean field approximations and
averaging procedures.

While very different in their original modelling of the problem, the
approaches outlined above, exhibit very interesting properties, such as
criticality and symmetries related to a re-scaling of time. The former is the
observation that in many solutions the difference between solutions where only
a very small number of individuals is concerned compared to those where a
significant fraction of the population becomes infected depends on certain
threshold parameters and the transition is rather sharp. This similarity to
phase transitions in physics has lead to further approaches using universality
classes of field theories, such as in the pioneering works [21, 22, 23, 24,
25, 26].

Furthermore, a more effective incorporation of large and short time-scale
invariance has been proposed in [27]. Here, the organisation of epidemic
curves around temporal symmetry principles was termed the epidemic
Renormalisation Group (eRG) framework [27]. It stems from the formal
identification of the running of the coupling strength in theories of
fundamental interactions [28, 29] with the time dependence of epidemiological
quantities, such as the cumulative number of infected [27]. After
demonstrating the power of reorganising the epidemic diffusion process around
time-scale invariances, the approach has been extended to account for human
interactions and mobility across different regions of the world in [30]. When
combined with mobility data provided by Google and Apple [31] as well as US
flight data [32], the framework was used to deduce the impact of lockdowns on
the global spread of the virus. This lead to the prediction, with few months
of advance, of the advent of a second pandemic wave that started in the fall
2020 in Europe [33]. The framework has been extended to contain quasi-fixed
points [34] to provide a first fully consistent mathematical description of
multiwave pandemics. The extended framework, dubbed Complex eRG, related the
interwave period to the timing of the insurgence of the next wave [35]. Last
but not least, a slight modification of the approach has shown effective to
incorporate the first impact of the US vaccination campaign [32].

In the comprehensive work of [36] we further showed how the eRG framework is
related to the traditional compartmental models of SIR type [37], and how they
emerge from microscopic percolation models which are stochastic in nature.

The investigations above were purely epidemiological in nature, however with
the advent of massive genome sequencing a new era commenced in which, as we
shall see, these data become integral part of refined epidemiological models
of the type discussed above.

In order to acquire intuition on the virus variant diffusion and dynamics, in
Section 2 we will use modified SIR-inspired compartmental models [13]. In
particular we will describe the competition between two virus variants in
terms of a SIIR model. Focussing on the temporal flow of the cumulative
numbers of infected, the analysis points to a more efficient description in
terms of an eRG system of equations. From a theoretical physics standpoint, we
find amusing the appearance of a degeneracy in the system, specifically in
terms of the asymptotic number of infected, which can be interpreted as the
emergence of a marginal operator regulating the end-point of the flow. These
features seem peculiar to the rather special systems and will grant a deeper
understanding at a more microscopic level. It will be interesting to study
similar features in other ’microscopic’ systems _e.g._ such as percolation
models described in [25, 26]. We therefore study the virus mutation version of
the eRG (MeRG) in Section 3. We show that the flow-equations can be
efficiently described in terms of a gradient flow diagrammatically. The flow
diagram reveals the existence of fixed points to be interpreted as initial and
final stages of the cumulative number of infected individuals by different
virus variants. We discover that virus mutations and their variant evolutions
can be efficiently represented in terms of theoretical physics concepts such
as critical surfaces and (quasi)-fixed points. We also discover that virus
mutations are at the heart of switching on relevant operators from the
perspective of fixed points controlling the final stage of a single wave
pandemic. These lend a natural and profound interpretation of the complex
fixed points eRG (CeRG) [34, 35], which mathematically models multi wave
pandemics.

To further substantiate our findings, in Section 4 we use the developed theory
to describe epidemiological data of COVID-19 for several regions of the world
where variants of concern have first emerged. The data on the genome
sequencing was extracted from the online repository GISAID (gisaid.org/), the
epidemiological data for the US states was taken from the NY times github and
the epidemiological data for the other countries on Ourworldindata. We first
used the GISAID data to extract the percentage of representation of the
variants among the sequences collected at each specific date. Then, by
multiplying this percentage to the epidemiological data on the number of new
cases within the country, we estimated the number of new cases related to each
variant. This method allowed us to generate the cumulative number of cases in
order to fit it with a logistic function and extract the associated infection
rate.

##  2 Compartmental Models as Mutation Templates

The eRG framework and its potential generalisations to several pathogens
(_e.g._ virus variants or different bacteria types) are effective descriptions
of the system, capturing its dynamics via certain key fields (typically the
cumulative number of infected individuals within a population). As such, other
(‘microscopic’) degrees of freedom of the system have been ‘integrated out’
(see [36] for more details on this point of view). While this approach is very
efficient in simulating infection numbers, it is non-trivial to model new
effects of the dynamics (such as the emergence of variants via mutations). To
this end, we find it convenient to first obtain intuition by studying more
basic models, of a compartmental nature, before using the obtained knowledge
to propose eRG-like models.

###  2.1 Simple SIIR Model

As a very simple ‘first-principle’ model, we consider a compartmental model to
describe the temporal evolution of two variants of a disease within a
population of size $N$. The latter is split into four different classes
(compartments):

  * •

_Susceptible individuals:_ we denote by $N\,S(t)$ the number of individuals
who can become infected with either variant of the disease.

  * •

_Infectious individuals:_ we denote by $N\,I_{1}(t)$ and $N\,I_{2}(t)$ the
number of individuals that are currently infected with the two variants of the
disease, respectively. Individuals in these two compartments can infect
susceptible individuals, with the same variant of the disease, if getting in
contact with them, with a well-defined constant rate. We assume that it is not
possible to be infected with both variants simultaneously.

  * •

_Removed individuals:_ we denote by $N\,R(t)$ the number of individuals that
can neither be infected themselves with any of the two variants, nor infect
susceptible individuals. While this removal may not only be due to recovery
from a previous infection with either of the two invariants, we assume that
the latter grants permanent immunity with respect to both variants.

We assume that $N$ is sufficiently large such that the _relative_ numbers of
susceptible $S(t)$, infectious $I_{1,2}(t)$ and removed individuals $R(t)$ are
continuously differentiable functions of time

$\displaystyle S\,,I_{1,2}\,,R:\,\mathbb{R}_{+}\longrightarrow[0,1]\,.$ (2.1)

Without loss of generality, we place the outbreak of the epidemic at time
$t=0$. Furthermore, we consider the population to be stable in time, _i.e._ we
impose the constraint

$\displaystyle S(t)+I_{1}(t)+I_{2}(t)+R(t)=1\,,$ $\displaystyle\forall
t\in\mathbb{R}_{+}\,.$ (2.2)

The dynamics of the epidemic is described by individuals being moved between
the four compartments introduced above with certain fixed777As we shall
explain below, this assumption is one of the reasons why this model is not
very suited to describe realistic epidemics and confront real world data.
However, it simplifies the analysis while still providing some intuition on
the dynamics of two competing variants of a disease. (_i.e._ time-independent)
rates:

  * •

$\gamma_{1,2}$ are the rates at which susceptible individuals turn to infected
once in contact with an individual infected by the two variants of the disease
respectively.

  * •

$\epsilon_{1,2}$ are the removal rates at which infectious individuals
(carrying variants 1 or 2 respectively) become non-infectious. This includes
recovery from the disease as well as other removal mechanisms (_e.g._ death of
the individual).

$N\,S$$N\,I_{1}$$N\,I_{2}$$N\,R$$N\,\gamma_{1}\,S\,I_{1}$$N\,\epsilon_{1}\,I_{1}$$N\,\gamma_{2}\,S\,I_{2}$$N\,\epsilon_{2}\,I_{2}$

Figure 1: Schematic flow between the four compartments of the SIIR model.

Schematically, the flow between the compartments $(S,I_{1},I_{2},R)$ is shown
in Figure 1, and it can be described by the following system of coupled first
order differential equations

$\displaystyle\frac{dS}{dt}(t)$
$\displaystyle=-(\gamma_{1}\,I_{1}(t)+\gamma_{2}\,I_{2}(t))\,S(t)\,,$
$\displaystyle\frac{dI_{1}}{dt}(t)$
$\displaystyle=\gamma_{1}\,I_{1}(t)\,S(t)-\epsilon_{1}\,I_{1}(t)\,,$
$\displaystyle\frac{dI_{2}}{dt}(t)$
$\displaystyle=\gamma_{2}\,I_{2}(t)\,S(t)-\epsilon_{2}\,I_{2}(t)\,,$
$\displaystyle\frac{dR}{dt}(t)$
$\displaystyle=\epsilon_{1}\,I_{1}(t)+\epsilon_{2}\,I_{2}(t)\,,$ (2.3)

together with the following initial conditions

$\displaystyle S(0)=S_{0}\,,$ $\displaystyle I_{1}(0)=I_{1,0}\,,$
$\displaystyle I_{2}(0)=I_{2,0}\,,$ $\displaystyle R(0)=0\,,$ with
$\displaystyle\begin{array}[]{l}0\leq S_{0},I_{1,0},I_{2,0}\leq
1\,,\\\\[4.0pt] S_{0}+I_{1,0}+I_{2,0}=1\,.\end{array}$ (2.6)

Concretely we shall consider variant 2 to have been created through a mutation
of variant 1. To quantify the effect of the second variant on the evolution of
the pandemic, we will study the evolution of the number of infected under the
variant 1, $I_{1}(t)$ for $I_{2,0}\neq 0$, and compare the results with the
control case, defined in absence of the second variant, _i.e._ $I_{2,0}=0$. To
better connect the results in this section with the eRG approach, we focus on
the cumulative number of infected $I_{\text{c},1}(t)$, defined as

$\displaystyle
I_{\text{c},i}(t)=N\,I_{0,i}+N\,\gamma_{i}\,\int_{0}^{t}I_{i}(t^{\prime})\,S(t^{\prime})\,dt^{\prime}\,,$
$\displaystyle\forall i=1,2\,.$ (2.7)

We remark that the model introduced above does not trace the origin of the
variant 2 of the virus as a mutation of the first one. The latter could be
encoded in a more evolutionary version of the 4-compartment model, where a
mutated variant could appear at some time $t_{0}>0$. To be concrete, the
evolutionary model can be described by the following set of differential
equations:

$\displaystyle\frac{dS}{dt}=-(\gamma_{1}\,I_{1}+\gamma_{2}\,I_{2})\,S\,,$ with
$\displaystyle S(0)=S_{0}\,,I_{1}(0)=I_{1,0}\,,I_{2}(0)=R(0)=0$
$\displaystyle\frac{dI_{1}}{dt}=\gamma_{1}\,I_{1}\,S-\epsilon_{1}\,I_{1}-\beta(t)\,I_{1}\,,$
$\displaystyle 0\leq S_{0},I_{1,0},\leq 1$
$\displaystyle\frac{dI_{2}}{dt}=\gamma_{2}\,I_{2}\,S-\epsilon_{2}\,I_{2}+\beta(t)\,I_{1}\,,$
$\displaystyle S_{0}+I_{1,0}=1$
$\displaystyle\frac{dR}{dt}=\epsilon_{1}\,I_{1}+\epsilon_{2}\,I_{2}\,,$ (2.8)

where $\beta(t)$ is a time-dependent rate with point-like support that
converts at time $t_{0}>0$ (the instant in which the variant appears) a
fraction of the infectious individuals of variant 1 into those of variant 2.
This function is defined as:

$\displaystyle\beta(t)=\left\\{\begin{array}[]{lcl}\beta_{0}&\text{if}&t=t_{0}\,,\\\
0&\text{if}&t\neq t_{0}\,.\end{array}\right.$ (2.11)

However, since we are excluding the possibility of a re-infection of removed
individuals (with either of the two variants) and we are only interested in
the impact of the appearance of the mutation on the dynamics of variant 1
(_i.e._ the evolution for $t>t_{0}$), the above SIIR model (2.1) is equivalent
to the simpler system (2.3) for a suitable choice of the boundary conditions
(2.6).

###  2.2 Impact on the Cumulative Number of Infected

While an analytic solution of the system (2.3) and (2.6) appears very
difficult, we can study numerical solutions using a simple Euler integration
method. Since to this end we have to consider fixed values for $\gamma_{1,2}$
and $\epsilon_{1,2}$, we first study the dependence of $I_{\text{c},1}(t)$ on
$\gamma_{2}$ and $\epsilon_{2}$. In analogy to the SIR model, it is convenient
to define the reproduction numbers of the two variants as follows:

$R_{0,1}\equiv\sigma_{1}=\frac{\gamma_{1}}{\epsilon_{1}}\,,\qquad
R_{0,2}\equiv\sigma_{2}=\frac{\gamma_{2}}{\epsilon_{2}}\,.$ (2.12)

As we will see, the time dependence of the cumulative number of infected
mainly depends on the two reproduction numbers instead of the detailed values
of infection and recovery rates, under certain conditions.

This point is illustrated in Figure 2, where we show solutions
$I_{\text{c},1}(t)$ for fixed $\sigma_{1}=1.2$ and for four values of
$\sigma_{2}$. The various curves in each panel correspond to different values
of $\gamma_{2}$ (where $\epsilon_{2}$ also varies accordingly to keep
$\sigma_{2}$ fixed). These results, along with further checks we have
performed, suggest that the cumulative number of infected by the original
variant only depends on $\sigma_{2}$ as long as its value is similar to that
of $\sigma_{1}$, i.e. if the two variants have similar reproduction numbers
(quantitatively, $\sigma_{2}$ can be at most $\sim 50\%$ larger than
$\sigma_{1}$). Conversely, for $\sigma_{2}\gg\sigma_{1}$, the solutions of the
SIIR model are sensitive to the specific values of $\gamma_{2}$ and
$\epsilon_{2}$. The results in Fig. 2 also highlight that the shape of
$I_{\text{c},1}(t)$, or equivalently its asymptotic value at large times, is
substantially affected by the presence of the second variant only is
$\sigma_{2}$ is significantly larger than $\sigma_{1}$.

Figure 2: Time evolution of $I_{\text{c},1}$ as a function of $\gamma_{2}$
with $\sigma_{2}=0.75$ (top left panel), $\sigma_{2}=1.25$ (top right panel),
$\sigma_{2}=1.5$ (bottom left panel) and $\sigma=2$ (bottom right panel). In
all cases we have used $\gamma_{1}=0.03$ and $\epsilon_{1}=0.025$ (such that
$\sigma_{1}=1.2$ ) as well as $S_{0}=0.9$, $I_{1,0}=0.099$ and
$I_{2,0}=0.001$.

To quantify the observations above, we can use the fact that the numerical
solutions for $I_{\text{c},1}(t)$ are fairly well approximated by a logistic
function of the form

$\displaystyle I_{\text{c},1}(t)\sim\frac{N\,A}{1+e^{-\lambda(t-\kappa)}}\,,$
with $\displaystyle\begin{array}[]{lcl}A\in[0,1]\,,\\\\[4.0pt]
\lambda,\kappa\in\mathbb{R}\,,\end{array}$ (2.15)

for different values of the parameters $(A,\lambda,\kappa)$. Here, $N\,A$ is
the asymptotic value of the cumulative number of infected

$\displaystyle N\,A=\lim_{t\to\infty}\,I_{\text{c},1}(t)\,,$ (2.16)

while $\lambda$ is an effective infection rate, _i.e._ a measure of how
infectious the variant is (and is generally related to $\sigma_{1}$), while
$\kappa$ is a parameter that allows to shift the outbreak of the spread of
variant 1. Exemplary fits of numerical solutions of $I_{\text{c},1}(t)$ for
different values of $(\sigma_{1},\sigma_{2})$ along with the numerical fitting
parameters are shown in Figure 3. We remark, however, that fitting the
cumulative number of infected individuals with logistic functions is not the
only viable approximation and we shall encounter other possibilities
(depending on specific cases) in Section 2.3.

Figure 3: Numerical solution (blue) and fitted logistic function (red) of
$I_{\text{c},1}$ as a function of time for two different combinations of
$\sigma_{1,2}$: top row:
$(\sigma_{1},\sigma_{2})=(0.75,0.5)\,,(0.75,1.25)\,,(0.75,2)$, middle row:
$(\sigma_{1},\sigma_{2})=(1.5,0.5)\,,(1.5,1.25)\,,(1.5,2)$, bottom row:
$(\sigma_{1},\sigma_{2})=(2,0.5)\,,(2,1.25)\,,(2,2)$. All plots use
$S_{0}=0.9$, $I_{1,0}=0.099$, $I_{2,0}=0.001$.

Finally, using the result that $I_{\text{c},1}$ can be approximated using a
logistic function, we can quantify more concretely the impact of the
appearance of variant 2. To this end, we plot in Figure 4 the values
$(A,\lambda,\kappa)$ of the approximations of $I_{\text{c},t}$ as functions of
$\sigma_{2}$. These examples suggest that an existing variant of a disease is
strongly impacted by the appearance of a new variant only if the latter is
significantly more infectious, _i.e._ if $\sigma_{2}$ is much larger than
$\sigma_{1}$ (in many cases approximately by a factor of 2).

Figure 4: Fit parameters $(A\,,\lambda\,,\kappa)$ for different values of
$\sigma_{1}$ as functions of $\sigma_{2}$. All plots are shown for
$S_{0}=0.9$, $I_{1,0}=0.099$, $I_{2,0}=0.001$.

We can summarise the basic findings suggested by the simple model (2.3) and
(2.6) through the following three points:

  * •

Despite not being an exact solution, logistic functions (sigmoids) are good
approximations to describe competing variants of a disease.

  * •

Key parameters of a variant in competition with a second variant that is not
significantly more infectious can be described, to first approximation, by the
reproduction numbers $\sigma_{1,2}$ alone (rather than individually
$\gamma_{1,2}$ and $\epsilon_{1,2}$).

  * •

In order to have a significant impact on an existing variant (_i.e._ to change
key parameters such as the asymptotic numbers of infected individuals, the
effective infection rates, _etc_.), a new variant needs to have a much higher
reproduction number, $\sigma_{2}\gg\sigma_{1}$.

###  2.3 Towards RG Flows

Having understood the impact of variants on the evolution of the infections in
a simple SIIR model, we now want to understand the dynamics of the system in
analogy to the renormalisation group in particle physics [28, 29]. This is a
preliminary step that will help us gear up towards the inclusion of variants
in the eRG framework. In practice, the solutions of the system of equations
(2.3) and (2.6) can be thought of as describing the flow of the physical
system (here represented by the cumulative number of infected by the two
variants) in the plane

$\displaystyle\mathbb{P}=\\{(I_{\text{c},1},I_{\text{c},2})\in[0,N]\times[0,N]|I_{\text{c},1}+I_{\text{c},2}\leq
N\\}\,.$ (2.17)

The trajectories are parameterised by time $t$. Following the eRG approach
[27], we shall try to interpret these trajectories as renormalisation group
flows from a repulsive fixed point at $(I_{\text{c},1},I_{\text{c},2})=(0,0)$
to a new fixed point $(I_{\text{c},1},I_{\text{c},2})\neq(0,0)$ in the plane
$\mathbb{P}$. In this regard, we shall consider the initial conditions (2.6)
as small deformations away from the initial fixed point. Keeping $S_{0}$
fixed, this gives a line (a co-dimension 1 surface) of starting points of
trajectories in the plane $\mathbb{P}$, which can be found by solving (2.3)
with the initial conditions satisfying

$\displaystyle S_{0}=1-I_{1,0}-I_{2,0}\,.$ (2.18)

We then plot $(I_{\text{c},1}(t),I_{\text{c},2}(t))$ for successive times
$t\in[0,\infty)$. The result is shown schematically, for different values of
$\sigma_{1,2}$, in Figure 5. We remark that here (and in the following) we
have assumed $\epsilon_{1}=\epsilon_{2}$, _i.e._ that the removal rates for
both variants are the same.888As we have seen in Section 2.2, choosing
different $\epsilon_{1,2}$ leads to only small deviations, as long as
$\sigma_{1}$ and $\sigma_{2}$ are not too different from one another.

Figure 5: Trajectories for different initial conditions in the $\mathbb{P}$
plane for $S_{0}=0.999$ and $(\sigma_{1},\sigma_{2})=(0.85,0.95)$ (left panel)
and $(\sigma_{1},\sigma_{2})=(1.075,1.1)$ (right panel).

Before analysing in more detail the cases of $\sigma_{1,2}<1$ and
$\sigma_{1,2}>1$ (which show certain qualitative differences), we first remark
that the flows, _i.e._ the lines in Figure 5, do not continue indefinitely but
end at determined points in the $\mathbb{P}$-plane. These can be considered
fixed points of the flow, as the evolution of the epidemic stops once these
points are reached. More concretely, we observe that the line of initial
conditions flows to another line of end-points. In fact, from the point of
view of the SIIR model, the endpoints are equivalent, and correspond to the
same final state of the system. To understand this, we recall that the
asymptotic solutions of the model (2.3) and (2.6) corresponds to the
eradication of both variants

$\displaystyle\lim_{t\to\infty}I_{1}(t)=0=\lim_{t\to\infty}I_{2}(t)\,.$ (2.19)

Therefore at the endpoints in Figure 5 (which correspond to $t\to\infty$),
only susceptible and removed individuals remain in the system and

$\displaystyle\lim_{t\to\infty}(I_{\text{c},1}+I_{\text{c},2})(t)=\lim_{t\to\infty}N\,R(t)\,.$
(2.20)

Hence, the co-dimension 1 line of end-points in Figure 5 reflects different
weighted re-distributions of the removed individuals into $I_{\text{c},1}$ and
$I_{\text{c},2}$. In the SIIR model, however, the removed individuals are
indistinguishable as they are collected in a single compartment $R$, and all
the final configurations in the end-point line correspond to the same type of
final state of the SIIR compartments (_i.e._ with $I_{1}=I_{2}=0$). In other
words, while the flow lines keep track of which individuals have been infected
with which variant, this distinction becomes irrelevant at the level of the
SIIR system when all previously infected individuals have become
(indistinguishable) removed individuals. We shall see that this redundancy of
the endpoints of the ’flow’ of the system has a very natural interpretation
from the perspective of the epidemic renormalisation group.

####  2.3.1 Flow for $\sigma_{1,2}<1$

We shall now analyse in more detail the trajectories of the system in the
$\mathbb{P}$-plane and we first focus on the case

$\displaystyle\sigma_{1}:=\frac{\gamma_{1}}{\epsilon_{1}}<1\,,$ and
$\displaystyle\sigma_{2}:=\frac{\gamma_{2}}{\epsilon_{2}}<1\,.$ (2.21)

In the simple SIR model [13], the values of $\sigma_{1,2}$ correspond to the
basic reproduction numbers of each variant. Physically, they tell the average
number of (susceptible) individuals who are infected by a single infectious
one during the period in which the latter is infectious. The condition
$\sigma_{i}<1$, therefore, implies that the number of infectious individuals
$I_{i}(t)$ is a monotonically decreasing function of time, since
$\frac{dI_{i}}{dt}(t)<0$ $\forall t>0$, as can be seen from eq. (2.3).

As is visible in the left panel of Figure 5, under the condition (2.21) the
end points of the trajectories of the system follow an approximate linear
relation in the $\mathbb{P}$-plane. Concretely, upon defining the positions of
the end-points as

$\displaystyle
I_{\text{c},i}^{\text{asym}}(I_{1,0},I_{2,0})=\frac{1}{N}\,\lim_{t\to\infty}I_{\text{c},i}(t)\,,$
$\displaystyle\forall i=1,2\,,$ (2.22)

which depend on the initial conditions $(S_{0},I_{1,0},I_{2,0})$ subject to
(2.18), we find to good approximation (see Figure 6 for an example) the
relation

$\displaystyle I_{\text{c},2}^{\text{asym}}(I_{1,0},I_{2,0})\sim
a\,I_{\text{c},1}^{\text{asym}}(I_{1,0},I_{2,0})+b\,,$ with
$\displaystyle\begin{array}[]{l}a=-\frac{1+\frac{1}{\sigma_{2}}\,W(-S_{0}\sigma_{2}\,e^{\sigma_{2}})}{1+\frac{1}{\sigma_{1}}\,W(-S_{0}\sigma_{1}\,e^{\sigma_{1}})}\,,\\\\[10.0pt]
b=1+\frac{1}{\sigma_{2}}\,W(-S_{0}\sigma_{2}\,e^{\sigma_{2}})\,,\end{array}$
(2.25)

where $W$ is the Lambert function. The coefficients $(a,b)$ can be inferred
from the points $(I_{\text{c},1}^{\text{asym}}(1-S_{0},0),0)$ and
$(0,I_{\text{c},2}^{\text{asym}}(0,1-S_{0}))$, which can be determined
analytically from the SIR model [13] with 3 compartments (_i.e._ with only one
variant of the disease).

Figure 6: Left panel: Approximation (2.25) of the codimension 1 surface
(dashed black line) that represents the end-points of the trajectories in the
$\mathbb{P}$-plane. Each red dot represents
$(I_{\text{c},1}^{\text{asym}},I_{\text{c},2}^{\text{asym}})$ at asymptotic
time $t\to\infty$ for different choices of $(I_{1,0},I_{2,0})$ satisfying
(2.18). Right panel: Numerical approximations (dashed black lines) for equal
time-slices $(I_{\text{c},1}(t),I_{\text{c},2}(t))$ for finite $t$ (red dots).
Both plots use $S_{0}=0.999$ and $(\sigma_{1},\sigma_{2})=(0.85,0.95)$.

In fact, even at finite time $t<\infty$, we have found numerical evidence for
a linear relation between $I_{\text{c},1}(t)$ and $I_{\text{c},2}(t)$ with
different initial conditions $(I_{1,0},I_{2,0})$: the right panel of Figure 6
shows equal time slices of the flows in the $\mathbb{P}$-plane, _i.e._
$(I_{\text{c},1}(t),I_{\text{c},2}(t))$ with different $(I_{1,0},I_{2,0})$
evaluated at the same time $t<\infty$ (red points), which are approximated by
linear curves (dashed black lines).

Next, in order to make contact with a possible eRG approach to describe the
flow of the system in the $\mathbb{P}$-plane, we consider the time derivatives
of $(I_{c,1}(t),I_{c,2}(t))$. According to the definition (2.7) they are given
by

$\displaystyle\frac{dI_{\text{c},i}}{dt}(t)=N\,\gamma_{i}\,I_{i}(t)\,S_{i}(t)\,,$
$\displaystyle\forall\,i=1,2\,.$ (2.26)

After dividing by the total size of the population for simplicity, from
studying the numerical solutions of (2.3) and (2.6) we observe that these time
derivatives can be well approximated by linear functions in $I_{\text{c},i}$

$\displaystyle-\beta_{i}(I_{\text{c},j}(t)):=\frac{1}{N}\,\frac{dI_{\text{c},i}}{dt}(t)\sim\lambda_{i}\left(1-\frac{I_{\text{c},i}(t)}{N\,A_{i}}\right)\,,$
$\displaystyle\forall\,i=1,2\,,$ (2.27)

where $(\lambda_{i},A_{i})$ are constants that are implicitly functions of the
initial conditions $(I_{1,0},I_{2,0})$. Indeed, examples of the approximation
are shown in Figure 7. Here, in order to better gauge the impact of the
initial conditions, we have parametrised them as

$\displaystyle I_{1,0}=(1-S_{0})\,r\,,$ $\displaystyle
I_{2,0}=(1-S_{0})(1-r)\,,$ with $\displaystyle r\in[0,1]\,.$ (2.28)

In fact, one can also represent the trajectory of the system (for fixed
initial conditions $(I_{1,0},I_{2,0})$) in the $\mathbb{P}$-plane as following
the vector field
$\frac{1}{N}(\tfrac{dI_{\text{c},1}}{dt},\tfrac{dI_{\text{c},2}}{dt})$, as is
schematically shown in Figure 8. This plot better illustrates the concept of
flow of the system in the two dimensional plane: the system slows down as it
approaches the fixed point. This is represented by the colour code of the flow
lines (arrows in Figure 8), matching the length of the derivative vector and
is also visible from the fact that the black dots of the trajectory (which
represent the system after equal time intervals) become denser as they
approach the fixed point.

Figure 7: Linear approximation of
$(I_{\text{c},1},\tfrac{dI_{\text{c},1}}{dt})$ (left panel) and
$(I_{\text{c},2},\tfrac{dI_{\text{c},2}}{dt})$ (right panel) as a function of
the parameter $r$ (defined in eq. (2.28)). The red dots represent numerical
solutions of $(I_{\text{c},i},\tfrac{dI_{\text{c},i}}{dt})$ for different
values of time $t$, grouped according to the initial conditions. All plots use
$S_{0}=0.999$ and $(\sigma_{1},\sigma_{2})=(0.85,0.95)$. Figure 8: Trajectory
of the system in the $\mathbb{P}$-plane for a fixed initial condition (black
line) following the vector field
$\frac{1}{N}(\tfrac{dI_{\text{c},1}}{dt},\tfrac{dI_{\text{c},2}}{dt})$ (the
colour of the vectors represents the length). We have used the numerical
parameters $r=\tfrac{1}{2}$, $S_{0}=0.999$ and
$(\sigma_{1},\sigma_{2})=(0.85,0.95)$.

We can express the approximated derivatives (2.27) in the following form:

$\displaystyle-\beta_{i}(I_{\text{c},j})=\frac{1}{N}\,\frac{dI_{\text{c},i}}{dt}(t)\sim\nabla_{i}\,\Phi(I_{\text{c},j}(t))\,,$
with
$\displaystyle\Phi(I_{\text{c},j})=\sum_{j=1}^{2}I_{\text{c},j}\,\lambda_{j}\left(1-\frac{I_{\text{c},j}}{2NA_{j}}\right)$
(2.29)

where we have used that $\beta_{i}$ in (2.27) only depends on $I_{\text{c},i}$
and we have defined the gradient operator in the $\mathbb{P}$-plane

$\displaystyle\nabla_{i}:=\frac{\partial}{\partial I_{\text{c},i}}\,,$
$\displaystyle\forall i=1,2\,.$ (2.30)

Hence, for fixed initial conditions, the trajectory of the system in the
$\mathbb{P}$-plane can (approximately) be described as a gradient flow. Notice
in this regard that $\Phi$ (through $(\lambda_{i},A_{i})$), also implicitly
depends on $(I_{1,0},I_{2,0})$.

Before moving on to the cases $\sigma_{1,2}>1$, in passing we make three more
remarks:

  * •

We find evidence for the fact that $(\lambda_{i},A_{i})$ depend also on
$\sigma_{i}$, but not on $\sigma_{j\neq i}$. Concretely, the numerical
solutions can be approximated by

$\displaystyle\lambda_{i}\sim\tilde{\lambda}_{i}\,r\,,$ $\displaystyle
A_{i}\sim\tilde{A}_{i}\,r\,\left(1+\frac{1}{\sigma_{i}}W(-S_{0}\,\sigma_{i}\,e^{-\sigma_{i}})\right)\,,$
(2.31)

where $(\tilde{\lambda}_{i},\tilde{A}_{i})$ are independent of $\sigma_{j\neq
i}$.

  * •

The differential equation (2.27)

$\displaystyle\frac{dI_{\text{c},i}}{dt}(t)=\lambda_{i}\,N\left(1-\frac{I_{\text{c},i}}{NA_{i}}\right)\,,$
with $\displaystyle I_{\text{c},i}(t=0)=I_{1,0}\,,$ (2.32)

can be solved analytically

$\displaystyle
I_{\text{c},i}(t)=A_{i}\,N+(I_{1,0}-N\,A_{i})\,e^{-\frac{\lambda_{i}}{A_{i}}\,t}\,.$
(2.33)

  * •

As we have explained, varying the initial conditions $(I_{1,0},I_{2,0})$ while
keeping $S_{0}$ fixed as in (2.18) leads to a codimension 1 surface of fixed
points. However, by changing at the same time $(\sigma_{1},\sigma_{2})$, it is
possible to find a family of trajectories that flow to a single point in the
$\mathbb{P}$-plane (an example is schematically shown in Figure 9). This can
be achieved for any choice of $\sigma_{1,2}$ and is not limited to
$\sigma_{1,2}<1$.

Figure 9: Modified trajectories where in addition to the initial condition
also $(\sigma_{1},\sigma_{2})$ have been modified in such a way to lead to a
single point in the $\mathbb{P}$-plane.

####  2.3.2 Flow for $\sigma_{1,2}>1$

We now move to analysing the case where both $\sigma_{1,2}>1$. The first
difference from the previous case is that neither the asymptotic end points
$(I_{\text{c},i}^{\text{asym}},I_{\text{c},2}^{\text{asym}})$ (see (2.22) for
the definition) of the trajectories in the $\mathbb{P}$-plane nor equal-time
slices for finite $t$ can be approximated by linear relations. As an
illustration, Figure 10 shows the corresponding plots for
$(\sigma_{1},\sigma_{2})=(1.2,1.3)$, highlighting visible deviations from a
linear regime. We notice that these deviations become more significant the
larger is the difference between $\sigma_{1}$ and $\sigma_{2}$.

Figure 10: Left panel: endpoints of the trajectories of the system in the
$\mathbb{P}$-plane for different initial conditions. Each red dot represents
$(I_{\text{c},1}^{\text{asym}},I_{\text{c},2}^{\text{asym}})$ at asymptotic
time $t\to\infty$ for different choices of $(I_{1,0},I_{2,0})$ satisfying
(2.18). The black lines correspond to (linear) approximations for initial
conditions leading to $I_{\text{c},1,2}^{\text{asym}}\ll 1$ respectively.
Right panel: Equal time-slices $(I_{\text{c},1}(t),I_{\text{c},2}(t))$ for
finite $t$ (red dots) show deviations from approximate linear relations
between $(I_{\text{c},1},I_{\text{c},2})$. Both plots use $S_{0}=0.999$ and
$(\sigma_{1},\sigma_{2})=(1.2,1.3)$.

We now turn to the differential equations (2.26) for $I_{\text{c},i}(t)$. The
main difference compared to the case $\sigma_{1,2}<1$ is that the factors
$N\,\gamma_{i}\,I_{i}(t)\,S_{i}(t)$ can no longer be approximated by linear
functions in $I_{\text{c},i}$. Instead, we find evidence that the time-
derivatives can be approximated by a quadratic function of the form

$\displaystyle-\beta_{i}(I_{\text{c},j}(t)):=\frac{1}{N}\,\frac{dI_{\text{c},i}}{dt}(t)\sim\lambda_{i}\,I_{\text{c},i}(t)\,\left(1-\frac{I_{\text{c},i}(t)}{N\,A_{i}}\right)+\delta_{i}\,,$
$\displaystyle\forall i=1,2\,,$ (2.34)

where $(\lambda_{i},A_{i},\delta_{i})$ are constants that implicitly depend on
the initial conditions $(I_{1,0},I_{2,0})$ (we remark that $\delta_{i}\ll 1$
is generically a rather small parameter, since it is related to the initial
conditions $I_{i,0}$). Indeed, examples of the approximations as functions of
the initial conditions parametrised by $r$, as defined in Eq. (2.28), are
shown in Figure 11. By varying $\sigma_{1,2}$, we find that the quadratic
approximation (2.34) becomes less satisfactory the larger is the difference
between $\sigma_{1,2}$.

Figure 11: Quadratic approximation of
$(I_{\text{c},1},\tfrac{dI_{\text{c},1}}{dt})$ (left panel) and
$(I_{\text{c},2},\tfrac{dI_{\text{c},2}}{dt})$ (right panel) as a function of
the parameter $r$ (defined in eq. (2.28)). The red dots represent numerical
solutions of $(I_{\text{c},i},\tfrac{dI_{\text{c},i}}{dt})$ for different
values of time $t$, grouped according to the initial conditions. All plots use
$S_{0}=0.999$ and $(\sigma_{1},\sigma_{2})=(1.2,1.3)$.

As before, we can also represent the trajectory of the system (for fixed
initial conditions $(I_{1,0},I_{2,0})$) in the $\mathbb{P}$-plane as following
the vector field
$\frac{1}{N}(\tfrac{dI_{\text{c},1}}{dt},\tfrac{dI_{\text{c},2}}{dt})$, as
schematically shown in Figure 12. Using (2.34), we can therefore write

$\displaystyle-\beta_{i}(I_{\text{c},j})=\frac{1}{N}\,\frac{dI_{\text{c},i}}{dt}(t)\sim\nabla_{i}\,\Phi(I_{\text{c},j}(t))\,,$
with
$\displaystyle\Phi(I_{\text{c},j})=\sum_{j=1}^{2}\left[I_{\text{c},j}^{2}\,\frac{\lambda_{j}}{2}\left(1-\frac{2I_{\text{c},j}}{3NA_{j}}\right)+\delta_{j}\,I_{\text{c},j}\right]\,,$
(2.35)

where we have used the fact that $\beta_{i}$ in (2.34) only depends on
$I_{\text{c},i}$. The gradient operator in the $\mathbb{P}$-plane $\nabla_{i}$
is defined in (2.30). Hence, for fixed initial conditions, the trajectory of
the system in the $\mathbb{P}$-plane can be (approximately) described as a
gradient flow. Notice in this regard that $\Phi$ also implicitly depends on
$(I_{1,0},I_{2,0})$ through $(\lambda_{i},A_{i})$.

Figure 12: Trajectory of the system in the $\mathbb{P}$-plane for a fixed
initial condition (black line) following the vector field
$\frac{1}{N}(\tfrac{dI_{\text{c},1}}{dt},\tfrac{dI_{\text{c},2}}{dt})$ (the
colour of the vectors represents the length). We have used the numerical
parameters $r=\tfrac{1}{2}$, $S_{0}=0.999$ and
$(\sigma_{1},\sigma_{2})=(1.2,1.3)$.

####  2.3.3 Flow for $\sigma_{1}<1$ and $\sigma_{2}>1$

We briefly comment on the case where one of the $\sigma_{i}$ (without loss of
generality, we can choose it to be $\sigma_{1}$) is smaller than 1 and the
other one larger than 1. This corresponds to a situation, where $I_{1}(t)$ is
a monotonically decreasing function of time, while $I_{2}(t)$ has a maximum
before tending to 0 for large $t$.

The situation is a combination of the cases discussed in the previous
subsections. Notably, the time evolution of $I_{\text{c},1}$ can be
approximated by an equation of the form (2.29) with a function $\Phi$ that is
quadratic in $I_{\text{c},1}$, while the time evolution of $I_{\text{c},2}$
can be approximated by an equation of the form (2.35) with a $\Phi$ that is
cubic in $I_{\text{c},2}$. We remark that, as before, the approximations are
less satisfactory the larger is the difference between $\sigma_{1}$ and
$\sigma_{2}$.

##  3 Mutation eRG (MeRG) Approach

The SIIR model discussed in the previous Section is a simple model that allows
to gain basic intuition about the dynamics of two competing variants of a
disease, without assuming too much about basic ‘microscopic’ processes that
govern its spread. However, due to the simplicity (notably the fact that
infection and recovery rates are assumed to be constant in time), this model
is not particularly useful to confront (or even predict) real world data. In
principle, it is possible to extend the model by allowing for time-dependent
$(\gamma_{i},\epsilon_{i})$ and/or adding additional compartments, however at
the price of complicating the analysis and loosing predicting power.

Hence, in the following we will use the intuition we have gained about the
pandemic with two variants to extend the eRG approach. The latter is an
economical effective approach that entails a high degree of predictivity in
terms of the time structure of the pandemic. As discussed in [27], the eRG
framework makes use of particular symmetries in the time evolution of an
epidemic to give a simplified description of certain key quantities (namely
the cumulative number of infected individuals), while more microscopic degrees
of freedom have already been taken into account (or ‘integrated out’). See
[36] for a more detailed description of this philosophy.

###  3.1 The Model

The eRG approach consists in defining $\beta$-functions that govern the time
evolution of the system at the global level. To this end, in analogy to
particle physics, for each variant we first define an epidemic ‘coupling
strength’ as a monotonic, differentiable and bijective, function
$\alpha_{i}=f_{i}(I_{\text{c},i})$ of the cumulative number of infected
individuals. Different choices of $f_{i}$ correspond to different
renormalisation schemes and we expect physical results in general not to
depend on the choice. The $\beta$-functions are then defined as

$\displaystyle-\beta_{i}=\frac{d\alpha_{i}}{dt}=\sum_{j=1}^{2}\frac{df_{i}}{dI_{\text{c},j}}\,\frac{dI_{\text{c},j}}{dt}(t)\,,$
$\displaystyle\forall i=1,2\,.$ (3.1)

The intuition obtained from the SIIR model suggests to model these
$\beta$-functions as polynomials in $I_{\text{c},i}$, which is either linear
(for variants with a basic reproduction number $\sigma_{i}<1$) or quadratic
(for variants with a basic reproduction number $\sigma_{i}>1$). In fact,
comparing the equations that emerged in the context of the SIIR model, _i.e._
(2.27) with (2.34) (and neglecting999As we have remarked above, in the SIIR
model the (small) parameters $\delta_{i}$ are related to the (finite) initial
conditions $I_{i,0}$ of the (active) number of infectious individuals which is
necessary to start the time evolution of the system. In the context of the
eRG, as we shall discuss in more detail below, we model the point
$(I_{\text{c},1},I_{\text{c},2})=(0,0)$ (_i.e._ the absence of the disease) as
a fixed point, albeit a repulsive one, which requires $\delta_{i}=0$. the
small coefficient $\delta_{i}$), we see that both can in fact be brought into
the same framework (3.1) simply by different choices of the functions $f_{i}$,
namely

$\displaystyle
f_{i}(I_{\text{c},i})=\left\\{\begin{array}[]{lcl}c_{i}\,\ln(I_{\text{c},i})&\text{for}&\sigma_{i}<1\,,\\\
c_{i}\,I_{\text{c},i}&\text{for}&\sigma_{i}>1\,,\end{array}\right.$ (3.4)

for some constant $c_{i}$. Moreover, (2.29) and (2.35) suggest to further
model the beta function in the form of a gradient equation, _i.e._

$\displaystyle-\beta_{i}=\frac{d\alpha_{i}}{dt}=\sum_{j=1}^{2}\frac{df_{i}}{dI_{\text{c},j}}\,\frac{dI_{\text{c},j}}{dt}(t)=\nabla_{i}\,\Phi(I_{\text{c},_{j}})\,.$
(3.5)

where the gradient operator is defined in (2.30) and $\Phi$ is a quadratic or
cubic function, respectively. In fact, taking these results together with the
liberty of scheme-redefinitions and comparing with the general form of the eRG
for single variants (see [34]) we are naturally led to choose for simplicity
$f_{i}(I_{\text{c},i})=I_{\text{c},i}$ and model $\Phi$ as

$\displaystyle\Phi(I_{\text{c},j})=\sum_{j=1}^{2}I_{\text{c},j}^{2}\,\frac{\lambda_{j}}{2}\left(1-\frac{2I_{\text{c},j}}{3NA_{j}}\right)\,.$
(3.6)

In practice, this means that we are modelling the time-evolution of each
variant as independent of the other, except for a potential change in the
parameters $(\lambda_{i},A_{i})$. From the intuition obtained from the SIIR
model, this is justified as long as the two variants do not have effective
reproduction numbers that are largely different from one another (see Section
2.2). We shall see that this assumption also leads to reasonable results
compared to real world data. Indeed, the system (3.5) for $\Phi$ given in
(3.6) allows for an analytic solution of $(I_{\text{c},1},I_{\text{c},2})(t)$
which can be written in terms of logistic functions

$\displaystyle
I_{\text{c},i}(t)=\frac{NA_{i}}{1+e^{-\lambda_{i}(t-\kappa_{i})}}\,,$
$\displaystyle\forall i=1,2\,,$ (3.7)

which well reproduce the data for the first wave of the COVID-19 pandemic
[27]. Here $\kappa_{1,2}$ are integration constants that determine which
trajectory the system follows in the $\mathbb{P}$-plane. These parameters
resemble the initial conditions $(I_{1,0},I_{2,0})$ from the perspective of
the SIIR model.

###  3.2 Structure of the $\beta$-Functions in the $\mathbb{P}$-plane

The vector field $(-\beta_{1}(I_{\text{c},i}),-\beta_{2}(I_{\text{c},i}))$ in
the $\mathbb{P}$-plane is schematically plotted in Figure 13. It has four
fixed points, _i.e._ points $(I_{\text{c},1}^{(k)},I_{\text{c},2}^{(k)})$ (for
$k=0,1,2,3$) where

$\displaystyle\beta_{i}(I_{\text{c},1}^{(k)},I_{\text{c},2}^{(k)})=0\,,$
$\displaystyle\forall i=1,2\,,$ (3.8)

explicitly given by

$\displaystyle P_{0}=(I_{\text{c},1}^{(0)},I_{\text{c},2}^{(0)})=(0,0)\,,$
$\displaystyle
P_{1}=(I_{\text{c},1}^{(1)},I_{\text{c},2}^{(1)})=(NA_{1},0)\,,$
$\displaystyle
P_{2}=(I_{\text{c},1}^{(2)},I_{\text{c},2}^{(2)})=(0,NA_{1})\,,$
$\displaystyle
P_{3}=(I_{\text{c},3}^{(k)},I_{\text{c},2}^{(3)})=(NA_{1},NA_{2})\,.$ (3.9)

Figure 13: Schematic structure of the RG-equations in the $\mathbb{P}$-plane
(the values used for the plot are $N\lambda_{1}=0.2$, $N\lambda_{2}=0.25$ and
$A_{1}=0.3$, $A_{2}=0.4$). The $\beta$-functions exhibit four fixed points:
$P_{0,1,2,3}$ of which $P_{0}$ is repulsive in all directions, $P_{1,2}$ have
one attractive direction and one repulsive one, while $P_{3}$ is attractive in
all directions. Different trajectories as solutions of (3.5) connecting the
fixed points are indicated by the dotted lines. The colouring of the vectors
indicates the norm $||(\beta_{1},\beta_{2})||$, which respectively leads to
smaller or larger distances between different points of the flow lines.

Among them, $P_{0}$ is repulsive in all direction (_i.e._ in Figure 13 all
arrows point away from it) and corresponds to the case where no disease is
(and never has been) present. In fact, moving away from this fixed point by
infecting even only a small number of individuals of the population (with
either of the two variants) causes the system to flow to one of the other
three fixed points. Among them, $P_{1,2}$ are repulsive in one direction, but
attractive in the other. Since they are characterised by $I_{\text{c},2}=0$ or
$I_{\text{c},1}=0$ respectively, they correspond to the endpoints of scenarios
in which variant 2 or variant 1 is never present in the population (_i.e._ all
infected individuals are infected with only one of the two variants). These
fixed points can be reached only by the flow lines represented in black in
Figure 13, which are initiated by a deformation away from $P_{0}$ along
$I_{\text{c},1}$ or $I_{\text{c},2}$ only, respectively.101010We shall discuss
in the following subsection 3.3 in more detail the case in which along one of
the black trajectories in Figure 13 infectious individuals of the other
variant appear. This scenario models for example the appearance of a mutation
of an already present variant. In this case the system will flow to the fixed
point $P_{3}$ (rather than continue towards $P_{1}$ or $P_{2}$.) Any
deformation that switches on both $I_{\text{c},1}$ and $I_{\text{c},2}$
(_i.e._ any scenario in which infected with both variants are present in the
population) causes the system to flow to fixed point $P_{3}$, an example of
such a flow is indicated in red in Figure 13. Which trajectory is realised
depends on the initial deformation, which is represented by the parameters
$\kappa_{1,2}$ in the solution (3.7) (and which is akin to the choice of
different initial conditions in the case of the SIIR model).

###  3.3 Critical Surfaces and Mutations

The assumption of a deformation away from the fixed point $P_{0}$ in Figure 13
(which represents the complete absence of the disease in the population) along
a generic direction is not realistic, since it would correspond to the
simultaneous appearance of infected individuals of both variants in the
population. A more likely scenario would be the appearance of one variant
first, while a second deformation at a latter stage introduces the second
variant. This dynamics can be understood from an RG perspective as the
switching on of a relevant operator.

To make this statement more precise, we first need to introduce the concept of
_critical surface_ associated with a fixed point of the $\beta$-functions. A
critical surface consists in all points in the $\mathbb{P}$-plane from where
the RG-flow leads to the fixed point in question. Concretely, for the fixed
points $P_{1,2,3}$, the critical surfaces are

$\displaystyle\mathcal{C}_{P_{1}}=\\{(I_{\text{c},1},I_{\text{c},2})\in\mathbb{P}|I_{\text{c},1}>0\text{
and }I_{\text{c},2}=0\\}\,,$
$\displaystyle\mathcal{C}_{P_{2}}=\\{(I_{\text{c},1},I_{\text{c},2})\in\mathbb{P}|I_{\text{c},1}=0\text{
and }I_{\text{c},2}>0\\}\,,$
$\displaystyle\mathcal{C}_{P_{3}}=\\{(I_{\text{c},1},I_{\text{c},2})\in\mathbb{P}|I_{\text{c},1}>0\text{
and }I_{\text{c},2}>0\\}\,.$ (3.10)

A _relevant operator_ (from the perspective of the fixed point in question),
corresponds to a direction that drives the theory away from the critical
surface, such that it flows to a new critical point. In the case at hand,
$\mathcal{C}_{P_{1,2}}$ have one critical direction orthogonal to it, which,
from an epidemiological perspective, precisely corresponds to the appearance
of the second variant. A small deformation at any point of
$\mathcal{C}_{P_{1,2}}$ (for example due to a relevant mutation of the virus)
causes the system to deviate from the critical surface and ultimately flow
towards $P_{3}$.

In a scenario with only two variants, the fixed point $P_{3}$ has no relevant
deformations and is attractive along all directions. Instead, small
fluctuations along trajectories leading towards $P_{3}$ (such as the red path
shown in Figure 13) can be interpreted as _irrelevant_ operators being
switched on. To make this statement more concrete, instead of
$(I_{\text{c},1},I_{\text{c},2})$ we can consider the following
functions111111SH would like to thank Michele Della Morte for useful exchanges
on the form of these functions as $O(2)$ transformations of
$(I_{\text{c},1},I_{\text{c},2})$.

$\displaystyle\mathcal{O}_{+}=\frac{A_{1}A_{2}}{A_{1}^{2}+A_{2}^{2}}\left(\frac{I_{\text{c},1}}{NA_{2}}+\frac{I_{\text{c},2}}{NA_{1}}\right)\,,$
and
$\displaystyle\mathcal{O}_{-}=\frac{A_{1}A_{2}}{A_{1}^{2}+A_{2}^{2}}\left(-\frac{I_{\text{c},1}}{NA_{1}}+\frac{I_{\text{c},2}}{NA_{2}}\right)\,.$
(3.11)

Notice that, at the fixed point $P_{3}$, we have

$\displaystyle\mathcal{O}_{-}(I_{\text{c},1}=NA_{1},I_{\text{c},2}=NA_{2})=0\,.$
(3.12)

Using the explicit solutions of $(I_{\text{c},1},I_{\text{c},2})(t)$ in (3.7)
we have schematically plotted $\mathcal{O}_{\pm}$ and its $\beta$-functions as
functions of time in Figure 14. The latter can be written in the form

$\displaystyle-\beta_{\pm}=\frac{d\mathcal{O}_{\pm}}{dt}=\frac{1}{N\,(A_{1}^{2}+A_{2}^{2})}\,\nabla_{\pm}\,\Phi(\mathcal{O}_{\pm})\,,\text{with}$
$\displaystyle\nabla_{\pm}=\frac{\partial}{\partial\mathcal{O}_{\pm}}\,.$
(3.13)

We can then interpret the time evolution of the system as a flow of
$\mathcal{O}_{+}$ from the repulsive fixed point $\mathcal{O}_{+}=0$ to the
attractive fixed point $\mathcal{O}_{+}=1$, while $\mathcal{O}_{-}$ can be
interpreted as an irrelevant operator that is switched on along the way, which
does not drive the system away from the fixed point
$(\mathcal{O}_{+},\mathcal{O}_{-})=(1,0)$.

Figure 14: Parametric Plot of the operators $\mathcal{O}_{\pm}$ (left panel)
and their $\beta$-functions (right panel) as functions of time. The plots use
$N\lambda_{1}=0.2$, $N\lambda_{2}=0.25$, $A_{1}=0.3$, $A_{2}=0.4$ and
$\kappa_{1}=16$, $\kappa_{2}=25$.

###  3.4 Multi Wave Dynamics and CeRG

We now consider in more detail a particular case of a flow along the critical
surface $\mathcal{C}_{P_{1}}$ with a relevant operator being switched on along
the way (_i.e._ the second variant appearing at some moment $t_{0}>0$) through
a small fluctuation. If $\sigma_{2}$ is not too large compared to
$\sigma_{1}$, for some time after the relevant deformation is turned on, the
theory stays close to the (initial) critical surface $\mathcal{C}_{P_{1}}$ and
flows towards the fixed point $P_{1}$. At some later time, the relevant
deformation along direction $I_{\text{c},2}$ becomes too large and the flow
runs significantly away from the critical surface $\mathcal{C}_{P_{1}}$. If
the deformation appears well after $I_{1}(t)$ has reached a maximum, the flow
can be described as a _crossover_ flow (see Figure 15 for a schematic
example): the RG-flow can be decomposed into a flow along the original
critical surface $\mathcal{C}_{P_{1}}$ followed by a flow perpendicular to it.
The latter drives the system from the proximity of the fixed point $P_{1}$ to
the new one $P_{3}$. From the perspective of the new fixed point, the second
phase of the flow looks like an RG flow from a UV fixed point $P_{1}$ to an IR
one $P_{3}$.

Figure 15: Left panel: Schematic example of a crossover flow. During the first
part, the system follows closely the critical surface $\mathcal{C}_{P_{1}}$
(_i.e._ parallel to the $I_{\text{c},1}$-axis) and comes close to the fixed
point $P_{1}$ without reaching it. After staying for some time in the vicinity
of $P_{1}$ the system enters into the second part of the flow to the fixed
point $P_{3}$. Right panel: Cumulative numbers of infected as functions of
time. The total number of infected shows a two-wave structure with a linear
growth phase at around $t=30$ corresponding to the time that the RG flow is in
proximity to $P_{1}$. The numerical values used for the plots are
$N\lambda_{1}=0.2$, $N\lambda_{2}=0.25$ and $A_{1}=0.3$, $A_{2}=0.4$.

Since during the first part of the flow, the number of active infected of the
second variant $I_{2}(t)$ is fairly small, this flow is very well approximated
by a usual eRG dynamics (see eq. (3.5)), which has been shown [27] to describe
real world data very well. Once the system reaches the vicinity of the fixed
point $P_{1}$, it will then enter into a quasi-linear growth phase, in which
the number of active infected with respect to both variants is small and
therefore the total number of infected $(I_{\text{c},1}+I_{\text{c},2})(t)$
only grows linearly (see Figure 15). However, after a certain time, the number
of infected $I_{\text{c},2}$ will grow exponentially (while the number of
infected with respect to the original variant remains small) and the system
enters into the crossover phase. Now, the $\beta$-function for
$I_{\text{c},2}$ is, once again, essentially modelled by a standard eRG
equation (see eq. (3.5)) describing the flow to $P_{3}$.

In this picture, the two-wave structure is explained as the (more or less
successive) appearance of two different variants of the disease. In
particular, the linear-growth phase (that has for example been observed in
real world data in the inter-wave period of the COVID-19 pandemic [35]) is
explained by the fact that the system comes close to a fixed point, which,
however, it cannot reach. It nevertheless spends significant time in its
proximity. A similar reasoning underlies the _Complex epidemic Renormalisation
Group_ (CeRG) approach [34]. The CeRG $\beta$-function of the following type
was proposed for the combined number of cumulative infected
$I_{\text{c},\text{tot}}=I_{\text{c},1}+I_{\text{c},2}$

$\displaystyle-\beta_{\text{CeRG}}=\frac{dI_{\text{c},\text{tot}}}{dt}(t)=\lambda\,I_{\text{c},\text{tot}}(t)\,\left[\left(1-\zeta\,\frac{I_{\text{c},\text{tot}}(t)}{A}\right)^{2}-\delta\right]^{p_{1}}\,\left(1-\frac{I_{\text{c},\text{tot}}(t)}{A}\right)^{p_{2}}\,,$
(3.14)

with $A$ the asymptotic number of infected, $\lambda$ the infection rate,
$\zeta>1$ and $\delta<0$. Indeed, besides the fixed points
$I_{\text{c},\text{tot}}=0$ and $I_{\text{c},\text{tot}}=A$, the beta-function
(3.14) also has the complex fixed points
$I_{\text{c},\text{tot}}=\frac{A}{\zeta}\left(1\pm i\sqrt{|\delta|}\right)$,
which cannot be reached by the flow, but are responsible for the linear-growth
phase.

To compare (3.5) and (3.6) with the $\beta$-function in (3.14), we assume that
$\sigma_{1}$ and $\sigma_{2}$ are not significantly different from one another
and that the mutation occurs significantly after the maximum number of
infected of the first variant (such that the condition of a crossover flow is
satisfied). Furthermore, we can write the following combined $\beta$-function
for the total cumulative number of infected

$\displaystyle-\beta_{\text{tot}}=$
$\displaystyle\frac{1}{N}\left[\frac{dI_{\text{c},1}}{dt}(t)+\frac{dI_{\text{c},2}}{dt}(t)\right]\sim\theta\left(\left(NA_{1}-I_{\text{c},1}\right)\right)\lambda_{1}\,I_{\text{c},1}\left(1-\frac{I_{\text{c},1}}{A_{1}}\right)$
$\displaystyle+\theta\left(\left(I_{\text{c},2}-\xi\right)\left(N(A_{1}+A_{2})-I_{\text{c},2}\right)\right)\lambda_{2}\,\left(I_{\text{c},2}-N\,A_{1}\right)\left(1-\frac{I_{\text{c},2}-NA_{1}}{NA_{2}}\right)\,,$
(3.15)

where $\theta$ is the Heaviside step-function. Furthermore, we can use the
solutions (3.7) of $(I_{\text{c},1},I_{\text{c},2})(t)$ as logistic functions
to schematically plot (3.15): the left panel of Figure 16 shows a parametric
plot of
$(I_{\text{c},1}(t)+I_{\text{c},2}(t),\tfrac{dI_{\text{c},1}+I_{\text{c},2}}{dt}(t))$
for different values of $t$. The latter are very well approximated by (3.15)
(shown by the thin black line), except for a small region around
$I_{\text{c},t}\sim A_{1}$, in which the beta-function does not in fact reach
zero, but interpolates between the two terms in (3.15). This region
corresponds to a non-trivial interaction between the variants and governs the
transition from the first part of the flow (close to the original critical
surface) to the crossover flow. It precisely corresponds to the linear growth
region in the context of the CeRG: Indeed, a similar shape of the beta-
function can also be achieved through a function of the form (3.14), as is
shown in the right panel of (16). The region around
$I_{\text{c},\text{tot}}\sim A_{1}$ corresponds to the RG-flow not quite
reaching a zero, thus leading to the quasi-linear growth phase.

Figure 16: Schematic plot of the beta-function of the total number of infected
computed from the solutions (3.7) (red dots). Left panel: comparison to
(3.15), right panel: fitting with a $\beta$-function of the type (3.14), with
$\kappa=0.2075$, $\delta=-0.0021$, $p_{1}=0.336$, $p_{2}=0.959$, $A=0.6996$
and $\frac{A}{\zeta}=0.2995$.

##  4 Explicit Examples from the COVID-19 pandemic

In the following we present selected examples of number of infected
individuals for different countries during the COVID-19 pandemic. These serve
merely as illustration for some of the points we have outlined in the previous
sections.

###  4.1 Second wave in California

As a first example, we consider the evolution of the SARS-CoV-2
Epsilon121212This variant comprises the variants lineage B.1.429 and B.1.427
under the Phylogenetic Assignment of Named Global Outbreak Lineages (pangolin)
tool and was first detected in California in July 2020. variant relative to
the other variants in California. The daily number of new cases is shown in
Figure 17 and we have also indicated changes in the social distancing and
lockdown measures imposed at the same time. To obtain the number of new
infections with the Epsilon variant relative to the others, we have used data
from GISAID [38] that provide the sequenced genomes of samples taken in
California from Sep/2020 until Apr/2021 (we refer to [1] for more details on
our methodology): while only a fraction of all positive test samples per day
is genetically analysed, we assume that the distribution of the Epsilon
variant in this subset is representative of the distribution of the variant
among all infected individuals in California. We therefore have first
calculated the percentage of each variant among the sequences analysed at each
specific date. By multiplying this percentage with the total number of new
cases (obtained from [39, 40]) within all of California, we are able to
extract an approximation of the number of new cases per day for each variant.
The statistical uncertainty inherent in this procedure has been estimated and
taken into account when fitting the cumulative number of cases (see Figure 18
below). Due to the large number of tests and genome sequencings performed in
California and the rather short duration of our study, the resulting
uncertainty is moderate. Indeed, the maximal number of new infections for both
curves in Figure 17 occur during a period of 2-3 months (Nov/2020 until
Jan/2021), during which furthermore the regional lockdown measures have stayed
unchanged. Therefore, in order to model the second wave of the COVID-19
pandemic in California, our implicit assumption that the parameters of the eRG
model (_i.e._ $\lambda_{1,2}$ and $A_{1,2}$) are time-independent seems a
reasonable starting point. Indeed, from the data in Figure 17 we can calculate
the cumulative number of infected individuals, which is shown in Figure 18
along with an approximation in terms of logistic functions (3.7) with constant
coefficients. Denoting the cumulative number of individuals infected with the
Epsilon variant and the other ones by $I_{\text{c},2}$ and $I_{\text{c},1}$
respectively, the parameters extracted from the fit are given by

Figure 17: Number of daily new infections in California as a function of time
since 08/Sep/2020. The epidemiological data have been extracted from [39, 40],
while the number of cases for the Epsilon variant are based on genome
sequencing data extracted from [38]. Finally, the information about lockdown
and social distancing measures is taken from [41]. Figure 18: Cumulative
number of infected individuals with respect to the different variants in
California since 08/Sep/2020 and their approximation in terms of logistic
functions. The latter have in fact been optimised to fit the daily new
infections taking into account their statistical errors.

$\displaystyle NA_{1}=(1.644\pm 0.021)\cdot 10^{6}\,,$
$\displaystyle\kappa_{1}=105.92\pm 0.40\,,$
$\displaystyle\lambda_{1}=0.0590\pm 0.0009\,,$ $\displaystyle NA_{2}=(1.090\pm
0.014)\cdot 10^{6}\,,$ $\displaystyle\kappa_{2}=123.00\pm 0.32\,,$
$\displaystyle\lambda_{2}=0.0799\pm 0.0007\,.$ (4.1)

These numbers are in fact determined by fitting the derivative of
$I_{\text{c},1,2}$ (_i.e._ the daily new cases), as shown in Figure 19. With
the help of these approximations, we can compute the function $\Phi$ in (3.6)
which in turn allows us to compute the $\beta$-functions (3.5). The
corresponding flow of the system in the
$(I_{\text{c},1},I_{\text{c},2})$-plane is shown in the left panel of Figure
20, where we have chosen the scheme $f_{i}(I_{\text{c},i})=I_{\text{c},i}$.
The right panel of that figure gives the relative deviation of the
$\beta$-functions from the time derivative of
$(I_{\text{c},1},I_{\text{c},2})$ as given by the number of daily new cases.
The plot also indicates changes imposed by the government on the social
distancing measures among the population, leading to a stronger deviation of
the modelled $\beta$-functions from the actual data.

Figure 19: Comparison of the time-dependence of the $\beta$-functions (3.5)
for $(I_{\text{c},1},I_{\text{c},2})$ with the numbers of daily new
infections.

Figure 20: Left panel: flow of the system in the
$(I_{\text{c},1},I_{\text{c},2})$-plane. The vectors indicate the
$\beta$-function as given by (3.5) and (3.6) (with
$f_{i}(I_{\text{c},i})=I_{\text{c},i}$ for $i=1,2$). Right panel: relative
difference of the modelled $\beta$-function (3.5) with the actual number of
daily infected along the flow of the system in the
$(I_{\text{c},1},I_{\text{c},2})$-plane.

From the left panel of Figure 20 we can see that the flow of the system for
the most part is not close to either of the two axes, and thus has not the
markings of a crossover flow. This is also clear from the fact that the maxima
of daily new infections of the variant B.1.427 and the remaining ones are not
very well separated in time (see Figure 17). Nevertheless, we expect that the
formulation of the flow in terms of the relevant operator $\mathcal{O}_{+}$
and the irrelevant operator $\mathcal{O}_{-}$ defined in Eq. (3.11) may still
give a viable description. The time dependence of $\mathcal{O}_{\pm}$ is shown
in the left panel of Figure 21, while the right panel of the same figure shows
an approximation of the corresponding $\beta$-functions. Finally, the flow of
the system in the $(\mathcal{O}_{+},\mathcal{O}_{-})$ plane (along with the
gradient vector field (3.13)) is shown in Figure 22.

Figure 21: Left panel: The functions $\mathcal{O}_{\pm}$ defined in eq. (3.11)
computed from the actual numbers of infected individuals as a function of
time. Right panel: The corresponding $\beta$-functions together with their
approximations (dashed black lines) implied by (4.1). Figure 22: Flow of the
system in the $(\mathcal{O}_{+},\mathcal{O}_{-})$-plane based on the data for
California.

We remark that the main reason for deviation of the flow from a gradient flow
(as is showcased in the right panel of Figure 20) can be attributed to the
fact that the active number of infected individuals is not zero at the
beginning and the end of the flow (see Figure 19). This is, in fact, because
we are not describing the entire pandemic in California, but only the period
of September 2020 to April 2021, when a second wave hit the state. Therefore,
we are not really describing the flow from one fixed point to another, but
rather the flow between two linear-growth phases of the system. As such, the
description in terms of the $\beta$-function (3.5) is only approximative,
which is particularly visible in the beginning and the end of the flow.

###  4.2 Second and third waves in the UK

As a next example, we consider the time evolution of the SARS-CoV-2
Alpha131313This variant is also called lineage B.1.1.7 under the pangolin tool
and was first found in Nov. 2020 (in a sample dating from September 2020) in
the UK. variant in the United Kingdom. The daily number of new cases starting
from 01/July/2020 is shown in Figure 23, where we have also indicated changes
in the social distancing and lockdown measures. The number of new infections
with the Alpha variant have been extracted using GISAID [38] combined with
epidemiological data from [39]. Due to the high number of PCR tests and genome
analysis performed in the UK, the inherent statistical uncertainty is rather
small (see _e.g._ Figure 24). In contrast to the evolution (of a single wave)
in California in the previous subsection, the time evolution studied in this
example spans two distinct waves lasting roughly 5 months (Oct/2020 until
Feb/2021). As Figure 23 indicates, during this period141414We also remark that
the time period includes the Christmas holidays, which traditionally leads to
increased social activity and travel among the population, the lockdown
measures have not remained constant, which potentially leads to additional
effects, as we shall remark later on. From Figure 23 we can compute the
cumulative number of infected individuals, which is shown in Figure 24. The
latter also shows approximations in terms of logistic functions: denoting the
cumulative number of individuals infected with the Alpha variant and the other
ones by $I_{\text{c},2}$ and $I_{\text{c},1}$ respectively, the parameters
extracted from the fit are given by

Figure 23: Number of daily new infections in the United Kingdom as a function
of time since 01/July/2020. The epidemiological data have been extracted from
[39], while the number of cases for the Alpha variant are based on genome
sequencing data extracted from [38]. Finally, the information about lockdown
and social distancing measures is taken from [42]. Figure 24: Cumulative
number of infected individuals with respect to the different variants in the
United Kingdom since 01/July/2020 and their approximation in terms of logistic
functions. The latter have in fact been optimised to fit the daily new
infections.

$\displaystyle NA_{1}=(1.895\pm 0.035)\cdot 10^{6}\,,$
$\displaystyle\kappa_{1}=136.06\pm 0.78\,,$
$\displaystyle\lambda_{1}=0.0447\pm 0.0009\,,$ $\displaystyle NA_{2}=(2.007\pm
0.030)\cdot 10^{6}\,,$ $\displaystyle\kappa_{2}=197.12\pm 0.33\,,$
$\displaystyle\lambda_{2}=0.0812\pm 0.0014\,.$ (4.2)

As in the case of California, these numbers are in fact determined by fitting
the derivative of $I_{\text{c},1,2}$ (_i.e._ the daily new cases), as
indicated in Figure 25. While the fit of $\beta_{2}$ correctly captures a
single peak, the fit of $\beta_{1}$ gives a single maximum rather than two.
The appearance of the second maximum in the red curve in the left part of
Figure 25 (stemming from the maximum of the orange curve in Figure 23 in the
end of Dec/2020) cannot be explained with statistical uncertainties inherent
in the way we have extracted the relative number of new infections with the
Alpha variant and may have different reasons:

  * _(i)_

change of the infection rate of all variants due to modified social behaviour
and/or travelling habits related to the Christmas holidays

  * _(ii)_

local geographic effects not captured by the data

  * _(iii)_

appearance of additional (subdominant) variants of the virus

  * _(iv)_

non-trivial interaction of the virus variants that are not captured by the
beta-functions (3.5)

Since the effect is rather small, we shall continue with the approximations in
Figure 25 and leave the analysis of this effect (notably the possibility
_(iv)_) for future work. With these approximations we can compute the function
$\Phi$ in (3.6) which in turn determines the $\beta$-functions (3.5). The
corresponding flow of the system in the
$(I_{\text{c},1},I_{\text{c},2})$-plane is shown in the left panel of Figure
26, where we have chosen the scheme $f_{i}(I_{\text{c},i})=I_{\text{c},i}$.
The right panel of that figure gives the relative deviation of the
$\beta$-functions from the time derivative of
$(I_{\text{c},1},I_{\text{c},2})$ as given by the number of daily new cases.

Figure 25: Comparison of the time-dependence of the $\beta$-functions (3.5)
for $(I_{\text{c},1},I_{\text{c},2})$ with the numbers of daily new infections
in the UK.

Figure 26: Left panel: Vector field and flow of the system in the
$(I_{\text{c},1},I_{\text{c},2})$-plane. Right panel: relative difference of
the modelled $\beta$-function (3.5) with the actual number of daily infected
along the flow of the system in the $(I_{\text{c},1},I_{\text{c},2})$-plane
for the UK. As mentioned before, the deviations from the beta-function
predicted by the eRG approach may be related to the changes in the lockdown
measures or may indicate of additional subleading effects not captured in the
current approach.

From the flow diagram in Figure 26 we can see that the system remains close to
the $I_{\text{c},1}$-axis. We can therefore try to see whether it can be
parametrised in a fashion resembling a crossover flow and whether it is
possible to explain the two-wave structure from the flow close to a fixed
point. To this end, we have plotted the time derivative of the total number of
infected individuals $I_{\text{c},\text{tot}}=I_{\text{c},1}+I_{\text{c},2}$
in Figure 27 along with approximations along the lines of eq. (3.15) and
(3.14) (along with its $0.999$ confidence interval). Indeed, we can see that
the function has a pronounced local minimum, which models the proximity of the
system to the fixed point near the $I_{\text{c},1}$-axis and which is
responsible for the short linear growth phase in the end of November/beginning
of December, as can be seen in Figure 23.

Figure 27: Approximation of the time derivative of
$I_{\text{c},\text{tot}}=I_{\text{c},1}=I_{\text{c},2}$ (red dots). Left
panel: $\beta_{\text{tot}}$ as defined in (3.15); Right panel:
$-\beta_{\text{CeRG}}$ as in eq. (3.14).

###  4.3 First and second waves in South Africa

As a final example, we consider the time evolution of the SARS-CoV-2
Beta151515This variant is also called lineage B.1.351 under the pangolin tool
and was first found in the Eastern Cape province of South Africa. variant in
South Africa. The daily number of new cases starting from 08/March/2020 is
shown in Figure 28. The number of new infections with the Beta variant have
been extracted using GISAID [38] combined with epidemiological data from [39].
Due to the rather small number of genome sequences of samples, the separation
between the Beta variant and others is afflicted with a rather large (time-
dependent) statistical uncertainty, which needs to be taken into account in
the following and which makes the interpretation of certain results delicate.
In Figure 29 we have plotted the cumulative number of infected individuals
alongside with their uncertainties. This figure also shows an approximation in
terms of logistic functions, which is weighted by the (time-dependent)
uncertainties: denoting the cumulative number of individuals infected with the
variant B.1.351 by $I_{\text{c},2}$ and the cumulative number of infected with
the other invariants by $I_{\text{c},1}$, the fit parameters along eq.(3.7)
are given by

Figure 28: Number of daily new infections in South Africa as a function of
time since 08/March/2020. The epidemiological data have been extracted from
[39], while the number of cases for the Beta variant are based on genome
sequencing data extracted from [38]. Finally, the information about lockdown
and social distancing measures is taken from [43].

Figure 29: Cumulative number of infected individuals for the variant B.1.351
(right panel) and all other variants (left panel) since 08/March/2020 and
their approximation in terms of logistic functions. The error bars represent
the (accumulated) statistical uncertainty following the rather low number of
genome sequencings compared to the number of (positive) tests.

$\displaystyle NA_{1}=(0.636\pm 0.009)\cdot 10^{6}\,,$
$\displaystyle\kappa_{1}=131.86\pm 0.35\,,$
$\displaystyle\lambda_{1}=0.0739\pm 0.0012\,,$ $\displaystyle NA_{2}=(0.852\pm
0.020)\cdot 10^{6}\,,$ $\displaystyle\kappa_{2}=298.91\pm 1.08\,,$
$\displaystyle\lambda_{2}=0.0424\pm 0.0007\,.$ (4.3)

These numbers are determined by fitting (part of) the derivative of
$I_{\text{c},1,2}$ (_i.e._ the daily new cases), as indicated in Figure 30.
With these approximations we may compute the function $\Phi$ in (3.6) which in
turn gives an approximation of the $\beta$-functions (3.5). The flow of the
system following the latter in the $(I_{\text{c},1},I_{\text{c},2})$-plane is
shown in Figure 31 (as before, we use the scheme
$f_{i}(I_{\text{c},i})=I_{\text{c},i}$). In the right panel, for better
visibility, we have combined the error bars of the (black) data points into a
gray region: while the black trajectory is not actually hitting the predicted
fixed point of the beta-function, the latter is within the error bars. It is
therefore difficult to say, whether this indicates additional subleading
effects in the time-evolution of the virus, or merely statistical uncertainty.
Apart from this effect, the flow strongly resembles a crossover flow, which
for the first part follows the $I_{\text{c},1}$-axis, coming close to a fixed
point at $(NA_{1},0)$. However, the appearance of the variant B.1.351 triggers
a crossover flow to a new fixed point.

Figure 30: Comparison of the time-dependence of the $\beta$-functions (3.5)
for $(I_{\text{c},1},I_{\text{c},2})$ with the numbers of daily new infections
in South Africa. The fits (along with their 0.99 confidence interval) are
weighted with respect to the statistical uncertainty of the data. This
explains for example, why the maximum of the red curve around day 300 in the
left panel is effectively neglected: it is afflicted with a very high
statistical uncertainty.

Figure 31: Vector field and flow of the system in the
$(I_{\text{c},1},I_{\text{c},2})$-plane. Left panel: the black dots represent
the actual numbers of cumulative infected (together with their error bars),
while the red dots represent the predicted eRG flow which lies within the
error bars. right panel: For better visibility we have fused the error bars
into the grey zone representing the statistical error of the cumulative
numbers of infected individuals.

The fact that the flow in Figure 31 stays close to the $I_{\text{c},1}$-axis
during the first part of the flow, allows us to model it in the form of a
crossover flow. To this end, we have plotted the time derivative of the total
number of infected individuals
$I_{\text{c},\text{tot}}=I_{\text{c},1}+I_{\text{c},2}$ in Figure 32 along
with approximations along the lines of eq. (3.15) and (3.14). Indeed, we can
see that the function has a pronounced local minimum, which models the
proximity of the system to the fixed point near the $I_{\text{c},1}$-axis and
which is responsible for the linear growth phase between the end of August
2020 and the end of November 2020 as can be seen in Figure 28.

Figure 32: Approximation of the time derivative of
$I_{\text{c},\text{tot}}=I_{\text{c},1}=I_{\text{c},2}$ (red dots). Left
panel: $\beta_{\text{tot}}$ as defined in (3.15); Right panel:
$-\beta_{\text{CeRG}}$ as in eq. (3.14) together with its 0.99 confidence
interval.

##  5 Conclusions

In this paper we have extended the epidemic Renormalisation Group (eRG)
approach to study the time evolution of several competing variants of a
disease. As a first step, in order to gain intuition, we have analysed a
simple compartmental model (termed SIIR) with two different groups of
infectious individuals (and different infection and recovery rates). Numerical
solutions of the model indicate that as long as the reproduction numbers of
the two variants are not too different, the cumulative number of infected
individuals of the two variants can be well approximated by independent
logistic functions (sigmoids). Moreover, we have approximated the dynamics in
terms of flow equations that describe the trajectories of the system in the
$\mathbb{P}$-plane keeping track of the cumulative number of infected. The
resulting equations (2.29) and (2.35) can be compactly formulated in terms of
the gradient of a single function $\Phi$, which is quadratic for variants with
a reproduction number $\sigma<1$ and cubic for $\sigma>1$. Furthermore, from a
theoretical perspective, we find it striking that the endpoints of the flows
are characterised by co-dimension 1 surfaces in the $\mathbb{P}$-plane that
represent equivalent system from the perspective of the SIIR model, akin to
surfaces of fixed points that are related by the action of marginal operators
in the context of conformal field theories.

As a second step, we have used the intuition gained from the simple SIIR model
to propose a generalisation of the eRG framework to include the dynamics of
multiple variants. Defining a beta-function (3.5) in the form of a gradient
equation, we have analysed its fixed points in the $\mathbb{P}$-plane, as well
as different trajectories connecting them. In particular, we have made contact
to the CeRG approach [34, 35] which has modelled the multi-wave structure of
epidemics with the help of complex fixed points: in this regard, quasi-linear
growth phases separating two waves are explained by the system coming close to
a complex fixed point, which it cannot reach. In the current paper we have
shown that such a behaviour can occur naturally through the appearance of a
new variant of the disease.

Finally, we have confronted our model with data from the spread of different
variants of SARS-CoV-2 in California, the United Kingdom and South Africa,
thus empirically validating our approach.

In the future it will be interesting to further generalise and extend the
approach developed here: on the one hand going beyond two competing variants
will lead to a richer structure of fixed points for the system, thus allowing
to model more complex multi-wave pandemics. On the other hand, so far our
analysis has not taken into account other factors that govern the time
evolution of pandemics, such as the impact of vaccines and non-pharmaceutical
interventions but also the possibility for re-infections (through only partial
immunity granted from recovery from a given variant). We hope to be able to
return to these points in the future. Finally, we shall use the model
developed in this paper for a computer-aided analysis of the spread of
different SARS-CoV-2 variants in Europe and the USA along the lines of the
work in our companion paper [1].

## References

  * [1] G. Cacciapaglia, F. Conventi, C. Cot, A. de Hoffer, A. Giannini, S. Hohenegger, F. Sannino, and S. Vatani, “Evolutionary model of COVID-19 diffusion via mutations applying Machine Learning to spike protein sequencing,” to appear.
  * [2] W. Hamer, “Age-incidence in relation with cycles of disease prevalence,” Trans. Epidem. Soc. London, vol. 15, pp. 64–77, 1896.
  * [3] W. Hamer, “Epidemic disease in England: The evidence of variability and of persistency of type; Lecture 1,” Lancet, pp. 569–574, March 1906.
  * [4] W. Hamer, “Epidemic disease in England: The evidence of variability and of persistency of type; Lecture 2,” Lancet, pp. 655–662, March 1906.
  * [5] W. Hamer, “Epidemic disease in England: The evidence of variability and of persistency of type; Lecture 3,” Lancet, pp. 733–739, March 1906.
  * [6] R. Ross, “The Prevention of Malaria,” second edition, John Murray, London, 1911.
  * [7] R. Ross, “An application of the theory of probabilities to the study of _a priori_ pathometry: Part I,” Proc. Roy. Soc. Lond. A, vol. 92, pp. 204–230, 1916.
  * [8] R. Ross and H. Hudson, “An application of the theory of probabilities to the study of _a priori_ pathometry: Part II,” Proc. Roy. Soc. Lond. A, vol. 93, pp. 212–225, 1916.
  * [9] R. Ross and H. Hudson, “An application of the theory of probabilities to the study of _a priori_ pathometry: Part III,” Proc. Roy. Soc. Lond. A, vol. 93, pp. 225–240, 1916.
  * [10] A. McKendrick, “The rise and fall of epidemics,” Paludism (Transactions of the Committee for the Study of Malaria in India), vol. 1, pp. 54–66, 1912.
  * [11] A. McKendrick, “Studies on the theory of continuous probabilities, with special reference to its bearing on natural phenomena of a progressive nature,” Proceedings of the London Mathematical Society, vol. 13, pp. 401–416, 1914.
  * [12] A. McKendrick, “Applications of mathematics to medical problems,” Proc. Edinburgh Math. Soc., vol. 44, pp. 98–130, 1926.
  * [13] W. O. Kermack, A. McKendrick, and G. T. Walker, “A contribution to the mathematical theory of epidemics,” Proceedings of the Royal Society A, vol. 115, pp. 700–721, 1927.
  * [14] M. Perc, J. J. Jordan, D. G. Rand, Z. Wang, S. Boccaletti, and A. Szolnoki, “Statistical physics of human cooperation,” Physics Reports, vol. 687, pp. 1 – 51, 2017.
  * [15] Z. Wang, M. A. Andrews, Z.-X. Wu, L. Wang, and C. T. Bauch, “Coupled disease–behavior dynamics on complex networks: A review,” Physics of Life Reviews, vol. 15, pp. 1 – 29, 2015.
  * [16] Z. Wang, C. T. Bauch, S. Bhattacharyya, A. d’Onofrio, P. Manfredi, M. Perc, N. Perra, M. Salathé, and D. Zhao, “Statistical physics of vaccination,” Physics Reports, vol. 664, pp. 1 – 113, 2016.
  * [17] H. W. Hethcote, “The mathematics of infectious diseases,” SIAM Review, vol. 42, no. 4, 2000.
  * [18] N. Bailey, The Mathematical Theory of Infectious Diseases, 2nd ed. Hafner, New York, 1975.
  * [19] J. W. Essam, “Percolation theory,” Rep. Prog. Phys., vol. 43, p. 833, 1980\.
  * [20] D. Stauffer, “Scaling theory of percolation clusters,” Phys. Rep. , vol. 54, pp. 1–74, 1979.
  * [21] C. Domb, “Fluctuation phenomena and stochastic processes,” Nature, vol. 184, pp. 509–12, 1959.
  * [22] L. Peliti, “Path integral approach to birth-death processes on a lattice,” J. Phys. France (Paris), vol. 46, pp. 1469–1483, 1985.
  * [23] M. Doi, “Second quantization representation for classical many-particle system,” J. Phys. A: Math. Gen., vol. 9, p. 1465, 1976.
  * [24] M. Doi, “Stochastic theory of diffusion-controlled reaction,” J. Phys. A: Math. Gen., vol. 9, p. 1479, 1976.
  * [25] J. L. Cardy and P. Grassberger, “Epidemic models and percolation,” Journal of Physics A: Mathematical and General, vol. 18, pp. L267–L271, apr 1985\.
  * [26] P. Grassberger, “On the critical behavior of the general epidemic process and dynamical percolation,” Mathematical Biosciences, vol. 63, no. 2, pp. 157 – 172, 1983.
  * [27] M. Della Morte, D. Orlando, and F. Sannino, “Renormalization Group Approach to Pandemics: The COVID-19 Case,” Front. in Phys., vol. 8, p. 144, 2020\.
  * [28] K. G. Wilson, “Renormalization group and critical phenomena. 1. Renormalization group and the Kadanoff scaling picture,” Phys. Rev. B, vol. 4, pp. 3174–3183, 1971.
  * [29] K. G. Wilson, “Renormalization group and critical phenomena. 2. Phase space cell analysis of critical behavior,” Phys. Rev. B, vol. 4, pp. 3184–3205, 1971.
  * [30] G. Cacciapaglia and F. Sannino, “Interplay of social distancing and border restrictions for pandemics (COVID-19) via the epidemic Renormalisation Group framework,” Sci Rep, vol. 10, p. 15828, 5 2020.
  * [31] G. Cacciapaglia, C. Cot, and F. Sannino, “Mining google and apple mobility data: Temporal anatomy for covid-19 social distancing,” Scientific Reports, vol. 11, p. 4150, 2021.
  * [32] G. Cacciapaglia, C. Cot, A. S. Islind, M. Óskarsdóttir, and F. Sannino, “Impact of us vaccination strategy on covid-19 wave dynamics,” Scientific Reports, vol. 11(1), pp. 1–11, 2021.
  * [33] G. Cacciapaglia, C. Cot, and F. Sannino, “Second wave covid-19 pandemics in europe: A temporal playbook,” Sci Rep, vol. 10, p. 15514, 2020.
  * [34] G. Cacciapaglia and F. Sannino, “Evidence for complex fixed points in pandemic data,” Front. Appl. Math. Stat., vol. 7, p. 659580, 2021.
  * [35] G. Cacciapaglia, C. Cot, and F. Sannino, “Multiwave pandemic dynamics explained: How to tame the next wave of infectious diseases,” Scientific Reports, vol. 11, p. 6638, 2021.
  * [36] G. Cacciapaglia, C. Cot, M. Della Morte, S. Hohenegger, F. Sannino, and S. Vatani, “The field theoretical ABC of epidemic dynamics,” 1 2021.
  * [37] M. Della Morte and F. Sannino, “Renormalisation Group approach to pandemics as a time-dependent SIR model,” Front. in Phys., vol. 8, p. 583, 2021\.
  * [38] “Gisaid.” https://www.gisaid.org.
  * [39] “Our world in data.” https://ourworldindata.org.
  * [40] “New york times github.” https://github.com/nytimes/covid-19-data/blob/master/us-states.csv.
  * [41] “Covid-19 pandemic in california.” https://en.wikipedia.org/wiki/COVID-19_pandemic_in_California.
  * [42] “Covid-19 pandemic in the united kingdom.” https://en.wikipedia.org/wiki/COVID-19_pandemic_in_the_United_Kingdom.
  * [43] “Covid-19 pandemic in south africa.” https://en.wikipedia.org/wiki/COVID-19_pandemic_in_South_Africa.