1. Introduction. In the course of battling COVID-19, public health policies sought to enforce non-pharmaceutical interventions to slow or halt the spread of the pandemic. Common policies included 'safer-at-home', 'social distancing' and 'mask wearing' mandates, which were seen as crucial during the early stages of the pandemic prior to the availability of vaccines. The timeline of COVID-19 globally and locally ([7, 26]) indicates that the evolution of policy affected the evolution of the pandemic and vice versa. For example, in the county of Los Angeles, social distancing was first mandated [9] on March 21, 2020, about a month after the first reported COVID-19 case in LA. Around that time, the Los Angeles Mayor's Office released the 'safer-at-home' policy [1]. One week later, beaches, hiking trails, dog parks, skate parks, and other public sites and facilities were temporarily closed. On April 15th, as infected cases continued to increase, facial coverings were mandated in many indoor places [19]. In hindsight, it is important to ask: Were policies that were enforced done so in an optimal way? What can we learn by using mathematical modeling to understand the interplay between policy and spread of disease? This paper introduces a policy model coupled to a susceptible-infected-recovered (SIR) epidemic model to study interactions between policy-making and the dynamics of epidemics. There have been several studies on the relationship between policies and epidemics [4,5,21,6,18]. In a study analyzing data from 16 US cities during the 1918 pandemic [5], Bootsma and Ferguson analyzed specific outcomes related to the impact of the delay of lockdown policies on the total deaths and also on the appearance of second waves of outbreaks due to reopening too early. The analysis was done fitting available data to an SEIR model. They also considered optimal control for the simpler SIR and the end-state of the pandemic, noting that there exists an optimal control level with fewer deaths and no second wave. More recently, Bliman et al. that can change continuously in time, which would imply, for example, the ability to shift in three successive instants between no restrictions, perfect "lockdown", and back to no restrictions. As observed during the COVID-19 pandemic, policies that change frequently in time cannot be easily followed. Moreover, policies must be relatively easy to interpret, with a small number of different intensity levels (see Fig. 3b). A practical implementation also requires a minimum time duration for a particular stage of the policy. These practical constraints can be modeled together as a piece-wise constant function of time with a minimum time interval for each well-defined policy level (i.e., not continuous). With this idea in mind, we aim to re-examine the optimal practical policy among all possible piece-wise constant policies with minimal time duration.

Second, Bliman et al. assume that that the only outcome to manage is the final epidemic size. This so-called "impact cost" is clearly a central concern (see below). However, as also seen during the pandemic, there are real trade-offs between decreased infections and the negative impact of strict policies on other aspects of society such as remote learning for young students, employment curtailment in certain job sectors, and lack of key services provided to the public. In the present work, we modify Bliman et al.'s model to take into account these other practical "implementation costs." Specifying a short minimal time interval during which policies must remain constant (e.g., one week), we find our results resemble Bliman et al.'s bang-bang controller [3] despite the more complex cost structure that includes both impact and implementation costs. With a larger minimal time interval during which time policies must remain constant (e.g., 28 days), optimal policies depart from the bang-bang solution.

Finally, Bliman et al. also assume a pandemic spreading in an single population pool overseen by a single policy-making entity. The reality of the COVID-19 pandemic is that there are policy makers at several (nested) hierarchical scales that oversee different population pools. For example, within the United States, policies may be set at Federal, State, County and local levels, not to mention finer-grained institutional and family scales. And populations at any one scale (e.g., counties) may interact to varying degrees. Inspired by the work of Jia et al. [14], we introduce a hierarchical version of Bliman et al.'s model with sequential (Stackleberg) policy-making. Specifically, levels higher in a jurisdictional hierarchy make policy decisions, while levels lower in the hierarchy make their decisions with full knowledge of the policy recommendations from above. We find that a hierarchical structure can make the policies converge in all regions using the right weight for a non-compliance cost.

The remainder of this paper is organized as follows. We first introduce the work in [3] and reproduce the results using our methods. We discuss how different optimal policies result from different parameter choices for model constraints and costs. Next, we discuss an empirical case study of the so-called "second wave" of the pandemic (November 6th 2020-May 12th 2021) in Los Angeles County, California.

Last, we use simulation to study optimal control of the pandemic in three counties with and without a governing state as an example of the multi-layer multi-regional model.

2. Policy model using optimal control. A policy function is a continuous function that has a range of [0, 1]. As the numerical value increases, the strictness of the policy decreases. The Numerical value 0 denotes a total lockdown and 1 denotes no control. We assume a policy u(t) directly influences the level of a lockdown, which affects the rate of the population transport from compartment S to I. We use the following policy-incorporated SIR:

Like the traditional SIR model, the reproduction number R 0 = β γ . Herd immunity occurs when a large proportion of the population has become immune to the infection. Mathematically, it is defined as the value of S below which the number of infected decrease and can be calculated as S herd = N R0 . In [3], a policy u(t) is assumed to belong to the admissible set U αmax,T0 defined by

The constant T 0 characterizes the duration of the policy, and α max its maximal intensity. In [3], Theorem 2.1 states that no finite time intervention is able to stop the epidemics before or exactly at the herd immunity. However, one may stop arbitrarily close to herd immunity by having a sufficiently long intervention of sufficient intensity. To determine the closest state S to this threshold attainable by control of maximal intensity α max on the interval [0, T 0 ], one is led to consider the following optimal control problem:

(2) Furthermore, Bliman et al. prove the existence and uniqueness of the optimal solution to problem 2 and that the solution is a bang-bang controller (a control that switches from one extreme to the other). More specifically, they have the following theorem:

with respect to α max and non-decreasing with respect to T 0 .

(ii) there exists a unique

1 [T0,+∞) (in particular, the optimal control is bang-bang).

3. Single region case. We use the same policy-incorporated SIR model for the epidemic dynamic as in [3]. Instead of minimizing the final epidemic size alone, we adopt a similar policy-making process as in [14] by using a cost function that takes into account the cost of implementing the policy, the impact of the infection and a penalty for being non-compliant. The latter cost only applies in hierarchical models where a lower-level unit can choose to not follow the policy recommendation of a higher-level unit.

We also consider practical implementation constraints, namely that the policy can only be implemented using a finite number of discrete levels of control and with a minimal time interval during which a policy must remain constant. As an example, consider the policy implementation in France during the year 2020 and 2021 shown in Fig. 1 ( [27]). Implemented policies were discrete both in terms of the small number of intervention types and the fixed time intervals of enforcement, the shortest of which was approximately 15 days in duration, with the longest lasting more than a year. A discrete policy model is realistic given the empirical pattern of real-world interventions. Such a model also simplifies the computation problem of optimal policy discovery by searching through a discrete set of potential policies rather than a continuum of policies.

3.1. The policy-incorporated SIR model. To model the evolution of the pandemic, we discretize the system of ODE using forward Euler's method with a time (i) the implementation cost, which represents the consequences of policies meant to curtail the pandemic on individuals and the broader economic and social systems.

(ii) the impact cost, which represents the consequences of people getting sick both on individuals and the broader economic and social systems.

(iii) the non-compliance cost, which is a penalty imposed by a policy-maker upon an agent within its jurisdiction for deviating from its recommendation (e.g., a fine or litigation costs).

The implementation cost is a non-increasing function of α and the impact cost function is a non-decreasing function of α. The coefficients κ, η, κ + η ∈ [0, 1]. The cost from time t 1 to t 2 is defined as the averaged integral of the cost function over a total time period T :

There are different ways to parameterize the cost function. In this paper, the cost function is parameterized in the following way:

where R t2 (u) is the fraction of the recovered population at time t 2 if policy u is adopted during [t 1 , t 2 ] and π(u) is the policy of the agent one level above. The parameterization of the implementation cost and the non-compliance cost are adopted from [14]. The impact cost is parameterized as the recovered population at time t 2 to approximate the impact on the medical system since a fraction of the recovered represents the hospitalized population. If the cost function u is fixed at constant value α over time interval [t 1 , t 2 ] , the cost can be written as:

An example of cost functions with different weights using the above parameterization is shown in Fig. 6. In our simulation for a single region, we use a averaged total cost over a time period T as the following:

If at time T , the SIR model has reached the equilibrium, we can use R T (α) to approximate R ∞ , the fraction of the final size of the recovered population. To find the optimal policy, we solve for the following optimization problem:

3.3. Algorithm. We discretize time by MPTI ∆t and the policy intensity into multiple levels. Let T be the total time and A be the set of possible policy intensities (e.g., A = {0, 0.5, 1}). We search for all the policies that lead to S final being close to S herd , i.e. S final > S herd -ϵ, for some sufficiently small ϵ using a depth-first search algorithm [23]. The depth-first search algorithm stores the cost up to the current time interval and reuses this result to obtain the total cost for each policy function through backtracking. Let N = T ∆t and N denote the number of stages of a policy.

In total, there are |A| N policies. We initialize the minimal cost c min to be 9999.

Assume the initial susceptible and infected population are S 0 , I 0 , respectively. For n-th time interval (n < N ), we choose a value from the set intensity levels A that has not been used before, calculate the cost for the policy intensity, add it to the previous cost, and calculate the susceptible and the infected at the end of n-th time interval using the chosen intensity. Then we move to (n + 1)-th time interval. If the end time interval is reached, we check if S final > S herd -ϵ. If yes, we calculate the cost for the final time interval and add it to the previous cost to get the current total cost c. If the total cost c is smaller than c min , we update c min with the total cost c, and the optimal policy u opt with u. Next, we go back to the previous time interval and repeat the same procedure. After searching over all policies, the policy with the lowest cost is the optimal policy. The detailed algorithm is presented in Alg. 1.

Algorithm 1 Single-region policy SIR 1: Input: Time T , initial infected population I 0 , initial susceptible population S 0 , intensity levels A, minimal policy time interval ∆t, policy end time T 0 , Tol ϵ 2: Initialize county policies, minimal cost c min = 9999, current cost c = 0

for intensity level α ∈ A do 10: the general cost function (8) reduces to the impact cost and is parameterized as the final epidemic size R ∞ . Bliman et al. assume that the paths considered all reach herd immunity. Therefore, in our search for the optimal policy, we exclude cases that do not reach herd immunity. Note that without this exclusion, the optimal solution is to adopt and hold the strictest possible policy starting from the beginning of the pandemic. This results in the least number of infections. For ease of computation, we consider three levels of policy intensity: 0, 0.5, 1 and fixed time intervals for the MPTI. We use the same set of parameters for the SIR model as in Bliman et al.

[3]: N = 6.7 × 10 7 , I 0 = 10 3 , S 0 = N -I 0 , R 0 = 2.9. Following [3], we also choose the policy end time T 0 as close as possible to 100, thus setting T 0 = 98 since the time interval needs to be a multiple of the MPTI of 7 days. We show the result our algorithm produces in Fig. 2a which we visually compare to the result from [3],

shown in Fig. 2b. Note that we normalized curves by the total population. Both solutions are bang-bang controllers. The solution using our model starts the control on day 63 (a multiple of 7) rather than day 61.9 (continuous). Slightly more people are infected under a policy that is forced to use seven day intervals compared with continuous time as used by Bliman et al.

Using a larger minimal policy time interval of 28 days and T 0 = 112, the optimal solution is no longer a bang-bang controller, as shown in Fig. 2c with a larger S ∞ = 0.32. The optimal policy starts with a looser "intermediate" policy phase followed by a stricter phase. Interestingly, in practice, during COVID-19 it was common for policies to start with the strictest restrictions followed by partial opening [27,9]. Thus, it is interesting to contrast the optimal policy with a policy in which the two stages are flipped in time, see Fig. 2d. The flipped policy is a sub-optimal solutionit results in a larger pandemic size and a second wave of infections, as was often seen during the first two years of the COVID-19 pandemic. Nevertheless, the policy in Fig. 2d, while infecting more people, divides the impacted population into two distinctive waves, which could decrease daily hospital demand over the course of the outbreak. Our policy model does not optimize for hospital demand. Since many public health agencies (including Los Angeles County) considered hospital demand when making policy decisions, it could be important to consider in future studies.

6.7 × 10 7 10 3 N -I 0 2.9 0.296 B 100 Not applicable 6.7 × 10 7 10 3 N -I 0 2.9 0.31 C 112 28 6.7 × 10 7 10 3 N -I 0 2.9 0.32 D 112 28 6.7 × 10 7 10 3 N -I 0 2.9 0.174

Table 1. Parameters.

1 u(i) represents the i-th entry of vector u.

has generated a wealth of models and results using mechanistic approaches taking explicitly into account the movement of individuals ( [13,15,22]). For example, in [22], the authors proposed a multi-regional compartmental model using medical geography theory (central place theory) and studied the effect of the travel of individuals (especially those infected and exposed) between regions on the global spread of severe acute respiratory syndrome (SARS). Another way to account for the interplay between regions is to use a cross excitation matrix [28]. This scheme assumes the a uniform mixing of the population across regions and the infected population in one region can trigger the infection in another. The entries of the matrix records the pair-wise cross excitation from one region to another. In this paper, we assume uniform mixing in the population and use an excitation matrix K = {K aa ′ } to model the travel and infections across counties. Our network-style SIR is the following:

For any county a, the rate of change from S a to I a triggered by I a ′ depends on

where κ f + η f = 1 and R f,T (α) is the number of the recovered of region f which is an aggregation of the the epidemic size of its leaf nodes. 4.1. Algorithms. The single region algorithm minimizes over all admissible piecewise functions, while the multiple-region algorithm only minimizes over every time interval. We assume there are up to three layers: federal government, the states, and the counties. At n-th time interval, we first determine the optimal policy intensity that minimizes the cost C f n∆t,(n+1)∆t for the federal layer. After obtaining the optimal federal policy, each state optimizes its own cost function C s n∆t,(n+1)∆t

for the period [n∆t, (n + 1)∆t] unilaterally, i.e., assuming other states follow their previous policies. Next, we choose the optimal policy intensity for the counties in the same manner. Note that the federal layer does not pay the non-compliance cost as it is not subject to any higher-level policy making. The states and counties may pay a non-compliance cost. The full details of the three-layer model is in Alg. 2. 4.2. Simulations. In this section, we present results for a three-county example of the multiple regions game and a three-county example with a state. First, we discuss when one layer exists (i.e., only counties).

The game between the counties 0.0 0.2 0.4 0.6 0.8 1.0 Policy intensity α 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Cost κ=0.33, η=0.33 κ=0.5, η=0.25 κ=0.25, η=0.75 κ=0.8, η=0.1 κ=0.99, η=0.01 Figure 6. Different cost functions vs policy intensity α. Algorithm 2 Game Policy SIR 1: Input: Time T , excitation matrix K, intensity levels A, time interval ∆t 2: Initialize state, county policies 3: Number of policy stages N = T ∆t , n = 1 4: while n ≤ N do 5: t = n∆t 6: while t < T do 7: for every state s do 8:

for every county a in state s do 9: update S a , I a , R a according to the current policy α a and the excitation matrix K: S a (t) = S a (t -1) -α a β a ′ K aa ′ I a ′ (t-1)Sa(t-1) Na 11:

R a (t) = R a (t -1) + γI a (t -1)

13: end for 14: end for 15: t += 1 16: end while 17:

for every state s do 19:

for every county a in state s do 21:

end for

end for

n += 1 25: end while is through cross excitation of infection among the counties. Next, we study the case 1 when a governing state is added.

We consider three interacting counties with the excitation matrix K:

We set the reproduction number R 0 = 2 and therefore, S herd = 0.5. Counties 1, 2, 3 have initial fractions of the infected population as i 0 = 0.2, 0.1, 0.1, respectively. This implies that county 1 has a bigger outbreak initially, and part of the infection in county 2 is excited from county 1 and part of the infection in county 3 is excited from county 2. The cost functions for all counties consist of an implementation cost and an impact cost with equal weights (η a = κ a = 1/2, for all a). The minimal policy time interval ∆ is set to be 7 (days).

The left column (Figs. 7a, 7c, 7e) are simulations for the counties without any intervention and the right column (Figs. 7b, 7d, 7f) are simulations with interventions. Without intervention, we see propagation of waves of infection from county 1 to county 2 and then to county 3. All of the counties reached herd immunity eventually. With interventions, policy restrictions started on day 7 and, for county 2 and 3, the infected curves decrease before reaching their peaks. With control, county 1 contained the pandemic and the final S ∞ is close herd immunity level S herd . With a fewer infected population to begin with, county 2 and 3 contained the pandemic before reaching herd immunity. Fig. 8 shows the results of adding a governing state on top of the county layer. We keep the ratio of the weights for the implementation cost and the impact cost to be 1:1, the same as in the no-state case in Fig. 7. The state has slightly different weights, with the ratio of the weights for the implementation cost and the impact cost being 1:2. Compared to Fig. 7, by adding a state, the three counties ended up with the same policy. In this case, the noncompliance cost results in each county choosing the same policy as the state rather than different policies. In the search for an optimal policy, we used a naive depth-first search algorithm for the one-region model. One can speed up the algorithm by removing some of the obvious non-optimal paths. In our model, the policy intensity α is a heuristic representation of the lockdown, social distancing and mask policy. It remains to be discussed how other policies, for example, vaccination policies, affects the spreading in the different stages of the pandemic. The model ignores some of the important features like the limitation of the hospital capacity [24], which could be added as constraints when minimizing the cost function. Fig. 3b shows the policy for the first wave is proactive while the one for the second wave is reactive. One possible effect is from fatigue of following policy, which increases in time and has a memory. So far, the model does not have the capability of modeling this fatigue. In the future, one could consider an adaptive term in the cost function to model it. The network example considered was rather simplistic, with just three counties within one state. One could consider more complex systems with multiple layers. The computational method here would likely