# Agent-based modeling of COVID-19 outbreaks for New York state, UK and
Novosibirsk region

Olga Krivorotko<EMAIL_ADDRESS>Mariia Sosnovskaia
<EMAIL_ADDRESS>Ivan Vashchenko<EMAIL_ADDRESS>Cliff Kerr
<EMAIL_ADDRESS>Daniel Lesnic<EMAIL_ADDRESS>Institute of Computational
Mathematics and Mathematical Geophysics Siberian Branch of the Russian Academy
of Sciences, 6 Prospect Akademika Lavrentieva Street, Novosibirsk, 630090,
Russia;  Novosibirsk State University, 2 Pirogova Street, Novosibirsk, 630090,
Russia;  Institute for Disease Modeling, Bill & Melinda Gates Foundation,
Seattle, USA.  University of Leeds, LS2 9JT, UK.

###### Abstract

This paper uses Covasim, an agent-based model (ABM) of COVID-19, to evaluate
and scenarios of epidemic spread in New York State (USA), the UK, and the
Novosibirsk region (Russia). Epidemiological parameters such as contagiousness
(virus transmission rate), initial number of infected people, and probability
of being tested depend on the region’s demographic and geographical features,
the containment measures introduced; they are calibrated to data about
COVID-19 spread in the region of interest. At the first stage of our study,
epidemiological data (numbers of people tested, diagnoses, critical cases,
hospitalizations, and deaths) for each of the mentioned regions were analyzed.
The data were characterized in terms of seasonality, stationarity, and
dependency spaces, and were extrapolated using machine learning techniques to
specify unknown epidemiological parameters of the model. At the second stage,
the Optuna optimizer based on the tree Parzen estimation method for objective
function minimization was applied to determine the model’s unknown parameters.
The model was validated with the historical data of 2020. The modeled results
of COVID-19 spread in New York State, the UK and the Novosibirsk region have
demonstrated that if the level of testing and containment measures is
preserved, the number of positive cases in New York State and the Novosibirsk
region will remain the same during March of 2021, while in the UK it will
reduce. Due to the features of the data for the Novosibirsk region (two
datasets are characterized as stationary series with probability of 1), the
forecast precision is relatively high for the number of hospitalizations, but
is lower for new cases of COVID-19.

###### keywords:

epidemiology, agent-based modeling, COVID-19, interventions analysis,
coronavirus data analysis, New York State, United Kingdom, Novosibirsk region,
mathematical modeling, forecasting scenarios, reproduction number,
optimization, OPTUNA software.

††journal: Journal of LaTeX Templates

1abstract

##  1 Introduction

In December 2019, the COVID-19 pandemic originated in the Wuhan province of
China. Since that time, more than 130 million people in 192 countries have
been infected with the disease, and more than 2.87 million people have died
after getting infected. During the year 2020, mankind mobilized its resources
to fight the pandemic. One of the useful tools in this struggle has been
mathematical modeling that uses known historical data to study different
scenarios of disease spread [1, 2].

The models including those for studying coronavirus infections can be divided
into two groups: compartmental and agent-based models. In compartmental
models, a population is divided into groups sharing similar features and
interacting with one another following the mass action law. Agent-based models
(ABMs), on the other hand, give each agent (people, social institutes, the
state, etc.) a set of features and determines the way the agents interact from
random graphs following disease spread principles. In other words, an agent’s
behavior is determined individually, and their joint behavior is described as
the interaction of multiple agents (bottom-top approach). Unlike compartmental
models, ABMs are capable of providing a detailed description of an
epidemiological situation, especially in a case of inaccurate and insufficient
data. ABMs account for the stochastic nature of epidemic spread, makes it
possible to estimate the likelihood of different epidemic scenarios, and
allows one to evaluate the risks of unfavorable events occurring due to policy
changes. The resulted data enable one to make conclusions about the duration,
severity and scale of an epidemic, evaluate the efficacy of the preventive and
quarantine measures, and assess its economic consequences.

During 2020, many models were developed to predict COVID-19 spread. In [3],
both compartmental and agent-based models are presented to study mass face-
mask wearing and predict its effect on COVID-19 spread. A graph-based ABM is
suggested in [4]. The paper considers a small population of 1000 agents and
graphs of different kinds such as fully-connected, Barabasi-Albert, Watts-
Strogatz ones, etc. Another approach to agents interactions is demonstrated in
[5]. In this paper, agents are initialized in a 2D space (so-called mesh) to
consider a distance between different agents, so such a factor as social
distancing can be explicitly accounted for. However, this approach leads to
excessive computation complexity, and for that reason, these models consider a
small number of agents (up to 500).

Another important direction of modeling has been a comparison of containment
measures and projecting the future for different scenarios. The paper
published by Chang et al. [6] combines an epidemiological SEIR model and the
hourly GPS data from the mobile phones of 98 million people in 10 US cities.
The model predicted that closing of the most crowded public places such as
restaurants and religious establishments would be a sufficient measure to
contain the pandemic unlike unilateral measures to limit people’s mobility. In
[7], the authors offer an ABM using SEIR agents to model COVID-19 spread
dynamics, where the agents imitate people, businesses and the government.
Using the model, they have analyzed seven social-distancing scenarios having
different epidemiological and economic effects. The paper has demonstrated
that the so-called vertical isolation has no positive effect. In [8]
macroeconomic epidemiological ABM been presented to study the economic effect
COVID-19 would have in different scenarios of pandemic containment such as
closing of educational and entertainment facilities. The model was calibrated
using the statistical data on country’s and business demography, households,
employment, profits and wages in Germany.

It is noteworthy that data collecting and processing is a very important step
in building an effective COVID-19 spread model. However, in the studies
mentioned above data pre-processing for the modeled regions was not performed.
Most of them concentrated on building the models and algorithms, whose
parameters were considered known either from literature or from experts’
estimations, so the issues of identifiability for unknown parameters have
remained unresolved, as has the need of devising a regularization algorithm
for solving the problem of epidemiological forecasting.

In this study, our focus is on the analysis of data, parameter identification,
and regularization algorithms. However, it is a known fact that epidemics
develop differently in different locations. To address this issue, in our
study, the epidemiological situations in New York State (USA), the UK, and the
Novosibirsk region (Russia) were compared and analyzed.

###  1.1 New York State

Nearly 2 million people were confirmed infected and more than 50000 people
died by April 6, 2021 in NY State [9]. In NY State, the pandemic spread
rapidly, reaching its peak in March-April of 2020 (see Section 2.1). The
healthcare system was overloaded in the very first months of the outbreak.
Thanks to the containment measures introduced, the number of infected people
reduced to a characteristic plateau that was followed by a second infection
wave several months later.

###  1.2 United Kingdom

Another pandemic scenario was observed in the UK. After the first infection
wave had been successfully suppressed by June, 2020, the second wave hit the
country hard and, due to the B.1.1.7 SARS-CoV-2 variant, the number of
infected people increased 10 times compared to the first wave (see Section
2.3). At the time of writing, this has reduced to 7000 cases a day, and the
pandemic is on wane.

###  1.3 Novosibirsk region

The region’s peculiar feature was the absence of the characteristic peaks that
rather looked like plateaus (see Section 2.3), so the lockdown introduced in
April 2020 had outlived its usefulness by the summer of 2020.

This study is organized as follows. First, open-source data of COVID-19 spread
in New York State, the UK and the Novosibirsk region were processed and
analyzed using the statistics and machine learning methods to find
interdependencies, study seasonality and predict possible future dynamics (see
Section 2). Second, we confirmed that the selected ABM met the identifiability
condition such as being sensitive to data errors and capable of unambiguous
determination of the unknown parameters of COVID-19 spread from additional
measurements (Section 3.2). The obtained space of identifiable parameters was
specified using the multilevel global-optimization method (Section 4).
Finally, scenarios of how the COVID-19 pandemic that could develop in New York
State, the UK and the Novosibirsk region were assessed concerning available
data and certain containment measures (Section 5).

##  2 Data analysis

Data analysis and data processing are important parts of forecasting modeling.
Before data processing begins, one has to understand the character of
available data and determine their features. Collecting such daily indicators
as the number of tests, diagnosed cases, ventilated COVID patients, etc. (see
Table 1) helps to overview the general picture for a considered region, while
anomalous time intervals may call for more scrupulous analysis.

Table 1: Indicators used in data analyses and their description. Indicator |  Description
---|---
New Tests |  Number of performed tests
New Diagnoses |  Number of diagnosed cases
New Deaths |  Number of deaths related to a positive COVID-19 diagnosis
Num Critical |  Number of ICU patients with suspected and positive COVID-19 diagnosis

###  2.1 New York State (USA)

To forecast the way the pandemic would develop, the data from the COVID
Tracking Project’s website [11] were used. The site contains information for
each state and for the country as a whole. The feature of the data in question
is the method they calculate New Diagnoses: positive cases (confirmed plus
probable) summing the total number of confirmed cases and the probable cases
of COVID-19 reported by the state or territory, ideally per the ”August 5,
2020 CSTE case definition”. Some states are following the older ”April 5th,
2020 CSTE case definition” or using their own custom definitions. The latter
method of data collection is more suitable for ABMs since it takes into
account the percentage of infected people that may have been simply neglected.

We can consider the COVID-19 spread in New York State in more detail, using
the epidemic time-series data. Any time series can be decomposed into the
following three elements:

$X(n)=T(n)+S(n)+N(n).$ (1)

Here, $X(n)$ is a time-series value, $T(n)$ is the value of the underlying
time-series component; $S(n)$ is the value of a seasonality component, $N(n)$
is the value of a noise component for the $n$-th day. When analyzing a time
series, we found it was most useful to analyze its trend, since this
determines an indicator’s behavior in time. For expansion by the expression 1,
the Hodrick-Prescott filter was applied [12] (see result in Fig. 1).

Figure 1: Graphs demonstrating the expansion results of the New Tests: $T(n)$
(top) and $S(n)+N(n)$ (bottom).

The non-smoothed graphs in Fig. 2 demonstrate widely dispersed points of
statistics, so these data were smoothed before using them in the model because
only the main trends of the curves were necessary for reaching an appropriate
result. The graphs demonstrate certain periodicity that is known to be time
series seasonality, which is clearly traced during summer in New Diagnoses.

Figure 2: Graphs of COVID-19 spread in New York State (top line): the number
of New Tests (left), New Diagnoses (middle) and New Deaths (right) and their
smoothed curves (bottom line).

Now, let us consider the average fraction of tests for each day of the week
that is calculated as:

$w_{i}=\sum\limits_{j=0}^{N}\frac{X(i+m\cdot j)}{S(i+m\cdot j)}.$

Here, $i$ is a week-day number, $N$ is the number of full weeks within a
considered time series, $S(i+m\cdot j)$ is the cumulative sum of the New Tests
performed within a week corresponding to index $j$, $m=7$. The results can be
seen in Table 2. It is apparent that the number of tests tends to its minimum
for Mondays and Sundays and reaches its maximum on Fridays and Saturdays,
which exemplifies the seasonality (i.e., periodicity) of the considered time
series.

Table 2: Average fraction from the number of tests and its daily distribution for a whole statistics-gathering period (NY State). Days of the week |  Average fraction
---|---
Monday |  0.113654
Tuesday |  0.127697
Wednesday |  0.134512
Thursday |  0.155779
Friday |  0.162017
Saturday |  0.163702
Sunday |  0.142640

COVID-19 pandemic data analysis has shown that there is a dependence between
New Tests and New Diagnoses datasets for different regions. In such cases, the
behavior of the second indicator is partly determined by that of the first. As
a matter of fact, the last day of the week does not mean an abrupt increase in
the number of the infected only because it is Friday. However, more people get
tested on Friday, so the number of positive tests may increase. For that
reason, analyzing the links between the number of infected people and that of
new tests becomes crucial. For a proper understanding of the situation, it is
not New Diagnoses, but their fraction from the number of New Tests that we
want to know. After all, if one tested every person in a region, they would
immediately indicate every one infected and their indicator would increase
abruptly while their percentage would remain the same. Based on this fact, one
can derive a time series that describes the percentage of newly tested people
with either positive or potentially positive COVID-19 test.

Figure 3: Fraction of New Diagnoses from New Tests in NY State.

Figure 3 shows that the ratio reached its peak in April 2020 and later started
to reduce due to an increased number of tests. In autumn, the second pandemic
wave began, so the number of daily confirmed cases increased (see Fig. 2), but
the ratio between New Diagnoses and New Tests remained at a low level. This
implies that testing behavior (i.e., the probability of testing with or
without COVID-19 symptoms) did not change. The relatively few number of deaths
over this period was due to the epidemic primarily spreading among younger age
groups (see Fig. 2).

To confirm change in testing during autumn and winter, the MACD indicator [13]
was used, whose histogram tracks a function’s rise and fall (in our case it
was the rise of New Tests). The graph demonstrates that the MACD indicator
increased significantly over this period (see Fig. 4).

Figure 4: MACD-indicator of the New Tests in NY State.

####  2.1.1 Seasonality

The above-mentioned term “seasonality” refers to the periodic fluctuations
observed in time series. In other words, if one takes a time-series space and
overlaps it against a neighboring space of the same size, their profiles are
going to coincide (or differ by a constant). The peak absolute values will
correlate with the same time points calculated from the beginning of the
space. As for our COVID-19 data, their seasonality can be traced in the New
Tests indicator, which has its logic since collecting this statistics involves
as a human factor as the features of the healthcare system. The seasonality of
a time series makes it possible to forecast the series behavior relative to
some average value.

Seasonality is commonly determined with an autocorrelation function (ACF). For
a discrete process $X(1),X(2),\ldots,X(n)$ its formula is written as:

$R(n)=\frac{1}{(M-n)\cdot\sigma^{2}}\sum\limits_{t=1}^{M-n}(X(t)-\mu)(X(t+n)-\mu).$

Here $\sigma$ is the standard deviation of a discrete process $X$, $\mu$ is
its average value, $M$ and $n$ are positive integers. Besides, to analyze
time-series seasonality, a partial correlation (PACF) was applied that removed
the linear dependence between shifted time series.

Figure 5: Graphical results of applying ACF (left) and PACF (right) to New
Tests (top) and New Diagnoses (bottom) in NY State.

The results confirmed New Tests really had weekly seasonality (7, 14 and 21st
days), while the New Diagnoses did not possess this property (see Fig. 5).

####  2.1.2 New Tests/New Diagnoses interrelation

We have also considered a percentage change for any current moment in relation
to the same moment a week before:

$pc(n)={X_{smoothed}(n)\over X_{smoothed}(n-7)}-1.$

This statistic demonstrates by how much an average indicator has changed in
fractions compared to a previous week. If compared for weekdays, the results
become smoother and easier to interpret.

Figure 6: New Diagnoses percent change (red line) and New Tests percent change
(blue line) time-series graphs in NY State.

The two time series in Fig. 6 have spaces where their trajectories almost
match. In other words, the percent change of one indicator differs from that
of the other one by a constant. In terms of tested/infected ratio, such spaces
confirm that within this period, the number of infected people grew owing to
the increased number of tests and not to a worsening pandemic situation in the
region.

For more specific analysis, a window of 28 days was correlated with the
previous 28 days for every day in the window (see Fig. 7).

Figure 7: New Diagnoses percent change and New Tests perecent change
correlations within a 28-day backward window with $\mid 0.8\mid$ and $\mid
0.9\mid$ correlation’s threshold lines for NY state.

In statistics, a significant linear dependence occurs when the absolute
correlation value exceeds 0.8. If it exceeds 0.9, it means there is a strong
dependence between two indicators. In NY State, with some deviations, New
Diagnoses strongly depended on New Tests. The deviations mean the pandemic
develops following its own scenario. For instance, in the end of May and the
beginning of June, the correlation coefficient exceeded 0.8 within a 28-day
backward window when the first pandemic wave in the region was considered
defeated and reached the so-called plateau that lasted till October when the
second correlation took place and the second pandemic wave began, which means
the second wave might have been triggered by the abrupt increase of New Tests.
However, the correlation coefficient got back to almost zero by December,
indicating that the rise and fall of New Diagnoses did not correlate with New
Tests, so the pandemic in the region at that time spread or reduced (depending
on graph direction). In February and March, the correlation became very close
to one, meaning the number of infected people depended only on a New Tests.

###  2.2 Novosibirsk region

Figure 8: COVID-19 spread graphs for the Novosibirsk region (top line): New
Tests (left), New Diagnoses (middle) and New Deaths (right) and their smoothed
curves (bottom line).

The data (see Fig. 8) for our research were taken from the RBC website that
published the daily statistics of COVID-19 spread in the Novosibirsk region
(the used data are given in Table 1). In March and May of 2020, these data had
a significant number of gaps. For instance, the statistic on the number of
performed tests was available only starting April 25, 2020, as the other
statistics were accumulated starting March 17, 2020. For that reason, the New
Tests was extrapolated in relation to the abrupt increase in the number of New
Diagnoses starting April 6, 2020. For filling the gaps, the following backward
extrapolation algorithm was iterated:

$X(n)=\frac{1}{m}\sum\limits_{j=1}^{m}X(i+j)\cdot C^{n}+E(n),\quad
i=L,\ldots,0.$

Here, $X(n)$ is the value of a time series of the $n$-th day; $m=7$, $C=1.03$,
$L$ is the number of gaps to be extrapolated; $E(n)$ is a random value of
normal distribution, i.e. $E(n)\in\mathcal{N}(0,{X(n)\over 3})$. The
extrapolation result for New Tests can be seen in Fig. 9.

Figure 9: Backward extrapolation results for the New Tests (right) compared
with the actual values of New Diagnoses time series in the Novosibirsk region
(left).

####  2.2.1 Seasonality

Here, the approach similar to that applied for the NY State data was
implemented to calculate the ACF and PACF for several time series (see Fig.
10), together with the fraction of tests performed for every day of the week
(see Table 3).

Figure 10: Graphical results of applying ACF (left) PACF (right) to the New Tests (top) and New Diagnoses (bottom) in the Novosibirsk region. Table 3: Average fraction from the number of tests and its daily distribution for a whole statistics-gathering period in the Novosibirsk region. Week days |  Average fraction
---|---
Monday |  0.064055
Tuesday |  0.176845
Wednesday |  0.148333
Thursday |  0.161924
Friday |  0.162173
Saturday |  0.183951
Sunday |  0.075146

As can be seen from Fig. 10, both functions confirmed the seasonality of New
Tests with the duration of the season to be one week. As it has been discussed
earlier, for this indicator this is absolutely normal. Table 3 demonstrates
that the lowest fraction of tests was registered on Sunday and Monday, and the
highest one – on Tuesday and Saturday.

####  2.2.2 Stationarity

A time series is regarded as stationary (or weakly stationary) if its expected
value and variance are time-independent, and its autocorrelation function
depends only on a difference of neighboring values. Let us consider two time
series for the New Diagnoses: from May 29 to July 1 and from October 24 to
December 7 (see Fig. 8). Having applied the augmented Dickey-Fuller test,
whose zero hypothesis is non-stationary time series (Unit root test), we find
that this hypothesis can be rejected for both time series intervals, since
their p-values are 0.0 and 0.001, respectively. If time series are stationary,
the considered data sets can probably be interpreted as a normal noise. Table
4 demonstrates the results of two statistical tests to verify the normality of
the considered time series.

Table 4: Results of the two statistical tests, whose zero hypothesis is the normal distribution of two stationary parts of New Diagnoses. Time interval |  Kolmogorov-Smirnov test ($p$-value) |  D’Agostino-Pearson test ($p$-value)
---|---|---
29.05.2020 - 1.07.2020 |  0.3842 |  0.1827
24.10.2020 - 7.12.2020 |  0.573 |  0.8193

According to the test results, the null hypothesis of the normality of the two
areas is not rejected. However, due to the small number of observations, it is
equally impossible to assert that the values in the selected areas are
normally distributed.

####  2.2.3 New Tests/New Diagnoses interrelation

Calculated New Diagnoses and New Tests percent change correlations (see Fig.
11) showed there was not even a single window of 28 days with a significant
New Tests/New Diagnoses linear dependence. In the middle of June, a negative
correlation peak was observed confirming the “more tests less positive cases”
tendency lasting from the end of June to the beginning of July. Another
interesting fact was the New Diagnoses, percentage that plummeted in August
while New Tests remained at the same level.

Figure 11: New Diagnoses to New Tests ratio (left); New Diagnoses (red)
percent change and New Tests (blue) percent change (center) and their
correlations within a 28-day backward window (right) with $\mid 0.8\mid$ and
$\mid 0.9\mid$ correlation’s threshold lines for the Novosibirsk region.

###  2.3 United Kingdom

To analyze COVID-19 spread in the UK, the data accumulated on an official
government web portal were used [15].

Figure 12:  COVID-19 spread graphs in the UK (top line): New Tests (left), New
Diagnoses (middle), New Deaths (right) and their smoothed curves (bottom
line).

Fig. 12 shows two infection waves: spring (from 2020.03.09 to 2020.06.01) and
fall (from 2020.09.01 to 2021.03.01) ones. The standard deviation of the time
series for the second wave was much bigger than for the first one. For that
reason, the standard deviation values were considered for three independent
time series of New Diagnoses and New Tests (see Table 5). The data for NY
State (NY) and the Novosibirsk region (NSK) are put in the table for
comparison.

Table 5: Standard deviation values for particular New Diagnoses and New Tests time series in the UK, NY State and the Novosibirsk region. Indicator |  First |  Second |  Third |  Fourth
---|---|---|---|---
(UK) New Tests std |  10626.7 |  16613.8 |  77382.5 |  …
(UK) New Diagnoses std |  455.8 |  201.1 |  5727.2 |  …
(NY) New Tests std |  5636.6 |  9543.9 |  31127 |  …
(NY) New Diagnoses std |  1137 |  114.2 |  1144 |  …
(NSK) New Tests std |  765.2 |  1271.2 |  1596 |  2492
(NSK) New Diagnoses std |  7 |  4.6 |  7 |  6

Unfortunately, the data did not allow us to conclude what was the exact reason
for the second wave’s higher variance. It could have been the week variance of
New Tests or something else since, in the UK, they started to register the
number of tests only on 2020.04.21. All we know is the standard deviation of
New Diagnoses increased with time as well as the variance of New Tests.
However, both in NY State and the Novosibirsk region, despite the growing
variance of New Tests, the standard deviation value of New Diagnoses remained
at the same level.

####  2.3.1 Seasonality

The ACF and PACF applied to the New Diagnoses and New Tests time series and
the fractions of tests calculated for every day of the week demonstrated the
weekly seasonality of New Tests confirmed by ACF (see Fig. 13 and Table 6).
However, no significant PACF delays were observed. Comparing the results
obtained with those from NY State and the Novosibirsk region demonstrated the
absence of linear dependence did not affect the time series’s seasonality.
However, a different situation was observed in the UK.

Figure 13: Graphical results of applying ACF (left) and PACF (right) to New Tests (top) and New Diagnoses (bottom) for the UK. Table 6: Average fraction from the number of tests and its daily distribution for a whole statistics-gathering period in the UK. Days of the week |  Average fraction
---|---
Monday |  0.125222
Tuesday |  0.135905
Wednesday |  0.159167
Thursday |  0.166208
Friday |  0.157837
Saturday |  0.131992
Sunday |  0.123668

In the next section, an approach to forecasting the New Tests time series will
be demonstrated. The approach relies upon several techniques, including
SARIMA, an algorithm requiring a significant correlation for the seasonal
delays (7, 14, 21, etc.) of a single parameter.

###  2.4 Forecasting New Tests

To draw forecast curves while modeling, one has to predict a number of certain
statistical data sets that are used as input parameters. In our model, such a
data set was the New Tests performed in the region since this indicator did
not depend on the others and had the highest value of seasonality that
determined the seasonality of the other indicators.

The first model considered for forecast purposes was SARIMA [16], a
modification of ARIMA [17] (AutoRegressive Integrated Moving Average), which
is able to support time-series data sets with a seasonal component. ARIMA is
an extension of the ARMA models for non-stationary time series that can be
converted into stationary ones through differencing of a certain order from an
initial time series (so-called integrated or differential-stationary time
series).

For a non-stationary time series $X(n)$, the ARIMA$(p,d,q)$ model is written
as:

$\triangle^{d}X(n)=c+\sum_{i=1}^{p}a_{i}\triangle^{d}X(n-i)+\sum_{j=1}^{p}b_{j}\epsilon(n-j)+\epsilon(n).$

Here, $\epsilon(n)$ is a stationary time series of white noise;
$c,a_{i},b_{j}$ – model parameters; $\triangle^{d}$ – difference operator of a
time series of order $d$ (sequential taking of $d$ times of first-order
differences: firstly, from a time series, then from obtained first-order
differences, then from second-order differences and so on).

Parameters $p$, $d$, $q$ and $P$, $D$, $Q$ are determined through cross-
validating minimization of the AIC-metrics. Here, $p$ is the number of the
last non-seasonal delay with significant PACF; $P$ is the number of the last
seasonal delay with significant PACF; $q$ is the number of the last non-
seasonal delay with significant ACF; $Q$ is the number of the last seasonal
delay with significant ACF; $d$ is a differentiation order with 1-day delay,
$D$ is a seasonal differentiation order with 7-day delay. The ARIMA model is
based on the following idea: one builds a model from the differences of a
value for $d$ sequential periods to obtain a stationary process. In a general
case, the difference order is limited to $d=2$, since taking the second-order
differences allows converting almost any non-stationary data series into
stationary ones.

As for the time-series forecast, the following algorithm was applied:

  1. 1.

Applying the Box-Cox transformation [18] to reduce the variance.

  2. 2.

Calculating the first-order seasonal difference (7-day shift).

  3. 3.

Calculating the difference (1-day shift) of the series obtained at step 2.

  4. 4.

Applying the Dickey-Fuller criterion to verify the stationarity of the series
obtained at step 3.

  5. 5.

Passing the hyperparameters corresponding to the actions performed and
selecting other ones using the minimized Akaike information criterion. As
data, the series from step 1 is passed.

  6. 6.

The model with adjusted hyperparameters was used for forecasting. The obtained
results were treated using the inverted Box-Cox transform.

The second model considered for the forecast was Holt-Winters [19], which
provides triple exponential smoothing of a time series. The idea behind the
model is to break a time series into such components as a level, trend and
seasonality. The last component explains repeating fluctuations around the
level and the trend and is characterized by its duration, i.e. a period after
which a fluctuation repeats itself. In this approach, each observation has its
own component, e.g. if a season’s duration is 7 days, there are 7 seasonal
components, one for every day of the week.

To build the model, parameters $\alpha$, $\beta$ and $\gamma$ had to be
determined. $\alpha$ is responsible for smoothing the series around the trend,
$\beta$ is for smoothing the trend itself, $\gamma$ is the seasonality
component. Minimizing mean absolute error (MAE) through time-series cross-
validation, we selected the model’s optimal parameters (iterative three-
parameter random search).

To extrapolate the New Tests time series from the Novosibirsk region, in
addition to the models described above, a linear regression model was used for
comparative analysis. In this respect, for different time intervals, SARIMA’s
forecast results were on average not better but more stable than those
produced by the other models (see Fig. 14).

Figure 14: Monthly forecast results for the New Tests time series (dark blue line) for 2021.01.18 obtained using the Holt-Winters (black dashed line), linear regression (blue dashed line) and SARIMA (green dashed line) models for the data from the Novosibirsk region. Table 7: Mean Absolute Error (MAE) between the New Tests time series and a corresponding prediction (from May 5, 2020 to January 18, 2020). Forecast model |  MAE
---|---
Holt-Winters model |  1646.49
SARIMA model |  1454.73
Linear Regression model |  1586.34

As a result, the SARIMA forecast gave the lowest MAE (see Table 7). This
algorithm also produced the most accurate seasonality and time-series
autoregression. For that reason, a decision was made to use SARIMA for
forecasting the future values of New Tests.

##  3 Agent-based mathematical model

This section presents the ABM devised to describe COVID-19 spread and
formulates a problem to identify (calibrate) the model’s unknown parameters as
objective functional minimization. It also presents a scheme for automatic
calibration of the parameters for time intervals exemplified by the number of
infected people in the Novosibirsk region.

###  3.1 ABM formulation

Within the framework of this research, stochastic ABMs for New York State, the
United Kingdom and the Novosibirsk region were devised. They were built using
the Covasim [20]. This package had been utilized to predict the number of
infected, dead and hospitalized people in the State of Oregon and become one
of the tools to make decisions about whether to relax or escalate COVID-19
containment measures [21, 22]. This library is written in Python to study non-
trivial COVID-19 dynamics. Its general algorithm is as follows: after all
necessary parameters and statistical data are uploaded, the package creates an
artificial population with account for age distribution. The model’s agent is
a person in a particular region. Then, the agents are united into contact
networks and the integration loop begins. At every time step (1 day), an
agent’s status is updated in relation to its contact network and the
containment measures relevant for this interval (self-isolation; closed access
to public places; wearing face masks, etc.). The agents can interact with one
another in particular networks. Depending on the network’s structure, both
full and random connectivity graphs are built for an agent to figure out how
the infection spreads. The average number of daily contacts is different for
every agent and every network. At any time moment, the agents distributed by
their age (bins of 0-9, 10-19, …, 90+) are found in their given state (see
fig. 16). Our ABM also accounted for tested agents according to the data on
COVID-19 tests performed in the Novosibirsk region. More details about the
structure of Covasim-based ABMs, their parameters and realization methods can
be found in [23].

Every agent has their set of properties and characteristics that can be
divided into 2 groups: constant (belong to each particular agent and do not
change while modeling) and time-dependent.

####  3.1.1 Time–independent agent characteristics

  * 1.

Age($t^{*}$). All the agents are subdivided into age groups of 10 years (0-9
years, 10-19, …, 90$+$. The age distribution depends on the demographic
situation in a studied region.

  * 2.

Social status (determined by an agent’s age $t^{*}$). Depending on their age,
agents contact one another in contact networks. All agents have contacts in
households and public places. Agents of 6-21 years old can also have contacts
in educational institutions with agents of their age. Agents of 22-65 years
old contact at work (see Fig. 15). Depending on a contact’s structure, the
transmission parameter $\beta$ is multiplied by corresponding constant
$w_{\beta}$ ($w_{\beta}=3$ for households, 0.6 – for educational institutions,
0.3 – for public places), i.e. the likelihood of virus transfer is different
for every network.

Figure 15: Agent contacts layers and their interactions in the COVID-19 spread
ABM.

  * 3.

Likelihood of disease progression (determined by an agent’s age $t^{*}$).
These parameters determine disease progression (see Fig. 16). Their
description is given in Table 8.

Table 8: Parameters of disease progression probability. Parameter |  Description
---|---
$p_{sym}(t^{*})$ |  Probability of developing symptoms
$p_{sev}(t^{*})$ |  Probability of developing severe symptoms (requires hospitalization)
$p_{crit}(t^{*})$ |  Probability of critical condition (requires ICU)
$p_{death}(t^{*})$ |  Probability of death

####  3.1.2 Time-dependent agent characteristics

  * 1.

Agent’s epidemiological status. Each agent may have one of the 10 stages of
the disease (Fig. 16) $\vec{X}=(S,E,I,A,Y,M,H,C,R,D)$.

  * 2.

Agent’s chance to be tested for COVID-19 ($\tilde{p}(X(t))$) that is
determined by the agent’s epidemiological status. The agents are tested daily,
the number of tests corresponds to the statistical data obtained in a region.
At every modeling step, the tests are distributed across the population, the
agents whose status is marked with an orange frame in a Fig. 16 can be given a
positive result. The agents whose test is recognized as positive are marked as
“confirmed” and included in New Diagnoses. The model assumes that the
likelihood for an agent to be tested as a symptomatic carrier is higher and
this chance ratio is controlled by parameter $\tilde{p}(X(t))$ at is restored
from solving the inverse problem (see the next section).

Figure 16: Agent state transition diagram in Covasim that is based on a SEIR-
type compartment model. The orange frame marks those agent states that can
give positive COVID-19 tests.

###  3.2 Parameter identification problem

The ABM developed was also characterized by unknown parameters vector
$\vec{q}$. To specify the model’s parameters a variational inverse problem
formulation was performed to minimize the misfit function:

$J(\vec{q})=\sum_{s}\sum_{t_{i}=1}^{T}{|X^{s}_{d}(t_{i})-X^{s}_{m}(t_{i},\vec{q})|\over
M_{s}}.$ (2)

Here, $s$ is the statistics used for comparison, e.g. New Diagnoses, New
Deaths, etc.; $X^{s}_{d}(t_{i})$ and $X^{s}_{m}(t_{i},\vec{q})$ are daily
statistical and model data for statistics $s$; $T$ – the number of modeled
days, $M_{s}=\max\limits_{{t_{i}}}\\{X^{s}_{d}(t_{i})\\}$ is a normalising
term.

###  3.3 Automatic parameter calibration

The model assumed that parameter $\beta$ was a piece-wise constant, so the
days of change and new values were determined by parameters $\vec{\beta}_{d}$
and $\vec{\beta}_{c}$. The longer was a considered time interval, the more
unknown parameters it included. Since every launch of the model’s calibration
algorithm was rather time-consuming, the time interval in question was divided
into periods of 1 month. For example, for NY State the first period was
2020.03-02 - 2020.04.01, the second – 2020.04.02 - 2020.05.0; for the UK –
2020.02.07 - 2020.03.08 and 2020.03.09 - 2020.04.07; for the Novosibirsk
region - 2020.03.12 - 2020.04.11 and 2020.04.12 - 2020.05.11 etc.

Each period was sequentially calibrated, so the parameters restored at a
previous step were used in the following iteration of the optimization
algorithm. Thus, for the initial period considered the unknown parameter
vector was

$\vec{q}_{1}=(E(0),\beta,\beta_{d}(1),\beta_{c}(1),\tilde{p}(X)),$

where $E(0)$ is the initial number of infected agents, $\beta$ is a
contagiousness parameter value, $\beta_{d}(1)$ is the day parameter $\beta$
changes, $\beta_{c}(1)$ is the value by which parameter $\beta$ changes on day
$\beta_{d}$, $\tilde{p}(X)$ is a test level parameter in relation to
statistical data. For all the following periods (second, third, etc.):

$\vec{q}_{i}=(\beta_{d}(i),\beta_{c}(i)).$

Figure 17 displays a phased restoration of the vector of unknown parameters
$\vec{q}$ for the Novosibirsk region exemplified by 4 out of 13 interim
calibration periods.

Figure 17: Phased restoration of parameters $\vec{q}$ using the New Diagnoses
dataset from the Novosibirsk region. The vertical dashed lines mark the days
when the parameter $\beta$ changed that were restored from solving an
optimization problem.

##  4 Methods and approaches

Identifying (calibrating) the parameters of an ABM, so that its outputs match
observed data, is quite a complex and computation-intensive task due to the
large number of the parameters involved. There are different approaches to the
problem. In most cases, the parameters are selected manually or one uses
averaged experimental results neglecting the specific features of a studied
region are used. Analysis of the papers describing epidemiological ABMs has
demonstrated that no algorithm could be considered superior for identification
of model parameters [24]. According to the paper “…it appears that calibrating
individual-based models in epidemiological studies of HIV, malaria and TB
transmission dynamics remains more of an art than a science."

In our model, the vector of unknown parameters $\vec{q}$ was calibrated using
the Optuna hyperparameter optimization software [25] to be one of the latest
optimizers designed to adjust hyperparameters in machine learning algorithms
and neural networks. The optimizer is based on the tree-structured Parzen
estimator (TPE) that in many ways is similar to the Bayesian optimizer [26].
However, unlike the Bayesian optimizer that calculates
$p(J(\vec{q})|\vec{q})$, TPE calculates $p(\vec{q}|J(\vec{q}))$ and
$p(J(\vec{q}))$ to determine the parameters domain to minimize functional $J$
by performing Parzen window density estimation, to generate two separate
distributions specifying the high - and low-quality regions of the input-space
respectively.

The TPE works by ranking the accumulated observations
$\mathcal{D}_{K}=\\{q_{k},J(q_{k}):k=1,...,K\\}$ based on the objective
function values. The algorithm selects where $J$ should be evaluated next by
splitting $\mathcal{D}_{K}$ into two sets, $\mathcal{D}^{l}_{K_{l}}$ and
$\mathcal{D}^{g}_{K_{g}}$ such that $\mathcal{D}^{l}_{K}$ contains the
$\gamma$-percentile of the highest quality (i.e. lowest function value) points
of $\mathcal{D}_{K}$, whereas $\mathcal{D}^{g}_{K}$ contains the remaining
points. With these definitions, the prior for $J$ is given by
$P(J<J^{\gamma})=\gamma$. Now, applying non-parametric adaptive parzen
windows, two probability distributions $l(\vec{q})$ and $g(\vec{q})$ are
estimated from $\mathcal{D}^{l}_{K_{l}}$ and $\mathcal{D}^{g}_{K_{g}}$,
respectively. Here, $l(\vec{q})$ is interpreted as representing the
probability of a region in the input space yielding a high-quality observation
while similarly, $g(\vec{q})$ represents low-quality regions. The likelihood
of $\vec{q}_{K+1}$ belonging to distribution $l$ and $g$ may now be expressed
by:

$p(\vec{q}_{K+1}|J(\vec{q}))=\begin{cases}l(\vec{q}),&J(\vec{q}_{K+1})<J^{\gamma},\\\
g(\vec{q}),&J(\vec{q}_{K+1})\geq J^{\gamma},\end{cases}$ (3)

where $J^{\gamma}$ is the lowest function value found in
$\mathcal{D}^{g}_{K_{g}}$.

The main optimization criterion in the TPE is an expected improvement (EI).
According to the definition [27], the EI is a value the function reduces by at
a given moment and can be written as:

$EI(\vec{q})=\left(\gamma+\frac{g(\vec{q})}{l(\vec{q})}(1-\gamma)\right)^{-1}$
(4)

Finally, in order to select a new point to evaluate the objective function at,
samples are drawn from $l(\vec{q})$, where the sample generating the highest
expected improvement is used to evaluate the objective function ($n_{samp}$).
A full TPE optimization procedure is described in Algorithm 1.

Algorithm 1 Tree-Parzen estimator optimization

1:Parameter values for $\gamma$, $n_{samp}$ and $max\\_iter$

2:Inintialize: accumulate initial observations

3:$\mathcal{D}_{init}=\\{\vec{q}_{k},J(\vec{q}_{k}),k=1,\ldots,n_{init}\\}$

4:for $m$=0 to $max\\_iter$ do

5: Split $\mathcal{D}_{n_{init}+m}$ to generate
$\mathcal{D}^{g}_{m_{g}},\mathcal{D}^{l}_{m_{l}}$

6: Estimate $l(\vec{q})$ from $\mathcal{D}^{l}_{n_{init}+m_{l}}$

7: Estimate $g(\vec{q})$ from $\mathcal{D}^{g}_{n_{init}+m_{g}}$

8: Draw $\vec{q}^{s}={\vec{q}^{s}_{k}}:k=1,\ldots,n_{samp}$, where
$\vec{q}^{s}_{k}\sim l(\vec{q})$

9: $\vec{q}_{m+1}=\mbox{argmax}EI(\vec{q})$

10: Evaluate $J(\vec{q}_{m+1})$

11: Augment set of observations
$\mathcal{D}_{n_{init}+m}\leftarrow\mathcal{D}_{n_{init}+m+1}$

##  5 Modeling and forecasting

In this section, we consider mathematical models and scenarios of COVID-19
spread in NY State, (Section 5.2), the UK (Section 5.3) and the Novosibirsk
Region (Section 5.4)

###  5.1 Initial datasets

To build and analyze the ABMs in the 3 considered regions, the following data
were used:

  1. 1.

Information on population’s age distribution according to the local government
statistics;

  2. 2.

Information on the average family size according to the UN data [28];

  3. 3.

Statistical data on the people infected with COVID-19, who recovered and died
including the number of tests performed that were collected from:

     * (a)

The COVID Tracking Project (New York State):

https://covidtracking.com/data;

     * (b)

The official UK Government website for data and insights on Coronavirus
(United Kingdom):

https://coronavirus.data.gov.uk/;

     * (c)

The RBC website and are available for downloading (Novosibirsk region):

http://covid19-modeling.ru/data/novosibirsk-region-data.csv.

For every region, the modeling results for the New Diagnoses, New Deaths and
Num Critical (require ICU) datasets were analyzed. In Sections 5.2-5.4 one can
find the graphs of 30-day forecasts validated with historical data. The
forecasts have an 80% confidence interval to characterize 10% and 90%
quantiles. In addition, effective reproduction numbers were calculated for
every region. These numbers indicate how many persons an infectious agent
infects on average during the time it has been infectious and calculated as:

$R(t)=\frac{I_{N}(t)\cdot f}{I_{C}(t)},$ (5)

where $I_{N}(t)$ is the number of new infections on day $t$, $I_{C}(t)$ is the
number of actively infectious people on day $t$ and $f$ is the average
duration of infectiousness. If $R(t)<1$, the pandemic is considered to stop
spreading and keeps spreading otherwise.

###  5.2 COVID-19 spread simulation in NY State

The mismatch (also called error or loss) function (2) for New York State was
as follows:

$J(\vec{q})=\sum_{t_{i}=1}^{T}\left(\frac{|Y_{d}(t_{i})-Y_{m}(t_{i},\vec{q})|}{M_{diag}}+\frac{|D_{d}(t_{i})-D_{m}(t_{i},\vec{q})|}{M_{death}}\right).$

Here, $Y_{d}(t_{i}),Y_{m}(t_{i},\vec{q})$ are smoothed New Diagnoses with
confirmed COVID-19, $D_{d}(t_{i}),D_{m}(t_{i},\vec{q})$ are smoothed New
Deaths.

After identification of the parameters, the model was validated with
historical data. The New Diagnoses dataset for NY State from 2020.03.02 to
2020.08.31 was taken as model data. The produced forecast covered 50 days and
was compared against statistical data for the given test period (2020.09.01 –-
2020.10.20). The modeling results can be seen in Fig. 18.

Figure 18:  Forecast for New Diagnoses validated with historical data for NY
State. The black line displays the modeled period, and the red dots – the New
Diagnoses dataset.

Despite the size of the confidence interval that started to increase from
2020.09.01, the increment rate of New Diagnoses was close to the test data.

Figure 19 presents the results of a 30-day forecast with restored vector of
unknown parameters for NY State, where the dots mark the real data accumulated
from 2020.03.02 to 2021.03.01. The forecast for New Diagnoses assumed the rate
of daily tests remained unchanged.

Figure 19: Model calibration results for 10 simulations + a 30-day forecast
for New Diagnoses (two top graphs) and New Deaths (two bottom graphs) in NY
State. The shaded areas are 10% and 90% quantiles, the solid line – the median
value of modeling result, and dots – real data. The dashed vertical lines are
COVID-19 containment measures.

The effective reproduction number $R(t)$ calculated by (5) showed the pandemic
began its downturn in NY State in January 2021.

Figure 20: Effective reproduction number $R(t)$ for the New York State. The
red dashed line marks $R(t)=1$.

###  5.3 COVID-19 spread simulation in the UK

The mismatch function (2) for United Kingdom was as follows:

$J(\vec{q})=\sum_{t_{i}=1}^{T}\frac{|Y_{d}(t_{i})-Y_{m}(t_{i},\vec{q})|}{M_{diag}}.$

Here $Y_{d}(t_{i}),Y_{m}(t_{i},\vec{q})$ are smoothed New Diagnoses with
confirmed COVID-19.

Firtly, the parameters that minimizing mismatch function were identified. The
New Diagnoses dataset for the UK from 2020.02.07 to 2020.08.31 was taken as
model data. The produced forecast covered 50 days and was compared against
observed data for the given test period (2020.09.01 – 2020.10.20). The
modeling results can be seen in Fig. 21. It is noteworthy that the prediction
error after comparing with real data was only 1.2%.

Figure 21: Forecast for New Diagnoses validated with historical data for UK.
The black line displays the modeled period, and the red dots – the New
Diagnoses dataset. The red dashed line marks the test period.

Figure. 22 presents the results of a 30-day forecast for the UK with restored
vector of unknown parameters $\vec{q}$ as is described in Section 4. The
forecast for New Diagnoses assumed the rate of daily tests remained unchanged.

As for New Deaths, the modeling results differed from the real data by 50%,
which was related to the target functional $J(\vec{q})$ that contained
information only about New Diagnoses, while the error for this indicator was
comparable to that of New Deaths. To reduce the error, weighed information
about the New Deaths indicator is to be added in the functional. The problem
of determining the weight coefficients needs further consideration.

Figure 22: Model calibration results for 10 simulations + a 30-day forecast
for New Diagnoses and New Deaths in the UK. The shaded areas are 10% and 90%
quantiles, the solid line – the median value of modeling result, and dots –
real data. The dashed vertical lines are COVID-19 containment measures.

The effective reproduction number $R(t)$ calculated by (5) for a time period
from 2020.02.07 to 2021.03.01 in the UK can be seen in Fig. 23. The graph
demonstrates the UK had a second infection wave in fall 2020 and in January
2021, after which the pandemic began its downturn.

Figure 23: Effective reproduction number $R(t)$ for the UK. The red dashed
line marks $R(t)=1$.

###  5.4 COVID-19 spread simulation in the Novosibirsk region

The mismatch function (2) for Novosibirsk region was as follows:

$J(\vec{q})=\sum_{t_{i}=1}^{T}\frac{|Y_{d}(t_{i})-Y_{m}(t_{i},\vec{q})|}{M_{diag}}.$

Here $Y_{d}(t_{i}),Y_{m}(t_{i},\vec{q})$ are smoothed New Diagnoses with
confirmed COVID-19.

The New Diagnoses dataset for validation model in the Novosibirsk region was
from 2020.03.12 to 2020.11.01. The produced forecast covered 45 days and was
compared against statistical data for the given test period (2020.11.02 –
2020.12.15). The modeling results shown in Fig. 24 demonstrated that
preserving the containment measures introduced in the region could have
reduced the spread by more than 25%. This difference between modeling and real
data was due to the data space’s stationarity (see Section 2.2.2) amid
unaccounted administrative measures.

Figure 24:  Forecast for New Diagnoses validated with historical data for the
Novosibirsk region. The black line displays the modeled period, and the red
dashed line – the test period. The red dashed line marks the test period. The
shaded areas are 10% and 90% quantiles for 10 simulations

Figure 25 presents the results of a 30-day forecast for the Novosibirsk region
with a restored vector of unknown parameters $\vec{q}$ as is described in
Section 4. The colored dots mark the real data (to 2021.02.15) that were used
to solve the inverse problem; the black ones (2021.02.16 – 2021.03.01) – to
test the forecast. The forecast for New Diagnoses assumed the rate of daily
tests remained unchanged. It is noteworthy that despite the forecast error for
the New Diagnoses indicator participating in the minimization of target
functional $J(\vec{q})$ was 10%, the forecast error for Num Critical did not
exceed 5%.

Figure 25: Model calibration results for 10 simulations + a 30-day forecast
for New Diagnoses and Num Critical in the Novosibirsk region. The shaded areas
are 10% and 90% quantiles, the solid line – the median value of modeling
result, and dots – real data. The dashed vertical lines are COVID-19
containment measures.

The effective reproduction number $R(t)$ calculated by (5) for a time period
from 2020.03.12 to 2021.02.01 in the Novosibirsk region can be seen in Fig.
26. The graph demonstrates the pandemic had been under control since the end
of 2020.

Figure 26: Effective reproduction number $R(t)$ for the Novosibirsk region.
The red dashed line marks $R(t)=1$.

##  6 Conclusions and discussions

Mathematical models are effective tools to deal with the time evolution and
patterns of disease outbreaks. They provide us with useful predictions in the
context of the impact of intervention in decreasing the number of infections
and deaths.

This paper provides a detailed analysis of the statistical data, including the
number of tested, positive, mortality, hospitalization and critical cases, on
COVID-19 spread in NY State, the UK and the Novosibirsk region and presents a
solution to the problem of identification of unknown epidemiological
parameters (transmissibility; the initial number of infected individuals;
probability of being tested, etc.) in an ABM. The problem has been considered
as the minimization of a target functional in relation to daily numbers of
tested, positive and mortality cases in the studied regions and become an
important modification of the Covasim package [21]. The minimization problem
has been solved using the Optuna hyperparameter optimization software and the
Parzen estimation method. Normal gradient descent methods do not work with
Covasim or other agent-based models, due to the stochastic variability between
model runs that makes the landscape very “bumpy” (i.e., many transient local
minima). One way of getting around this is to perform many different runs and
take the average. However, averaging over many runs is computationally
expensive, since running $N$ simulations of ABM will only reduce the noise by
$\sqrt{N}$ [23].

The results of data analysis in every studied region showed the weekly
seasonality of New Tests, which helped us forecast the future values of this
time series. It worth noting that some countries have their own rules of
statistical analysis that have to be accounted for when carrying out modeling,
e.g. in the USA, the New Diagnoses indicator contains a certain percentage of
probable cases, while New Deaths in the UK accounts for all the deaths that
have occurred within 28 days since a positive COVID-19 test, even if such
death has not been provoked by the virus.

Due to the high sensitivity of the transmissibility parameter $\beta$, the
accuracy of its identification becomes crucial for uncovering the pattern of
COVID-19 spread in an investigated region. However, not all containment
measures affect the pattern. For that reason, when developing an ABM
calibration algorithm based on epidemiological data, we paid special attention
to transmissibility and the times this parameter changed while modeling. These
characteristics were determined from solving the minimization problem as a
piecewise - constant function, while solving the inverse problem restored the
parameter together with its times of change.

The devised ABM has been validated with historical data from 2020. The
modeling results for the three regions in question have demonstrated that
preserving the introduced containment measures would have sustained New
Diagnoses in NY State and the Novosibirsk region during March 2021 and would
have reduced them in the UK.

Due to the specific features of the observed data in the Novosibirsk region,
where two time series are stationary with probability of 1, the model’s
forecast accuracy has been lower for new cases of infection and more accurate
for cases of hospitalization, e.g., the confidence interval values for the new
case of infection registered daily was [70; 270] individuals and for
hospitalizations - [90; 300] individuals, both including as model solutions as
real data.

The proposed agent-based model has the following limitations. We do not fix
population number changing during model year, consider waning of immunity to
coronavirus, or the possibility of re-infection. We also use more simplistic
contact structures than in real life.

Our future plans are investigation of model identifiability to real data and
sensitivity analysis. It will also be necessary to investigate the influence
of vaccination on COVID-19 propagation.

##  7 Competing interests statement

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

##  8 Acknowledgements

The data analysis part (section 2) is supported by the Russian Foundation for
Basic Research and Royal Society (project no. 21-51-10003). The agent-based
mathematical model construction and analysis of numerical results (sections 3,
4, 5) was supported by the Russian Science Foundation (project no.
18-71-10044) and the Royal Society $IEC\backslash R2\backslash 202020$ –
International Exchanges 2020 Cost Share between UK and Russia. Authors would
like to thank Professor Sergey Kabanikhin and Doctor Alexey Romanyukha for the
problem analysis and fruitful discussions. Also we would like to thank Yan
Reznichenko for language help.

## References

  * [1] I.Cooper, A.Mondal and C.G. Antonopoulos: A SIR model assumption for the spread of COVID-19 in different communities. Chaos, Solitons and Fractals 139(110057) (2020)
  * [2] F. Ndairoua I. Area, J.J.Nieto and D.F.M.Torresa: Mathematical modeling of COVID-19 transmission dynamics with a case study of Wuhan. . Chaos, Solitons and Fractals 135(109846) (2020)
  * [3] Kai D., Goldstein G.F., Morgunov A., Nangalia V., Rotkirch A.: Universal masking is urgent in the COVID-19 pandemic: SEIR and agent-based models, empirical validation, policy recommendations. arXiv preprint arXiv:2004.13553. (2020). 10.13140/RG.2.2.21662.08001
  * [4] Wolfram C.: An Agent-Based Model of COVID-19. Complex Systems 29(1), 87-105 (2020). 10.25088/ComplexSystems.29.1.87.
  * [5] Cuevas E.: An agent-based model to evaluate the COVID-19 transmission risks in facilities. Computers in biology and medicine 121(103827) (2020). https://doi.org/10.1016/j.compbiomed.2020.103827.
  * [6] Chang S., Pierson E., Koh P.W., Gerardin J., Redbird B., Grusky D., Leskovec J.: Mobility network models of COVID-19 explain inequities and inform reopening Nature 589, 82–87 (2021). https://doi.org/10.1038/s41586-020-2923-3.
  * [7] Silva P.C.L., Batista P.V.C., Lima H.S., Alves M.A., Guimarães F.G., Silva R.C.P.: COVID-ABS: An agent-based model of COVID-19 epidemic to simulate health and economic effects of social distancing interventions. Chaos, Solitons and Fractals 139(110088), (2020). https://doi.org/10.1016/j.chaos.2020.110088
  * [8] Mellacher P.: COVID-Town: An integrated economic-epidemiological agent-based model. Munich Personal RePEc Archive Paper No. 103661, (2020).
  * [9] Tracking Coronavirus in New York: Latest Map and Case Count https://www.nytimes.com/interactive/2021/us/new-york-covid-cases.html
  * [10] RBC site - coronavirus statistics in Novosibirsk region: https://nsk.rbc.ru/
  * [11] The COVID Tracking Project in USA: https://covidtracking.com/data
  * [12] Hodrick R.J., and Edward C. P.:“Postwar U.S. Business Cycles: An Empirical Investigation.” Journal of Money, Credit and Banking, vol. 29, no. 1, pp. 1–16. (1997)
  * [13] Bartolucci F., Cardinali A., Pennoni F.: A Generalized Moving Average Convergence/Divergence for Testing Semi-strong Market Efficiency. Mathematical and Statistical Methods for Actuarial Sciences and Finance. 10.1007/978-3-319-89824-7_18. (2018).
  * [14] Grimstad, B. and Knudsen, B.R. Mathematical programming formulations for piecewise polynomial functions. J Glob Optim 77, 455–486 (2020).
  * [15] Coronavirus (COVID-19) in the UK: https://coronavirus.data.gov.uk/
  * [16] Dabral P.P. and Murry, M.Z: Modelling and forecasting of rainfall time series using SARIMA. Environ. Process. 4, 399–419 (2017).
  * [17] Balasmeh O., Babbar R. and Karmaker T.: Trend analysis and ARIMA modeling for forecasting precipitation pattern in Wadi Shueib catchment area in Jordan. Arab J Geosci (2019).
  * [18] Box G.E.P. and Cox D.R.: An analysis of transformations. (With discussion). J. Roy. Statist. Soc. Ser. B 26, 211–252 (1964).
  * [19] Williams, T. Adaptive Holt—Winters Forecasting. J Oper Res Soc 38, 553–560 (1987).
  * [20] Kerr C., Stuart R., Mistry D., Abeysuriya R., Hart G., Rosenfeld K., Selvaraj P., Nunez R., Hagedorn B., George L., Izzo A., Palmer A., Delport D., Bennette C., Wagner B., Chang S., Cohen J., Panovska-Griffiths J., Jastrzebski M., Oron A., Wenger E., Famulare M. and Klein D.: Covasim: an agent-based model of COVID-19 dynamics and interventions. Medverix (2020). 10.1101/2020.05.10.20097469.
  * [21] Kerr C., Rosenfeld K., Hagedorn B., Mistry D. and Klein D.: COVID-19 trends in Oregon: Preparing for opening up. Institute for Disease Modeling (2020).
  * [22] Kerr C., Hagedorn B., Mistry D. and Klein D.: COVID-19 trends in Oregon: Implications for interventions. Institute for Disease Modeling (2020).
  * [23] Covasim documentation: https://docs.idmod.org/projects/covasim/en/latest/index.html
  * [24] C.M. Hazelbag, J. Dushoff, E.M. Dominic, Z.E. Mthombothi and W.Delva: Calibration of individual-based models to epidemiological data: A systematic review. PLOS Computational Biology (2020).
  * [25] OPTUNA: hyperparameter optimization framework: https://optuna.org/
  * [26] Rasmussen C.E. and Williams C.K.I.: Gaussian Processes for Machine Learning (2006).
  * [27] J.Bergstra and R.Bardenet, Y.Bengio et al.Advances in Neural Information Processing Systems. Curran Associates Incorporation. 24 (2011)
  * [28] Household Size, 2019, UN. https://population.un.org/Household/#/countries/840