key: cord-0705999-ggr6alga
authors: Hritonenko, Natali; Yatsenko, Olga; Yatsenko, Yuri
title: Model with transmission delays for COVID‐19 control: Theory and empirical assessment
date: 2021-11-23
journal: J Public Econ Theory
DOI: 10.1111/jpet.12554
sha: 969fc6760b653a08398a423c8677fc3dc187faf9
doc_id: 705999
cord_uid: ggr6alga

The paper focuses on modeling of public health measures to control the COVID‐19 pandemic. The authors suggest a flexible integral model with distributed lags, which realistically describes COVID‐19 infectiousness period from clinical data. It contains susceptible–infectious–recovered (SIR), susceptible–exposed–infectious–recovered (SEIR), and other epidemic models as special cases. The model is used for assessing how government decisions to lockdown and reopen the economy affect epidemic spread. The authors demonstrate essential differences in transition and asymptotic dynamics of the integral model and the SIR model after lockdown. The provided simulation on real data accurately describes several waves of the COVID‐19 epidemic in the United States and is in good correspondence with government actions to curb the epidemic.

satisfactorily explain the dynamics of many real-world epidemics (see Getz & Dougherty, 2018; Lloyd, 2001 , and the references therein). To describe realistic infectiousness distributions, multicompartment epidemiologic models have been suggested and analyzed (Arino & Portet, 2020; Champredon et al., 2018; Lloyd, 2001) . These models give more flexibility in describing infection dynamics, though they are described by larger systems of ordinary differential equations (ODEs). At the best of the authors' knowledge, no systematic research has been done to analyze the time-varying transmission rate in such models.

An alternative way to accurately describe the realistic infectiousness distributions in continuous time is provided by integral and integrodifferential equations (Brauer et al., 2019; Breda et al., 2012; Kermack & McKendrick, 1927) . Integral dynamic models have been known for many decades and successfully used in various applications (Anita, 2000; Boucekkine et al., 1997; Hritonenko & Yatsenko, 2005 Jovanovic & Yatsenko, 2012; Volterra, 1930) . They open new possibilities to describe and optimize the economic control of epidemics. Such models regularly appear in mathematical epidemiology, but their study is mostly focused on stability and bifurcations (see, e.g., Cooke et al., 1996; Hethcote, 2000; Hritonenko & Yatsenko, 2005; McCluskey, 2010) . Breda et al. (2012) and Champredon et al. (2018) compare integral epidemiological models to ODE-based compartmental models and demonstrate that they produce the same dynamics in special cases. All known integral epidemic models are theoretical and not approbated on real data. Epidemiologists commonly use discrete statistical techniques (Cori et al., 2013; Flaxman et al., 2020 , and others) that incorporate clinical infectiousness distributions of specific diseases, which are difficult to include in ODE-based models (Koyama et al., 2021) but easy in intergal models. As shown below, the integral models are more general and allow for a deeper analysis and control of related processes.

The goal of this paper is a thorough study of connections between government actions to control COVID-19 spread and resulting dynamics of the epidemic. To do this, we offer a new, flexible epidemic model of COVID-19 spread and mitigation. The constructed integral model with distributed lags depicts various patterns of infectivity and latency, including those obtained in clinical data. We prove that the constructed model contains SIR, SEIR, and other wellknown epidemic models as its special cases. It accurately describes delays in the variable transmission rate to capture government control of COVID-19.

Next, we focus on the formal description of government controls (lockdown and reopening the economy) as changes in the transmission rate. A review of current modeling efforts to optimize controlled COVID-19 economics demonstrates that those models oversimplify a critical element, the reaction function, that connects the transmission rate to the effectiveness of government actions to mitigate epidemic. Specifically, recent economic models of COVID-19 control (Acemoglu et al., 2020; Atkeson, 2020; Garriga et al., 2020; Gori et al., 2021, and others) assume that the transmission rate possesses a discrete jump at the instant when a lockdown begins. The delayed nature of the infection process becomes critical for adequate modeling when the transmission rate sharply changes. We prove that both asymptotic and transition dynamics of the integral and SIR models differ after the lockdown time.

To get insight from real data, we calculate the time-varying transmission rate in the constructed model on COVID-19 data in the United States over 2020-2021. It leads to an ill-posed Volterra integral equation (VIE) of the first kind, which is solved using regularization techniques. Our simulation demonstrates that the calculated transmission rate does not have a jump at lockdown times as Acemoglu et al. (2020) , Atkeson (2020) , and Garriga et al. (2020) presume. The observed changes in the transmission rate and reproduction number match the timeline of US public health actions (including three lockdowns) to mitigate COVID-19 spread.

The analysis of simulation outcomes allows us to assess the effectiveness of previous lockdowns but does not provide enough data to construct the reaction function.

We do not pursue an optimal lockdown problem in the constructed model with distributed lags, because its solution may have a little policy value in the absence of reliable reaction function and cost estimates.

Original theoretic contributions of the paper include: (1) a novel integral model with realistic infectiousness distributions, (2) a proof of the model equivalency to the SIR, SEIR, and multicompartment epidemic models in special cases, and (3) theoretic analysis of transmission rates in the integral model and the standard SIR model after a lockdown.

The paper is as follows. Section 2 introduces and analyzes the integral epidemic model with an accurate COVID-19 infectiousness distribution and its connections to known epidemic models. Section 3 discusses formal descriptions of lockdown as a controlled transmission rate, calculates this rate for the US COVID-19 epidemic in 2020-2021, and uses simulation results to evaluate government actions to curb the epidemic. Section 4 concludes.

As in the classic SIR model of Kermack and McKendrick (1927) , we divide a population of N individuals into three stages: S (susceptible), I (infectious), and R (recovered or dead). The sizes S (t), I(t), and R(t) of those stages at time 0 ≤ t < ∞ are governed by the following VIEs with a finite distributed delay:

where T is the maximal duration of the infectiveness period (the memory of a dynamic process). The given model functions are the transmission rate b(u, t − u), t − T < u < t, at which one individual infected at time u contaminates others at time t, and the fraction of individuals f(t − u), u < t, infected at time u and still infectious at time t. These functions reflect two different processes:

• The infectiousness intensity b(u, s) ≥ 0 that depends on the contact time u and the time s passed after the infection occurs, b(u, s) = 0 for s > T. • The survivor distribution, that is, the probability f(s) ≥ 0 that an infectious individual remains in the stage I after s days being infected, f(0) = 1, f(s) = 0 for s > T.

The function b determines how suspected individuals become infected, while the function f defines how the infected individuals recover. Equation (2) defines the number of suspected individuals that become infected (moved to the stage I) at time t. Equation (1) means that the percentage f(t − u) of the individuals infected at time u remains in the stage I at time t, while 1 − f(t − u) is moved to the stage R.

Using standard techniques for VIEs, one can see that Equations (1)-(3) have a unique nonnegative solution I(t), S(t), R(t), 0 < t < ∞, at the given initial conditions

As in most epidemic models (Arino, 2020; Arino & Portet, 2020) , we suggest N = const and do not consider birth or natural death, only the death by removal from the infectious stage (in contrast to endemic models). All recovered individuals are assumed to become immune, so, the only susceptible fraction S/N of the population can become infected. The fraction 0 < δ(t) < 1 of "recovered or dead" is assumed to die, so, the total number of deaths is

The death rate δ is considered as a constant fixed parameter in many models, but it fluctuates in reality (Our World in Data, 2020).

Integral dynamic models have been commonly used for the description of population dynamics, including epidemiology, recall, for example, Volterra and Lotka integral models. Cooke et al. (1996) , Hethcote (2000) , McCluskey (2010) , and others have studied the global stability of such models with a constant transmission rate. Breda et al. (2012) discuss the integral epidemic model

with respect to S(t) and the force of infection F(t − s) that corresponds to b(t − s, s)I(t − s) in (1)-(3). The model (6) does not consider the dependence of F on the infection time s. Integral models with finite memory have been used in other applications (Boucekkine et al., 1997; Hritonenko & Yatsenko, 2013 . As we will show, the finite memory T < ∞ makes the epidemic model (1)-(3) more realistic. Models with finite memory have also appeared in epidemiologic modeling for the purpose of global stability analysis.

The unique features of the model (1)-(3) are:

• it clearly distinguishes two important delay processes: varying infectiousness b(u, t − u) and the recovery of infected f(t − u) after the infection occurs, • it describes any realistic distribution of infectious periods (Lloyd, 2001) , • it contains other known epidemic models as special cases.

In this section, we temporarily assume T = ∞ to compare our model to its well-known ODE alternatives. Then, the integral model (1)-(3) can be rewritten as

is equivalent to the integral model (7)-(9) at

Proof. The proof follows from differentiating the integral equations (7)-(9) at the assumptions (14). Indeed, taking the derivative of (7) leads to

the derivative of (8) is given by (12), and the derivative of (9) is

The theorem is proven. □ Theorem 1 reflects the well-known fact that infectious individuals in the SIR model (10)-(13) move to the stage R at a constant Poisson rate γ > 0 (Hethcote, 2000) . This model describes the so-called memoryless dynamic system, in which the probability to stay in stage I is the same regardless of how long the individual was infected. This means that some individuals will remain in stage I at any distant time in the future, which is impossible. In reality, all individuals exit the stage I in a finite time after becoming infected.

Compared with the SIR model, the SEIR model includes the separate stage E (exposed), which contains individuals in a latent state of infection (that are infected but not yet infectious). The exposed individuals E eventually become sick and move to the Infectious stage at the exponential rate σ > 0. The model is described by four nonlinear ODEs (Hethcote, 2000) :

Proof. To prove this result, we first solve the linear ODE (16) and obtain an analytic expression for E(t). Next, substituting E(t) to the ODE (17) and solving it, we obtain the expression for I(t) that coincides with the integral equation (1) 

The theorem is proven. □

The SEIR or susceptible-latent-infected-removed (SLIR) model 2 was a step ahead in modeling of diseases with latent periods, such as COVID-19. However, later comprehensive research (Arino & Portet, 2020; Getz & Dougherty, 2018; Lloyd, 2001) concludes that the SEIR model does not satisfactorily match epidemiological data.

To describe more realistic infection distributions, multicompartment epidemiological models have been constructed and analyzed (Champredon et al., 2018; Lloyd, 2001) . Such models employ the so-called "linear chain trick" (using additional compartments that do not exist in reality) to obtain better epidemiologic models in the ODE framework (Champredon et al., 2018; Getz & Dougherty, 2018; Lloyd, 2001 ).

Following Arino and Portet (2020) , Champredon et al. (2018) , Getz and Dougherty (2018) , and Lloyd (2001) , a multicompartment modification of the standard SIR model (10)-(13) divides the "Infectious" stage into n substages I 1 , …, I n and is described by the system of n + 2 nonlinear ODEs:

In this model, newly infected individuals enter the first substage I 1 , pass through each successive substage, and move to the "Recovered or dead" stage as they leave the nth substage I n . The variables S, R and parameters β, γ have the same meanings as in the SIR model (10)-(13). Champredon et al. (2018) show that the probability of individuals to remain in the stage I in the model (21)-(25) is described by the Erlang distribution:

where F(τ; n, γ) is the regularized gamma function with parameters n and γ.

The distribution (26) is determined by its rate γ and the integer shape parameter n. It is close to the normal distribution, except for small values of n. As the number n of compartments increases, the distribution becomes more closely centered around its mean γ. The limit of (26) at n → ∞ is the Dirac delta distribution δ(γ), which corresponds to all individuals having the same infectiveness duration 1/γ. Arino and Portet (2020) construct a SLIAR model with additional compartments for COVID-19 epidemics.

Theorem 3. The multicompartment SIR model (21)-(25) is equivalent to the integral model (7)

Proof. To prove this result, we sequentially solve the linear ODEs (21)-(25) starting with i = 1. □ Champredon et al. (2018) analyze a multicompartment Erlang SEIR model, which divides both E and I stages of the SEIR model (15)-(18) into substages, and numerically demonstrate that it produces the same epidemic dynamics as the epidemic integral model

with respect to S and the incidence i = bSI. However, their analytic expressions for g(s) are extremely complex, so, we do not analyze them here. Instead, we focus on other patterns of infectivity and latency not covered by ODE-based models. As we will see in Section 2.2, such patterns for COVID-19 are neither exponential nor gamma distributed.

The key natural feature of the integral model (1)-(3) is the assumption T < ∞ about the finite delay (a finite length of process prehistory). In reality, all epidemic processes have a finite length that depends on the nature of disease. Ignoring this fact in both ways (assuming no memory or infinite memory) makes models less accurate. All infected individuals recover or die, which makes the infinite memory assumption ridiculous. The finite memory of processes also leads to new phenomena in modeled dynamics. The integral equations (1)-(3) with finite delay have a unique solution I(t), S(t), R(t), 0 < t < ∞, at the given initial conditions (4). The initial problems (1)-(4) are solved sequentially on the intervals [(k − 1)T, kT], k = 1, 2, 3, …. The delay equations possess a richer dynamic picture than the ODEs. In particular, depending on given functions I˜0 and S 0 in (4), the delay equations (1)-(3) may have oscillating and discontinuous solutions. To consider only continuous solutions of (1)-(3), here and thereafter we assume that I˜0 satisfies (1) at t = 0.

The distributions b and f are given key components of the model (1)-(3), which define its flexibility and applicability to real epidemics. The correct choice of b and f for COVID-19 spread depends on fundamental features of a region or country, such as the population size and density, economy, transportation, habits, and others.

The distribution b(u, s) in (1) describes the infection intensity depending on two independent instants: the contact time u and the time s passed after the infection occurs. It is a critical assumption of the model (1)-(3). The dependence of the infection intensity b on the time-sinceinfection s is a major concern in modern epidemic models (Arino, 2020; Martcheva, 2015) . All comprehensive epidemic models with additional compartments (see Section 2.1.3) were designed to describe more accurate infectiousness period distributions from clinical practice.

Let us assume a fixed infection instant u = u 0 and denote β(s) = b(u 0 , s), where s = t − u 0 ≥ 0 is the time-since-infection. The dependence of b(u, s) on the infection time u will be explored in Section 3. In theory, the infectiousness period distribution β(s) may be any distribution, from exponential to normal distribution or uniform distribution (Lloyd, 2001) . In practice, it should be determined on a finite interval [0, T] from detailed clinical data. Following recent clinical research (van Kampen et al., 2021) , the dependence β(s) for COVID-19 is approximately described by the solid black curve in Figure 1 and is used in Section 3 simulation. Such a shape of the dependence of infectivity on the time-since-infection is common for influenza-like diseases (see Figure 13 .1 in Martcheva, 2015) . For comparison, Figure 1 contains the SEIR distribution (26) and the Erlang distribution (27) at n = 2 and 5. The distribution β(s) in Figure 1 reflects the latent period of COVID-19 pandemic: symptoms become visible 4-5-days after infection. In our simulation, we assume the infectiousness interval of T = 14 days, which includes both latent and infectious periods of infected individuals.

The most important feature of the integral model (1)-(4) is that it can use any arbitrary infectiousness distribution β(s). We shall note that all published epidemic models (until 2020) have not seriously analyzed the dependence of b(u, s) on the infection time u. It appears to be very relevant in Section 3 for description of current public efforts to control the pandemic.

The survivor distribution f(s), s ∈ [0, T], describes the fraction of individuals that remain infectious s days after becoming infected. It should be also based on detailed clinical data, if available. It is natural to assume that f(s) decreases from f(0) = 1 to f(T) = 0. The integral model (1)-(3) can use an arbitrary distribution f. Below we provide several special cases of the survivor distribution f with a finite memory T > 0.

Triangular survivor distribution

means that the constant part 1/T of infected individuals is released to the stage R each day.

Exponential partial distribution

Dependence of infectivity on time-since-infection for COVID-19 clinical data, SEIR distribution, and Erlang distribution at n = 2 and 5

In this extreme case, all infected individuals remain infectious exactly T days after becoming infected and then they immediately move to the stage R.

combines the triangular and uniform distributions. All infected individuals remain infectious T 1 days after becoming infected (reflecting the latent period) and then gradually move to the stage R.

We shall note that the further analysis and simulation in Section 3 demonstrate that the shape of the distribution f affects modeling outcomes less seriously as compared with the infectiousness period distribution β.

Since 2020, epidemiologic models from the SIR family have been intensively used to describe the economic control of the COVID-19 epidemic without addressing the models' suitability (Acemoglu et al., 2020; Alvarez et al., 2020; Atkeson, 2020; Caccavo, 2020; Fernández-Villaverde & Jones, 2020; Gollier, 2020; Gori et al., 2021) . We argue that the specifics of COVID-19 epidemic control require more accurate analysis and modeling.

This section analyzes formal ways to incorporate government actions to control epidemics in economic-epidemiologic models. A necessary step is to assess how the government's actions and preventive measures affect the epidemic's spread. We consider the period of March 2020-March 2021, when widespread vaccination was not available and, therefore, public health measures to reduce transmission of COVID-19 were restricted to so-called nonpharmaceutical interventions (NPIs; see Peak et al., 2017 , and the references therein). NPIs are supported by government actions and include compulsory stay-at-home policies (lockdowns), bans on public gatherings, closing/opening schools and nonessential businesses, face mask orders, quarantine and sanitary zones, and similar. The effectiveness of NPIs has been studied in historical data about the influenza pandemic (Fong et al., 2020; Hatchett et al., 2007) . The consensus of such studies is that NPIs, quickly implemented after detecting a new contagious virus, can essentially reduce the virus's transmission.

In this section, we attempt to identify and formalize the effectiveness of the most common NPIs, the lockdown and its recall, in terms of the decreasing numbers of infected people. Obviously, NPIs carry high costs and cause significant economic losses that must be properly assessed in more detailed economic models.

Numerous economic modeling efforts of COVID-19 epidemic control started in 2020. The focus of those efforts is to describe how the public health NPIs and government actions change the transmission rate in the ODE-based epidemic models. Thus, epidemic-economic models (Acemoglu et al., 2020; Alvarez et al., 2020; Atkeson, 2020; Caccavo, 2020; Fernández-Villaverde & Jones, 2020; Gori et al., 2021) use the time-varying transmission rate β(t) in standard SIR and SEIR models to capture some public health NPIs, specifically, lockdown and related behavioral changes.

To highlight the time-dependent transmission rate in the integral model (1)-(3), here and thereafter we assume that the infectiousness intensity b(u, s) is the product of two

where u is the time when infection occurs, s is the time passed since the infection occurs. The function β is described in Section 2.2.1. Then, Equation (1) becomes

Reproduction number: The most publicly known epidemiologic parameter is the reproduction number R t , which is the expected number of new infections generated by one infected individual. This parameter is frequently used in public health decision-making and by top government officials. Thus, the UK Prime Minister Boris Johnson regularly refers to the reproduction number to describe COVID-related goals and actions of the UK government (BBC News, 2020).

The number R t is directly related to the transmission coefficient β but calculated differently in different models. In SIR model (10)-(13), it equals

At R t > 1, the number of infected in the SIR model increases, otherwise, it decreases and eventually extinguishes. The basic reproduction number

describes the expected number of new infections generated by one infected person in a totally susceptible large population, that is, when I(t) « N and S(t)~N.

The reproduction number in the integral population model (1)

(see Hritonenko & Yatsenko, 2005 Iannelli & Milner, 2017) . At the beginning of epidemic, when I(t) « N, the basic reproduction number is

Using (33), we can show that then the number of infected in the model (1)-(3) increases at R t > 1 and decreases at R t < 1, provided that initial functions I˜0 and S 0 in (4) are monotonic.

Several SIR-based economic models of COVID-19 (Acemoglu et al., 2020; Alvarez et al., 2020; Garriga et al., 2020; Gollier, 2020) suggest that β(t) possesses a discrete jump at the time t l when the government introduces a shelter-in-place order to lock down a fraction of susceptible and infected individuals. 3 Namely, the controlled transmission rate is described by the reaction function

where the function 0 ≤ L(t) < 1 is the lockdown intensity, and a fixed parameter 0 < θ < 1 is the measure of lockdown effectiveness (Alvarez et al., 2020; Stock, 2020) . In these models, the lockdown control jumped instantly (overnight) from 0 to a certain value L, which can be described by the step function

Exponential parameterization of the variable transmission rate is the alternative hypothesis to (37). Fernández-Villaverde and Jones (2020) applied the SIR model (10)-(13) to extensive data about COVID-19 deaths in different countries and found that discrete jumps at the lockdown time fit the data worse than a smooth evolution of β(t). Caccavo (2020) uses the exponential parametrization of β(t) after the lockdown time:

characterized by its initial value β 0 > 0, a final value β 1 > 0 at infinity, and the decay rate λ > 0. Fernández-Villaverde and Jones (2020) used the "exponential decay model" (38) in the initial version of their paper but dropped it later in the favor of recovering the time-varying β(t) from observations (what we do in Section 3.4). Atkeson (2020) uses a combination of two exponents to obtain a U-shaped function β(t) (to describe both lockdown and reopening the economy). However, any such parametrization has obvious weaknesses and cannot describe more detailed features of the process. In general, the variable transmission rate leads to specific delays in disease propagation that cannot be adequately described in the ODE framework. Section 3.3 proves that such delays essentially affect the transition and asymptotic dynamics of epidemics after a change in β(t). 

We compare two alternative models of the time-varying transmission rate, jump (37) and exponential parameterization (38), in the integral model (1)-(3) and standard SIR model (10)-(13). For analytic clarity, we restrict our analysis to the initial epidemic stage, when the proportion of Infected is small: I(t) « N. Then, both models are reduced to linearized equations for I.

If I(t) « N, then S(t) = N − o(N), so, the SIR Equation (11) can be formally written as I'(t) = β(t)I (t)(1 − o(N)/N) − γI(t), which is well approximated by the equation

The linear ODE (39) with the jump (36)-(37) at t = t l has a piecewise exponential solution I t I e t t I t e t t ( ) = at 0 < , ( ) at < < .

Typical dynamics of (40) is presented in Figure 2 with a dashed line. The solution of the linear ODE (39) at the exponential β(t) defined by (38) is

F I G U R E 2 Dynamics of I(t) in linearized models (1)-(3) and (10)-(13) with lockdown HRITONENKO ET AL.

The curve (41) is shown in Figure 2 with a gray solid line and is smoother than (40).

At I(t) « N, the delay equation (30) is well approximated by the linear VIE of the second kind

known as the renewal equation (Champredon et al., 2018; Hritonenko & Yatsenko, 2013 . For simplicity, we assume the exponential kernel f(t) = exp(−st), β = 1. Then, the integral equation (42) reduces to the linear ODE with lump delay

with the initial condition I(t) = I 0 (t) at −T ≤ t ≤ 0. Analytic solution of Equation (42) ).

x s T −( + )

The rate x has the following properties:

The linearized SIR equation (39) at β(t) = const has the exponential solution I t I e ( ) =

The proof is in the appendix. Thus, both linearized models at a constant β possess exponential solutions though with different growth rates and reproduction numbers. Next, we prove that the lockdown (37) causes quantitatively and qualitatively different dynamics in both models. 

and have the following properties

The proof is in the appendix.

Following Acemoglu et al. (2020) , Alvarez et al. (2020) , and Fernández-Villaverde and Jones (2020), we take the values β 0 = 0.2 and γ = 0.067 in the SIR model (10)-(13). Let the lockdown (36) occur at t l = 11. Choosing L = 0.7 and θ = 0.75 as in Acemoglu et al. (2020) and Alvarez et al. (2020) leads to the post-lockdown β 1 = 0.045 in the SIR model (10)-(13). To adjust the exponential solution of delayed Equation (43) to the SIR solution I(t) = I 0 exp((β -γ)t) before lockdown, we select a smaller β 0 = 0.155 in (43). Then, the rate x ≈ 0.133 in Lemma 1 is close to β 0 − γ. For clarity, we take the uniform f(t) ≡ 1 in (42). Figure 2 compares the dynamics of three linearized models: (1)-(3), SIR (10)-(13), and the SIR with β parametrization (38). The infection curve I(t) in all three cases is exponential before the lockdown t l = 11 and possesses a corner point at t l = 11. However, the further dynamics of I (t) is quite different in those models.

To analyze further dynamics of I(t), we rewrite (43) at s = 0 as initially, but it decreases and becomes negative β 1 − β 0 < 0 at t* such that t Τ t Ι( *− ) = Ι( *). Later, I t I t ′( )/ ( ) jumps up from β 1 − β 0 Ι(t − Τ)/Ι(t) < 0 to = β 1 (1 − Ι(t − Τ)/Ι(t)) < 0 at t = t l + T. So, the (43) solution has the second corner point at t = t l + T = 26, clearly visible in Figure 2 . It describes a new phenomenon in the integral model (1)-(4), the infection echo caused by a sharp change in the transmission coefficient β at the lockdown time t l = 11. Such effects are common in delay equations arising in economics, population biology, and other applied sciences (Boucekkine et al., 1997; Hritonenko & Yatsenko, 2013 .

As Theorem 4 and the above simulation reveal, the difference of the rates y and β 1 − γ of both models is much larger after lockdown. Indeed, the rate y = −0.076 is 3.5 times smaller than β γ − = −0.22 1 . This fact greatly affects the asymptotic of upcoming stages of an epidemic. Thus, I(t) decreases three times faster in the model (1)-(3) at t > t l + T as compared with the SIR model. By Figure 2 , I(t) in the integral model is 10-15 times smaller after 60 days compared with the SIR model. Thus, a distributed delayed nature of the infection becomes critical for adequate modeling when the transmission rate β(t) sharply changes.

The goal of this section is to determine the unknown transmission rate β(t) from the integral model (1)-(3), assuming that the variables S(t), I(t), and R(t) are known over a long enough interval t T T [ − , ] 0 m a x . In the integral model (1)-(3), such analysis requires solving a VIE of the first kind to identify the variable transmission coefficient. We develop a regularization algorithm for solving this identification problem and approbate it on the COVID-19 data for the United States in 2020-2021. It allows us to assess the epidemic response to related US government and public health actions of that period.

Data set description: The COVID-19 data for the United States over the time horizon of February 2020-March 2021 are taken from the reliable source (Our World in Data, 2020). "Daily new confirmed COVID-19 cases per million people (the rolling 7-day average)" is used for I(t) and shown in Figure 3 . The infectiousness interval of T = 14 days includes both latent (3 days) and infectious periods. We use the infectiousness distribution β s( ) for COVID-19 and the uniform survivor distribution f(s) from Section 2.2. The 2-week delayed "Total confirmed COVID-19 cases (per 1 M)" minus "Daily confirmed COVID-19 deaths (per 1 M)" is used as a proxy for R(t).

Because of the COVID-19 latent period, β s( ) = 0 at 0 ≤ s ≤ τ, where τ = 3. Then, the integral equation (30) for β(t) can be rewritten as

The linear VIE of the first kind (47) belongs to the category of ill-posed problems, which means that small variations of given data can lead to large changes in unknown variables. There are several proven approaches how to solve ill-posed problems, known as regularization algorithms (Kalitkin, 1978; Lamm, 2000) , that could be used for solving Equation (47).

The shape β s( ) creates additional challenges in solving the ill-posed identification problem (47). Namely, a classic technique to numerically solve VIEs of the first kind is to differentiate the equation, which usually leads to a well-posed VIE of the second kind (Hritonenko & Yatsenko, 2013; Volterra, 1930) . However, our Equation (47) still remains the VIE of the first kind after differentiation:

because the given β τ( ) = 0. Equation (48) contains the derivatives I t ′( ) and β s f s (˜( ) ( ))′, which makes its numeric solution less accurate than (47), see more in Lamm (2000) about this case. For this reason, we use a numeric algorithm based on direct discretization of Equation (47) using a proper numeric integration rule. The algorithm includes two parts.

Part 1. Discretization and regularization of Equation (47) The discrete functions S(t), I(t), and R(t) are known over the interval t T T [ − , ] 0 m a x with the discretization step h = 1 day. Assuming that the independent variable t is an integer and approximating the integral equation (47) using the trapezoid rule, we obtain the system of linear algebraic equations

where C 0 = C T = 1/2, C i = 1 for i = 2, 3, …, T − 1, and N = 10 6 (according to the used data). By Figure 1 , β s ( ) = 0 at s = 0, 1, 2, 3 and β s( ) > 0 at s ≥ τ + 1 = 4. Introducing the function K s β s τ f s τ ( ) =~( + ) ( + ), s = 1, …, T − τ, and replacing the independent integer variable

From (50), we obtain the explicit formula for calculating the unknown β(t):

where t = t 0 -τ + 1, …, T max , and α > 0 is a regularizing parameter. In the recurrent formula (51), β(t) depends on β(u), u = t − 1, …, t -T + τ. After finding β(t), we use the discretized equations (2) and (3) to calculate R(t) and S(t).

The discretized formula (51) is still ill-defined as its continuous analogue (48). The value K(0) > 0 in the denominator is small. First numeric calculations by (51) at α = 0, showed that the obtained β(t) alternated between large negative and positive values because of oscillations in the given I(t). To obtain a meaningful solution, we have combined three smoothing and regularization techniques: (i) adding a small regularizing term αβ(t) to the LHS of the integral equation (50), (ii) restricting calculated values β(t) to an interval [B min , B max ]: 0 < B min ≤ β (t) ≤ B max , and (iii) using 7-day central moving averages. Varying the values α, B min , and B max as regularization parameters of the identification problem, we have obtained an acceptable solution for β(t), shown in Figure 3 . The criterion for choosing regularization parameters (Kalitkin, 1978) is the visual control of differences between I(t), t = t 0 + T, …, T max , calculated below in Part 2, and the observed curve I. The corresponding reproduction number R t is calculated using (33) and shown in Figure 3 with solid black line.

Part 2. Verification via the original integral model Because of the regularization used in Part 1, we must verify that the obtained transmission rate β satisfies the original model (1)-(3). To do so, we substitute the solution β(t) at t = t 0 , …, T max obtained in Part 1 and shown in Figure 3 to the discretized equations (1)-(3) and calculate the discrete functions S(t), I(t), and R(t) for t = t 0 + T + 1, …, T max , where t 0 = March 4, 2020 and T max = March 4, 2021. In doing so, we use the known values of S(t), I(t), and R(t) at t = t 0 , …, t 0 + T, from Part 1 as given initial conditions on the prehistory of the process. Next, we compare the calculated S, I, and R with the given S, I, and R of Part 1.

Because of the presence of the latent period, the formulas (1)-(3) are explicit. Nevertheless, this calculation also appears to be unstable because the given initial I(t) is much smaller at the offset of epidemic. Specifically, I(t) increases 40 times from 2.5 to 95 (per 1M) during 20 days from March 4 to March 20, 2020. As a result, if t 0 is March 20, then the calculated values of I are twice larger than those in the data set. However, solving the verification problem with a delayed t 0 = May 1 demonstrates a good qualitative agreement with original data, shown in Figure 4 , with the average difference of 3%-5% between the given and calculated values of I.

This section relates changes of the reproduction number R t to corresponding actions of the US local and federal governments. A nationwide lockdown is not possible in the United States because of its multilevel government structure. Nonetheless, many US states have taken extreme public health measures to control COVID-19 spread during three waves of COVID-19 in 2020-2021 (Mervosh & Lee, 2020; USA Today, 2021; Washington Post, 2020 , and similar sources), which are described below.

Uncontrolled growth of COVID-19 infection occurred at the beginning of 2020. The number I(t) of infected increased 40 times from 2.5 to 95 (per million) during 20 days from March 20 to April 18. In the model (1)-(3), the identified reproduction number R t reached 10 on March 24 (such large initial values of R t were also reported for other countries, specifically, India, China, and France).

On March 15, the US Center for Disease Control and Prevention (CDC) issued the first guidance for the 8-week period. Most US schools were closed starting March 16 through at least April 20. Bars and restaurants were also closed. Starting March 18, most states had implemented stay-at-home orders (on average, from March 28 to May 9). However, five rural states never imposed stay-at-home orders favoring other restrictions. All those and many other local actions led to a sluggish decrease of the Infected number and the corresponding reproduction number R t < 1 starting April 5 (Figure 4 ).

In May, many states lifted stay-at-home orders, reopened nonessential businesses, and eased social distancing measures intended to curb the virus' spread. It was too early and caused a spike in the number of Infected, which happened immediately after the calculated reproduction number R t became greater than 1 at the end of June. This relation is clearly visible in Figure 5 that shows the period May 2020-April 2021 on a larger scale.

Responding to spikes in coronavirus cases, state governments took urgent actions. Beaches and bars reopened for Memorial Day were shut back down for the fourth of July. Masks were required in many states and on all airlines. On July 25, at least 14 states imposed new restrictions in places that serve food and drink or revived those recently lifted. The number of infected reached a (local) maximum of 200 per million on July 25 (when R t became 1) and decreased immediately after that. Some states jumped ahead of the CDC guidelines on safe reopening. So, the number I(t) of infected quickly increased in October 2020. The reproduction number was around 1 during October 4-11 and reached 1.5 on November 2. To combat the second wave, several states imposed statewide restrictions in November, closed bars and restaurants for indoor dining, limited or prohibited indoor social gatherings, closed museums, casinos and movie theaters, and limited the capacity of most stores to 25%-50%. Most colleges and schools returned to allremote instruction. Many states adopted curfews, closing businesses at certain hours.

As a result, the number I(t) reached its local maximum on December 19, briefly went down, and increased again after Christmas break. The reproduction number was 1 or above from December 26 to January 15, 2021. In early January, indoor social gatherings were banned, indoor dining was prohibited, and indoor entertainment establishments were limited in many states and cities. Some state governments extended stay-at-home advisory. The reproduction number went down to 1 on January 15 and decreased fast to 0.7 later on. The number of Infected was more than 700 per million on January 5-12, 2021, and then decreased four times in 1 month.

Simulation based on the model (1)-(3) satisfactorily describes the dynamics of COVID-19 in the United States that included three waves of infection over the 1-year period. It also demonstrates sophisticated dynamics of the time-varying transmission coefficient β(t) and reproduction number R t . Observed changes in the reproduction number R t are explained by related government actions.

F I G U R E 5 The calculated reproduction number (a bold black curve) R t > 1 when the observed number I(t) of infected (in gray) increases, and R t < 1 otherwise. The transmission rate β(t) is shown for comparison with a thin curve Formally, our final goal is finding a reaction function for the controlled transmission rate β, which would describe β dependence on government controls and preventive measures. A simple example of such solution is the dependence of β(L) = β 0 (1 − θL) 2 on the lockdown effectiveness L, given by (35). Unfortunately, this dependence is not very realistic. As seen in Figures 3-5 , the time-varying transmission coefficient β(t) does not follow the step function (37).

Three approximate periods with intensive epidemic mitigation because of related public health measures (mostly, lockdowns) are highlighted in Figures 3-5 . Let us denote the effectiveness (severity) of those three lockdowns as 0 < L i < 1, i = 1, 2, 3. In Figures 4 and 5 , both β(t) and R t fluctuate due to systematic data reporting errors and biases (Bergman et al., 2020) and the ill-posed nature of the identification problem 4 . After the removal of short-term fluctuations associated with weekly data errors, the dependences β(L i ) can be approximated by decreasing linear functions.

Provided simulation demonstrates the applicability of the integral epidemic model (1)-(3) for description and control of COVID-19 pandemic. The authors aim to apply the developed analytic and simulation technique of COVID-19 epidemic to the UK and other countries. It will provide additional information for more accurate formalization of the dependence of controlled transmission rate β on government controls and preventive measures.

Dozens of recent publications estimate the time-varying reproduction number of COVID epidemics in the United States, France, China, and other countries, based on different models (e.g., Caccavo, 2020; Fernández-Villaverde & Jones, 2020; Fitzgibbon et al., 2020; Pang et al., 2020) . Some of them relate the dynamics of COVID-19 spread to economic actions and public health measures. Below we compare simulated dynamics of the reproduction number R t (33) with estimates by two distinct modeling approaches, Chudik et al. (2021) and Koyama et al. (2021) . The resulting curves are shown in Figure 6 , where the observed number I(t) of newly infected is shown with a thin line. Chudik et al. (2021) empirically estimate the effects of social distancing combining the classic SIR model (10)-(13) with agent-based modeling. They compute the time-varying transmission rate β(t) minimizing the 2-week-averaged squared difference between the observed stochastic rate and its exponential approximation. This technique eliminates short-term fluctuations. As shown in Figure 6 , the obtained dynamics of R t over March 2020-March 2021 is similar to the one for model (10)-(13). However, unlike our integral model, it does not accurately match the observed dynamics of I(t) (the rolling 7-day average). Namely, the calculated R t > 1 on some intervals, such as April 15-May 10, August 1-14, where the observed I(t) decreases. Similarly, R t < 1 on September 17-30 and August 1-14, where the real I(t) increases. The likely cause is too fat right tail of the survivor distribution f u e ( ) = γu − in the SIR model (see Theorem 1). We expect that the application of Chudik et al. (2021) algorithm to the integral model (1)-(3) would produce more precise results.

Well-known discrete statistical techniques such as Cori et al. (2013) and Flaxman et al. (2020) estimate the time-varying reproduction number using clinical infectiousness distributions. Compared with ODE-based models, the integral model (1)-(4) is closer to this approach because of the finite length T of infectiousness period. Flaxman et al. (2020) and Koyama et al. (2021) analyze the effects of NPIs on COVID-19 spread using a Bayesian mechanistic model with a finite infection cycle. However, they impose too strict assumption that the timevarying reproduction number R t is a piecewise-constant function that varies only when an intervention occurs. As result, their calculated dynamics of R t for the United States over 2021 in Koyama et al. (2021) does not satisfactorily describe changes in the observed number I(t) shown in Figure 6 . It is quite different from our Figures 3-5 and other sources.

An ultimate modeling goal is to identify effective governmental actions to mitigate COVID-19 and support economy. All COVID-related papers attempt to address this goal to a certain degree. Current modeling efforts to optimize controlled COVID-19 economics (Acemoglu et al., 2020; Atkeson, 2020; Garriga et al., 2020) assume the objective of minimizing the total costs that include production losses, cost of treatment, loss of life, and other damages. A simplified mathematical outline is as follows. Two major losses from the epidemic can be described as:

(a) the economic loss due to the lockdown L(t) at time t 

where χ is the relative weight of the objective (53) compared with (52). Justifications of the cost of death parameter χ in (54) are offered by Alvarez et al. (2020) and Acemoglu et al. (2020) , who relate it to the well-known statistical value of life. Next, the optimal lockdown control L(t), t ∈ [0, T), that minimizes (54), is estimated and interpreted. Numerous enhancements of the above economic outline appeared in 2021, particularly, in the JME special issue (Boucekkine et al., 2021) . Such traditional economic framework can lead to interesting theoretic insights and is potentially beneficial for COVID-related policy analysis. However, it should be essentially improved before applying to real data. The use of infinite horizon H = ∞ in the objective (54) looks inappropriate because of the speed and severity of economic, health, and social changes during the real time of emergency. It is arguable whether the number of deaths can be discounted in (53). In reality, recovered individuals are also included in the economic description (52) of lockdown cost. More similar issues can be identified.

Still, a major drawback is an oversimplified impact of economic controls on epidemic dynamics. Indeed, Acemoglu et al. (2020) , Atkeson (2020) , and Garriga et al. (2020) consider a simple government control, the lockdown (36), and assume that it causes an immediate jump (37) in the transmission rate. In contrast, our simulation on real US data demonstrates more complex delayed impact of lockdown on the transmission rate, which cannot be approximated by a single jump.

Modeling of the optimal government lockdown policies has a little practical value unless they are feasible. It requires a preliminary economic assessment and mathematical formalization of specific control actions, such as border shutdowns, travel restrictions, quarantine, self-isolation, lockdown, and social distance. Robust methods need to be developed for an assessment of the effectiveness of governmental controls and preventive measures in pandemic conditions. Actual governmental policies should be explored before formalizing them in mathematical models. The future research should extend identification techniques to advanced control actions, such as reopening businesses, repeated lockdowns, travel restrictions, border shutdowns, self-quarantine, contact tracing, and so forth.

Corresponding models should link epidemic spread and control, economic dynamics, population health, and social habits (Avery et al., 2020; Bloom et al., 2018; Eichenbaum et al., 2020; Gori et al., 2021) . Capturing specific characteristics of a region/country such as the level of government centralization, population density, or public transportation is also important (Fernández-Villaverde & Jones, 2020). The authors are extending a similar analysis to COVID-19 epidemics in the UK and France that have more centralized governments and where lockdowns were implemented on a federal level as opposed to the United States.

Optimal targeted lockdowns in a multigroup SIR model

A simple planning problem for COVID-19 lockdown (NBER Working Paper 26981)

Analysis and control of age-dependent population dynamics

Mathematical epidemiology in a data-rich world

A final size relation for epidemic models

A simple model for COVID-19

What will be the economic impact of COVID-19 in the US?

Policy implications of models of the spread of coronavirus: Perspectives and opportunities for economists

R number for UK below 1 for first time since

Oscillations in U.S. COVID-19 incidence and mortality data reflect diagnostic and reporting factors. mSystems

Epidemics and economics: New and resurgent infectious diseases can have far-reaching economic repercussions

The economics of epidemics and contagious diseases: An introduction

Replacement echoes in the vintage capital growth model

Mathematical models in epidemiology

On the formulation of epidemic models (an appraisal of Kermack and McKendrick)

Chinese and Italian COVID-19 outbreaks can be correctly described by a modified SIRD model

Equivalence of the Erlang-distributed SEIR epidemic model and the renewal equation

Covid-19 time-varying reproduction numbers worldwide: An empirical analysis of mandatory and voluntary social distancing

Analysis of an SEIRS epidemic model with two delays

A new framework and software to estimate timevarying reproduction numbers during epidemics

The macroeconomics of epidemics

Estimating and simulating a SIRD model of COVID-19 for many countries, states, and cities (NBER Working Paper

Analysis of a reaction-diffusion epidemic model with asymptomatic transmission

Estimating the effects of non-pharmaceutical interventions on COVID-19 in Europe

Nonpharmaceutical measures for pandemic influenza in nonhealthcare settings-Social distancing measures

Optimal management of an epidemic: Lockdown, vaccine, and value of life

Discrete stochastic analogs of Erlang epidemic models

Cost-benefit analysis of age-specific deconfinement strategies

COVID-19 epidemic and mitigation policies: Positive and normative analyses in a neoclassical growth model

Public health interventions and epidemic intensity during the 1918 influenza pandemic

The mathematics of infectious diseases

Optimization of harvesting age in integral age-dependent model of population dynamics

Bifurcations in nonlinear integral models of biological systems

Mathematical modeling in economics, ecology and the environment

Nonlinear integral models with delays: Recent developments and applications

The basic approach to age-structured population dynamics. Models, methods and numerics

Investment in vintage capital

A contribution to the mathematical theory of epidemics

Estimating the time-varying reproduction number of COVID-19 with a state-space method

Optimal control of prevention and treatment in a basic macroeconomic-epidemiological model

A survey of regularization methods for first-kind Volterra equations

Realistic distributions of infectious periods in epidemic models: Changing patterns of persistence and dynamics

An introduction to mathematical epidemiology

Complete global stability for an SIR epidemic model with delay-distributed or discrete

See which states are reopening and which are still shut down. The New York Times

Daily confirmed COVID-19 deaths, rolling 3-day average

Transmission dynamics and control strategies of COVID-19 in Wuhan

Comparing nonpharmaceutical interventions for containing emerging epidemics

Coronavirus data gaps and the policy response to the novel coronavirus (Discussion Paper 20-82)

Map of COVID-19 case trends, restrictions and mobility

Duration and key determinants of infectious virus shedding in hospitalized patients

Theory of functionals and of integral and integral-differential equations

Model with transmission delays for COVID-19 control: Theory and empirical assessment

The authors are thankful to Associate Editor Dr. Raouf Boucekkine for valuable comments that essentially improved the clarity of the paper. The paper is partially supported by PVAMU FEP and the Ministry of Education and Science of Kazakhstan (Grant AP09261118).

The authors declare that there are no conflict of interests.

The data that support the findings of this study are available from the corresponding author upon reasonable request. 

To justify that Equation (A1) has an exponential solution, which is the same over the prehistory [−T, 0] and the interval [0, ∞), we substitute I˜(t) = I 0 exp(xt) to (A1) and obtain the nonlinear equation (44) for the rate x. Next, let us rewrite Equation (44) in the nondimensional form aswith respect to the new unknown variable z = (x + s)/β and parameter ρ = βΤ = β/γ, and introduce the notationis concave downward, monotonically increasing from −∞ to 1, and passing the point (0, 0). Thus, Equation (A2) has the trivial solution z = 0 and may have a nontrivial solution z*. Applying the Taylor series to (A2) in the neighborhood of (0, 0),, which implies that if ρ > 1, the straight line y = z is below the curve y e = 1 − zρ − at small 0 < z « 1 and intersects the curve later because F z ( ) 1 → . It means that the nontrivial solution in this case is positive z* > 0, or, in the original notations, x + s > 0 at β > γ. Analogously, if 0 < ρ < 1, the line y = z is below the curve y e = 1 − zρ − for small z < 0. Thus, z* < 0, or x + s < 0 at β < γ. A straightforward analysis of (A2) shows that z → 0 at ρ → 1, z → 1 at ρ → ∞, and z → − ∞ at ρ → 0. In terms x and s, we obtain the properties x + s → 0 at γ → β, x → β − s at γ → 0, and x → − ∞ at γ → ∞.The lemma is proven.

At conditions of the theorem, β(t) = β 0 for 0 ≤ t ≤ t l . Therefore, the delay differential equation (43) on the interval 0 ≤ t ≤ t l before lockdown is the ODETo find an exponential solution to (A4), which is the same on the prehistory [−T, 0] and the interval [0, t l ), we substitute I˜(t) = I 0 exp(xt) to (A4) and obtain the nonlinear equation (20) for the rate x. The properties of its solution x follow from Lemma 1.The delay differential equation (43) at t l ≤ t ≤ t l + τ becomes the ODEUsing the formula for general solution of the linear ODE, the solution to (A5) is the combination of two exponents:At t > t l close to t l , I˜(t) increases but I˜′(t) becomes 0 at some t max , and then it decreases while t < t l + T.Finally, the delay differential equation (43) at t l + τ < t < ∞ becomesThe influence of the given prehistory diminishes for large t » t l + T and the solution of the linear equation (A7) asymptotically tends to an exponent. Substituting I˜(t) = I 1 exp(yt) to (A7) at large t, we obtain the nonlinear equation (45) for the rate y. The properties (a) and (b) are proved similarly to the analysis of the nonlinear equation (A2) in the proof of Lemma 1.The theorem is proven.