Optimal control of pandemics via a sociodemographic model of non-pharmaceutical interventions

1.   Introduction

The human population is no stranger to a modern epidemic, with resurgences of Dengue fever throughout hundreds of years of history ^[47], repeated Ebola outbreaks since 1976 ^[21], the AIDS epidemic beginning in 1981 ^[19], and more recently the COVID-19 pandemic, which took center stage in the year 2020 ^[42] and has been a unifying motivation of researchers across all scientific fields. With each surge of virus comes a desire to understand both its biology and transmission. Different viruses can be transmitted differently; via direct physical contact, airborne droplets produced from coughing, sneezing, etc ^[32], and in some cases via a nonhuman vector population ^[44]. Understanding the transmission route is imperative to developing a mitigation strategy. For example, in the case of the West Nile Virus, it is more important to control the human interactions with the vector population than to control human-to-human contact ^[14]. However, in the case of a virus like COVID-19, it is very important to understand human-to-human contact to make forecasts and ultimately policy decisions that are best for both a population and its economy^[1,18].

In this paper, we introduce a general model for studying interactions between people of different sociodemographic categories in the case of a virus where transmission occurs during a human-to-human interaction. We start by developing the model which uses an interaction matrix to govern interaction levels based on a sociodemographic characteristic. Next, an optimal control problem is considered with regards to optimizing demographic-targeted lockdowns, followed by simulations to test the model and, lastly, a conclusion with future steps.

1.1.   Epidemiological modeling of populations and optimal control

Epidemiological models representing a population in terms of interacting susceptible/infected sub-populations date back to Kermack and McKendrick who subdivided populations into compartments and provided differential equations driven by infection and recovery rates ^[27]. The model type that was born from this work is often called a "SIR" model or "SEIR" model standing for Susceptible, Exposed, Infected, and Recovered which are the labels representing the compartments of the model. The equations of such a model read $\dot{S} = -\frac{\beta S I}{N}$, $\dot{E} = \frac{\beta S I}{N}$ $-\delta E, \dot{I} = \delta E - \gamma I$, $\dot{R} = \gamma I$, where $\beta$ and $\gamma$ are the infection rate, and recovery rate respectively, with $\delta$ being the latent period of the virus, or the time between exposure and infectivity. The key parameter $R_0 = \frac{\beta}{\gamma}$ ^[15] represents the "basic reproduction number", that is, the number of secondary cases in which one case would produce in a completely susceptible population.

This basic model has been adapted in many ways by increasing the number and complexity of compartments to: capture disease progression ^[20], consider undetectable infections or asymptomatic infections ^[11,36], consider social parameters such as age-structure and spatial distribution of populations ^[8,26,53], control vaccination and other mitigation strategies ^[13,35], consider in-host dynamics ^[4,43], and more.

SIR models are used as a basic logic for many modeling efforts, including those that do not necessarily use differential equations to govern movement. With the increasing computing power available to researchers, Agent-Based Modeling (ABM), otherwise known as individual-based modeling, has gained quite a bit of interest ^[5]. In ^[22], for example, an ABM's ability to easily model spatial parameters is employed to study indoor non-pharmaceutical interventions of contagious respiratory diseases. These types of models also allow for easy parameterization of a unique population and the characteristics identified as important to viral transmission.

Model results, including forecasts, are improved by more descriptive parameters. In the case of a SEIR model, once such population-specific information is provided, and the researcher can answer questions and make forecasts based on model results. In ^[52], optimal control is implemented to measure the impact that awareness has on a population being plagued by AIDS. In ^[6] Pontryagin's Maximum Principle (PMP) is applied to a control problem for disease in livestock where the controls involve vaccination, isolation, culling, and reduction of transmission. The optimal lockdown schedule was shown in ^[49] to depend strongly on the mutation of the virus itself over a time period. In fact, when viral mutation is considered, waves of virus can match very closely to the waves seen during the COVID-19 pandemic.

1.2.   Paper's contributions

In this paper, we propose a general SEIR model that uses a sociodemographic interaction matrix to capture the interaction levels between subpopulations. The additional parameterization of subpopulation allows us to track more specifically the disease progression and allows the description of targeted non-pharmaceutical measures. The availability of specific data is usually limited, thus we focus on age-brackets to subdivide the population in groups which match well the US census data. The mathematical approach is easily adaptable to other sociodemographic variables.

Then, we formulate an optimal control problem, where the control variables are the level of severity of a finite set of possible lockdown interventions. We choose three specific interventions: general lockdown, school closures, and restriction to elderly interactions. In practice, the latter can be achieved via food delivery services, home care, and retirement homes policies. With a virus pandemic in mind, as the recent COVID-19 one, the selection of the three interventions allows to capture the importance of protecting vulnerable subpopulations (elderly had a much higher death rate) as well as understanding the role played by intense social interactions (as those happening in schools). The cost function is designed to capture the economic cost of interventions. More precisely, a first component accounts for the cost of interventions, due to production loss and increased costs, while a second component accounts for cost of deaths.

Once the model is detailed, we perform standard analysis using mathematical control theory. The existence of optimal control follows from classic results and we point out that the PMP suggests that optimal control is likely bang-bang, i.e., for most of the time, the intervention will be either not performed or performed to its maximum extent. More conclusive results are impeded by the high-dimensionality and nonlinearity of the dynamics, thus we resort to numerical methods.

Extensive simulations allows us to study both the importance of various factors on the efficiency of interventions, as well as the characteristics of optimal policies. We first compare the cost of no intervention with the optimal one for $R_0$ = 1.7 (similar to COVID-19 in some of the pandemic phases), showing an improvement of $90\%$ of the cost. Even for more aggressive viruses ($R_0$ = 2.5) the cost reduction is of $80\%$. Interestingly, the optimal lockdown policies include an initial general lockdown, restrict the elderly interactions for most of the time horizon, and include alternating closure and opening of schools. Subsequently, we show how the timing of interventions is more critical than their intensity. As expected, if limited interventions are allowed, then the optimal policy tends to extend them over time.

The cost per person of each intervention is estimated from literature, thus we include an investigation of the impact of daily costs over optimal policies. The total cost here is the sum of the costs of lockdown with the cost of death. This represents a realized economic cost of dealing with the pandemic. The daily cost per person represents the cost, each day, of locking down that person. This estimate includes things like not going to work, not using public resources, and not shopping/eating at public establishments. An estimated cost of future deaths has been included in order to deal with poor behavior of the optimizer in the final few days since the simulation ends before the final infected can pass to recovered and subsequently be included in deaths calculations. The general conclusion is that the optimal policies are very consistent for different daily costs, with optimal policies always consisting of initial general lockdowns limited in time, lengthy restrictions to elderly, and intermittent school closures. We finally look at the different interventions separately. School closures are the least effective when used as only measures, while restrictions for elderly and general lockdowns provide significant cost reductions. We also observe that general lockdowns may happen to be intermittent.

Our main conclusion is that models including sociodemographic information can be of paramount importance in highlighting key factors in pandemic management. We compute optimal policies for non-pharmaceutical interventions, and estimate the economic and social impact on the population. Future development will include the exploration of different sociodemographic variables and a wider set of interventions.

It is worth noting that when attempting to model the lockdown of a population via specific sociodemographic characteristics using equation-based modeling, things can get more complicated when overlap exists between two subpopulations. For example, if rather than age we chose to parameterize based on "beach goers" and "restaurant goers", a nonzero lockdown in both could have a double counting effect on a person that eats at restaurants and goes to the beach. Naive attempts at addressing such cases may result in unrealistic dependencies such as a model that can only limit attendance of either restaurants or beach visits, but not both at the same time. For this reason, the decisions of which variables to consider when labeling populations as well as how to incorporate such labels into the modeling effort are extremely important.

2.   A sociodemographic SEIR model

We incorporate sociodemographic variables into a classical SEIR model in such a way that a more complex mitigation strategy can be found. We augment the basic SEIR model to differentiate between subgroups of the population who may have unique interactions, as for instance, those due to school activities, and focus on age groups. However, our method may be adapted to other types of sociodemographic variables, such as geographic location, income, health access, and others, which plays an important role in pandemic progression. Our model will be fit to census data for the state of New Jersey (but could be easily adapted to other states and countries), thus we consider 7 different groups corresponding to age brackets of census data denoted by a subscript, e.g., $S_j$, $j = 1, \ldots, 7$. See Section 11 for the specific age brackets.

The interactions among different age groups will be encoded by a matrix ${l} = (l_{k, j})$, where $l_{k, j}$ quantifies the interactions that group $k$ is having with group $j$. In particular, $l_{k, j}$ will be used to weight the infections resulting from susceptibles in group $j$, i.e., $S_j$, interacting with the infected in group $k$, i.e., $I_k$. For the sake of generality, we will state our problem with a number $h\in {\mathbb N}$ of groups, while simulations will be performed using census data, thus with $h = 7$. Let us denote by $N_j = S_j+E_j+I_j+R_j$ the total population in the $j$-th age group (which is constant), and set $\widetilde{N}_j = \sum_{k = 1}^{{h}} l_{k, j}\ N_k$, which represents the total number of individuals the $j$-th susceptibles are interacting with. Then, we set:

for $j = 1, \ldots, {h}$, where $1-u$ is the lockdown rate, $\beta$ the infection rate, $l$ the interaction matrix, $\delta$ the latent period, and $\gamma$ the recovery rate. Often, as in the case of the simulations here, interaction information is available for the general population but not exactly available for the interactions between susceptible and infected people. A more severe symptomatic virus would have fewer interactions due to the realization that the infectious person is infected but a higher infection rate per interaction, while a less severe virus might have more interactions but less infectivity. One method that could be employed to capture such nuances is separating asymptomatic versus symptomatic infection into two separate compartments ^[33]. Here, we introduce the model in its most basic state, so the assumption is that the parameter choice of our infection rate and interaction matrix capture these nuances.

To study how the control of interactions amongst different groups could affect the system, we sub-divide the interaction matrix $l = (l_{k, j})$ into a well-known socially driven subset of uniquely interacting groups. We define a set of $l^i$ to be $n$ interaction matrices where $l^{i}_{k, j}$ encodes some portion of the interaction between groups $k$ and $j$ due to the $i$-th subset of interactions. For example, if there are four sociodemographic variables that are known to drive the interactions between age groups, to capture the different interactions, $n = 4$. Notice that while $l^i$ is the same size as $l$, each element is either a portion of the corresponding element of $l$ or $0$. It holds that

for every $k, j$ and, thus, $l = \sum_{i = 1}^{n}l^i$. Each new interaction matrix brings the ability to encode the regulation of that subgroup of interactions, and we will denote by $u_i$ the "amount of allowed interactions" such that the severity of lockdown for each group of interactions is $1-u_i$. Therefore, our final model reads

with $j = 1, \ldots, {h}$ and other parameters set as in Eq (2.1).

2.1.   Optimization of lockdown policies

We consider an optimal control problem for the system (2.3) representing the interactions between the various age groups. For this, we define a cost functional with two main goals; minimizing the economic cost of locking down a population or subpopulation, and minimizing the deaths that the whole population will incur over a given time horizon. Because this model has no reinfection considered, the recovered compartment is monotonically increasing and accounts for all people who have experienced the full tenure of the virus. We then can compute deaths simply by multiplying death rates with the total recovered at the end of a simulation. This total deaths is represented in the following by the variable "D".

As is customary in optimal control, we will consider a running cost ${\mathcal L}$ to be integrated over the time horizon and a final cost $\psi$ depending on the state at the end of the time horizon. Specifically, we set

with $P_i$ as the total population isolated due to the $i$-th control action (lockdown), $c_i$ is the cost of the $i-th$ lockdown per person, usually ranging between USD 10 and USD 100, and $c_{n+1}$ the cost per death, usually set to USD $1,500,000$. Additionally, $I_{\textbf{end}}$ represents the total infected at the end of the simulation. This part of the cost function takes into account the cost of deaths that will happen beyond the simulation horizon. The weight $c_{n+2}$ needs to take into account future deaths due to current and future infected and will be calibrated in simulations.

The values of $c_i$ are typically estimated in terms of the economic impact of the lockdown measure for lost production output or increased costs. For instance, isolating workers will likely decrease the production rate, while school closure will increase childcare costs. A more complex phenomenon may be measured as the loss in educational outcomes due to online school activities compared to in-class ones; see ^[16]. Many studies have shown that disparities in education between vulnerable and non-vulnerable populations grow substantially during events like school closure events ^[17]. Therefore, similar to assigning a cost of mortality, the true cost of school closure must incorporate not only a financial value but also a social value (which in the end may result in financial deficits as well).

3.   Model analysis

The proposed optimal control problem can be stated as follows:

where $f$ is the drift term, $f_i$ the controlled fields, $u = (u_1, \ldots, u_n)$ the control vector, ${\mathcal U}$ the set of admissible controls, $T$ the time horizon, $J(u)$ the cost functional, ${\mathcal L}$ the Lagrangian or running cost, and $\psi$ the final cost.

Referring to model (2.3), the system can be written in the general form Eq (3.1) with the following choices: $x = (S_j, E_j, I_j, R_j)_{1\leq j\leq 7}$, $f$ and $f_i$ are defined by Eq (2.3), ${\mathcal U}$ is the set of measurable functions $u:[0, T]\to ([0, 1])^n$, and ${\mathcal L}$ and $\psi$ are given in Eq (2.4).

Due to the specific properties of Eqs (2.3) and (2.4), we get the following:

Theorem 1. Consider the optimal control problem (3.1) with dynamics given by (2.3) and cost function by (2.4). There exists an optimal control for every time horizon $T > 0$ and initial conditions $x(0) = x_0$.

Proof. For every $i = 1, \ldots, N$, we have $u_i\in [0, 1]$, thus the control set $U = ([0, 1])^{N}$ is compact and convex. Moreover, the dynamics (2.3) is linear in the control, thus the associated multifunction $F(x) = \{f(x)+\sum_{i = 1}^n u_i f_i(x):u_i\in [0, 1]\}$ has compact and convex values.

The dynamics (2.3) is smooth in all variables and admits the natural invariant region $\{x = (S_j, E_j, I_j, R_j)_{1\leq j\leq 7} : S_j\geq 0, E_j\geq 0$, $I_j\geq 0, R_j\geq 0, S_j+E_j+I_j+R_j = N_j\}$ where $N_j > 0$; thus, for every initial condition, the support of trajectories on $[0, T]$ is uniformly bounded. Moreover, the cost function (2.4) is smooth in all variables.

We can apply Theorem 5.2.1 of ^[7], obtaining the existence of an optimal control. □

Being that the optimal control problem is linear in the controls, one may attempt to use necessary conditions, such as the PMP, (see Theorem 6.5.1 of ^[7]) to find optimal controls. However, the number of variables (7 age groups with 4 populations each) and the nonlinearity w.r.t. state variables (e.g., terms $SI$ to describe infections) render the problem too complex to obtain detailed information. Therefore, to compute optimal controls we will resort to numerical methods described in the next section. However, we provide some information derived from PMP for the general case and some results for a special case. Let us start recalling that, from the PMP, there exists a covector $(\lambda_0, \lambda)$ so that the optimal control minimizes the Hamiltonian:

Since both the dynamics (2.3) and the cost function (2.4) are linear in $u$, this condition implies $u_i\in \{0, 1\}$ unless the corresponding switching function vanishes:

where we use the fact that we can write ${\mathcal L} = {\mathcal L}_0+\sum_i u_i {\mathcal L}_i$ being the running cost linear in the control. This is indeed confirmed by simulations of Section 4, where we observe $u_i\in \{0, 1\}$ at all times. However, the number of switchings, i.e., changes from the value $u(t) = 1$ to the value $u(t) = 0$ and vice versa, is difficult to bound as observed in Figure 2. The resulting number of oscillations could be too frequent to be considered as actual policy changes over a time horizon.

Singular controls, i.e., in $]0, 1[$, may appear if the functions $\psi_i$ vanish identically on an interval. While a general analysis is out of scope, we provide a result for the simple case of one control and one age group, thus we consider the dynamics:

where we omit the dynamics of $R$ (not influencing the system). The system can be written compactly as $\dot{x} = f(x)+u_1 f_1(x)$ with $x = (S, E, I)$ and $f$, $f_1$ defined by (3.3). Singular controls may appear only if trajectories run on singular points, i.e., points $x = (S, E, I)$ such that the vectors $f(x)$, $f_1(x)$ and $[f, f_1](x)$ do not span the whole space ^[31,34]. Here, as usual, $[f, f_1]$ indicates the Lie bracket of $f$ and $f_1$; see ^[7] for details. We have the following:

Theorem 2. Consider the optimal control problem (3.1) with dynamics given by (3.3). Then, the only singular points in the first orthant satisfy either $S = 0$ or $I = 0$.

Proof. We can easily compute:

If the vectors are linearly dependent, then the matrix formed with columns $f$, $f_1$, and $[f, f_1]$ must have determinant equal to zero. A simple computation shows that the determinant is given by a constant times $S^2I^3$, thus we conclude. □

The consequence of Theorem 2 is that trajectories corresponding to singular controls, also called singular trajectories, may only run on the coordinate axes $S = 0$ or $I = 0$. Now, if $S\equiv 0$, then we are in the case of whole population exposed, infected, or recovered, thus the control policy is not much relevant. If $I\equiv 0$ on a time interval, then, necessarily, $E\equiv 0$, thus we are in the case of no disease. Concluding in all relevant cases, singular trajectories do not appear.

4.   Simulation of optimal control policies

In this Section, we provide simulation results obtained using the numerical algorithms described in Section 4.3. The values of parameters are reported in the Appendix A.

4.1.   Description of lockdown policies

We now look at the case where $n = 3$. Our focus is on the three possible lockdown strategies:

(1) general lockdown with home mandate for the whole population;

(2) school closure affecting school-age population;

(3) restriction of elderly interaction (for instance, by retirement homes policies, food delivery, home care, and others.).

The second strategy corresponds to control $u_2$ and an interaction matrix $l^2$ defined as the $7 \times 7$ matrix with nonzero elements in the upper left $3 {\times} 3$ square representing interactions between groups 1–3. This group represents "school interactions".

The third strategy corresponds to control $u_3$ and the matrix $l^3$ is defined as the $7\times 7$ matrix with nonzero elements in the final column and row, representing all interactions the elderly have with each other and the other age groups.

It was clear during the COVID-19 pandemic that curfews and social distancing measures specifically targeted at the elderly communities were crucial in slowing the virus spread and mortality ^[3]. Interestingly, the elderly represent the most vulnerable population while the school-age children represent the population with the strongest interaction coefficients. We assume that not all interactions are affected by the second and third lockdown measure, thus will use parameters $\tau_i$, $i = 2, 3$, to represent the amount of interactions cancelled by the lockdown. We set $l^2_{k, j} = \tau_2 l_{k, j}$ for the indices $k, j = 1, \ldots, 3$ (and zero elsewhere), and $l^3_{k, j} = \tau_3 l_{k, j}$ for the indices $k, j$ such that $k = 7$ or $j = 7$ (and zero elsewhere). Finally, because of Eq (2.2), we have for every $k$ and $j$,

Thus, $l^1$ is the matrix representing the remaining interactions not affected by the second and third strategy. We choose $\tau_2 = .66$ in our simulations. In the case of a lockdown, children will likely still see certain friends and family in their age ranges as well as their siblings. We choose $\tau_3 = .8$. This is because a lockdown of the elderly will result in continued interactions with the elderly person's family or the staff at any institution providing care.

4.2.   Description of analyzed scenarios

In Section 4.4 we study how the optimal lockdown changes depending on the delayed response to the virus. We start with immediate response and then take a look at how the optimal changes if it takes some months to identify the virus and pass legislation to begin a lockdown. We also explore how a more severe $R_0$ or a less intense maximum lockdown affect the optimal lockdown's plan.

In Section 4.5, we then study the impact that the cost per person has on the policy surrounding each of our populations being locked down. There are varying opinions about how much a lockdown at different levels costs a society, so we explore how these shifting numbers might change optimal strategies.

Lastly, in Section 4.6, we perform a cost analysis of optimal lockdown strategies for each of our three populations asking the questions "If only one type of lockdown is employed, how is the cost affected?" and "Which lockdown on its own impacts the population in the most positive way? "

All simulations are parametrized and assumed to be starting on December 20 2019. All simulation time horizons are discretized on the order of days.

4.3.   Numerical methods for optimization

To compute the optimal controls for this problem, we use CasADi, an open-source tool for nonlinear optimization and algorithmic differentiation with a large degree of flexibility ^[2]. The dynamics of (2.3) are implemented with controls using an implicit Runge-Kutta 2nd order method and the cost functional with the trapezoidal rule. All controls are initialized with a "best guess" of the maximum control over the time horizon where each time-step represents one day. We then use $\texttt{IpOpt}$ ^[46] (the interior-point optimization routine) to solve the optimization problem. The discretization of the controls is chosen to be half of the number of days, meaning over $180$ days, and the control is updated 90 times. This is more than sufficient for our purposes (no government will update policy more than every other day) and reduces run-time from a model that would choose optimal controls daily.

4.4.   Optimal lockdown varying pandemic stages, infectivity, intervention delays, and lockdown intensity

We first study the difference in costs, lockdown dynamics, and virus trajectories when allowing the lockdown to begin at different initial times. For new types of viruses, the initial time is affected by recognizing the new infection and devise screening methods.

To start the analysis, we report in Figure 1 a simulation of the age-based SEIR model (2.3) with no lockdown intervention, a replication rate, $R_0$, of 1.7 and initial S, E, I, and R values taken directly from COVID-19 data for the state of New Jersey. We take the total population of New Jersey from the 2021 census, and develop our exposed and infected as approximations of 14 days of COVID-19 data. Specifically, we begin simulations with a total population of 9.267 million people with initial infected across all age groups totaling 13736. We see that essentially all people in age groups one and three are infected over the time horizon of the virus. These dynamics mimic closely to that of a classic SIR model with a wave of cases emerging from few initial infections and then the virus dying out. Under the curtain, this model is enriched with the interactions between the various subpopulations at each time step, driven by our interaction matrix. In this case, the corresponding total cost consists only of the cost of deaths, being the lockdown controls set to zero. This cost represents an upper limit for any control policy representing non-pharmaceutical intervention and is around USD 183 billions.

In each of the following simulations, the lockdown cost is $＄70$ per person for a general lockdown, $＄40$ per person to lockdown the elderly, and $＄25$ per student to lockdown schools ^[9,12,28,38], unless otherwise stated. There is a wide range of suggested costs in the literature where a more severe lockdown costs more per person than a less severe. We consider a lockdown of only the elderly to be a first step and, thus, a less severe mitigation strategy than locking down the entire population. Therefore the cost per person should be on the lower side of the estimates. Quantifying the "cost of death" in an economic sense can be tricky. This value is assumed in the literature to be anywhere from hundreds of thousands of dollars to the assumed VSL (value per statistical life) of $＄10$ million dollars ^[24,29]. Because most COVID-19 deaths occur in the elderly, the VSL may be an inappropriately high measure of the average cost of a COVID-19 death. Therefore, while we explore the parameter in this section, we choose to follow the estimated value of $＄1.5$ million dollars to be the average cost of a COVID-19 death in the United States ^[10]. We incorporate a "cost of future deaths" to dampen end behavior effects of a lockdown fully dismissing with an influx if infected who do not have time (on the chosen time horizon) to be moved to recovered. This cost is exactly $C_{n+2} = 100C_{n+1}$. For brevity, some figures of dynamics are omitted for each simulation. For example, because the exposed mimics so closely to the infected, in future figures, the exposed compartment will be omitted.

In Figure 2, we report the best case scenario: A situation where every type of lockdown is available at a 100% level at any instant of time. Notice that this is an idealized situation for various reasons: (1) A 100% level lockdown may be not feasible; (2) usually lockdowns need to be limited in time; (3) different types of lockdowns may need coordination; and (4) change in lockdown policies cannot occur too frequently. Interestingly, if we can start our lockdown policies early, the optimal case is to have intermittent lockdowns for school activities, keeping the exponential growth of a full wave of virus under control from the beginning. This, in a sense, was adopted by many states and school districts who chose to creatively reduce school interactions by splitting students into groups and having certain groups attend on certain days, taking certain days off to do remote instruction, etc. We also see that the optimal solution reduces interactions with the elderly as much as possible. We assume the ability to lockdown 100% of each population, which may be unrealistically effective, to create a lower bound for the realized cost to the population, while we consider differing levels of lockdown in future simulations. The total cost for New Jersey is now reduced to USD 31.4 billion, with a reduction of more than $80\%$ of the cost without intervention.

In Figure 3, we see that when the virus is more aggressive ($R_0 = 2.5$), the optimal lockdown of schools and the elderly remains fairly similar while the general public must get more severe restrictions with recurring lockdowns. Even so, there is still almost an $80\%$ reduction in cost compared to a less aggressive virus with no mitigation strategy employed. This highlights the important point that a lockdown is crucial as a first mitigation; even a much more aggressive virus has a fraction of the impact of a less aggressive virus if a lockdown is employed for the former.

To further investigate real scenarios, we consider the case where a lockdown is a response to an already propagating virus. It is known that a delay in the initial lockdown most likely results in a higher number of cases ^[25]. In Figure 4, we consider a 50 days delay in interventions with an $R_0$ of 1.7, and, as expected, the lockdowns start immediately after the delay. The COVID-19 pandemic started in the United States with the first reported case in January, while the first lockdown of the state of New Jersey was not until March, so a $50$ day delay in response is a reasonable choice ^[50]. Additionally, the World Health Organization estimated COVID-19's initial $R_0$ to be between 1.4–2.4, so a choice of 1.7 is in line with a real virus ^[40]. Moreover, there is a greater strain due to the virus's ability to propagate, thus the locking down of the most vulnerable population (the elderly) and school activities are prolonged for almost the entire optimization horizon. We also see that the cost is much larger when mitigation strategies are not employed early. The total cost is USD $89.8$ billion. The main message is that delays in policy activation are much more impactful than the virus infectivity.

Additionally, for logistical, economic, or legal reasons, most populations may be unable to perform a 100% lockdown. When we reduce the severity of possible lockdown from 100% to 60%, the optimal policies correspond to lockdowns extended over time. In Figure 5, we see that even a reduced lockdown rate can greatly curb the cost of a virus to a population with a total cost of USD 47.8 billion for New Jersey. The optimal solution here is to lock down the elderly and school activities as much as possible for as long as possible with the rest of the public able to go back to no lockdown earlier. However, with less severe lockdowns, the trade-off is that there must be infrastructure in place to accommodate a longer lockdown for certain populations. In the case of school children, it is known that a prolonged loss of in-person education can be detrimental, so a short but strong lockdown is likely more beneficial ^[23]. However, it may be easier to accommodate a long-term but less severe lockdown for the elderly. A key takeaway here is that early intervention may outweigh more intense intervention in cost reduction. In fact, with the 50 day delay in Figure 4, there are 16,322 deaths, whereas with our 60% lockdown allowed but with no delay, there are only 8,298 deaths, as seen in Table 1.

Interestingly, by the end of 2020, the state of New Jersey had seen approximately $19,200$ COVID-19 deaths, with over 16,000 after the first 6 months ^[51]. This suggests that a well-tuned delay in response is important in parameterizing an epidemiological model. Earlier intervention almost assuredly would have reduced deaths in the first year of the pandemic in a significant way. This is reflected in the fact that when mimicking New Jersey parameters, the total deaths is very close to the realized COVID-19 deaths, while if we force earlier intervention in the model, the deaths are much smaller.

4.5.   Impact of cost of lockdown

We now study how varying the cost of lockdown per person, differentiating among different type of lockdowns, will affect the optimal control policy. In each simulation, our initial conditions are employed to mimic a pandemic that, if left with no intervention, is on its ascent to a peak in cases. Notice that without reinfection or more complex ways to choose parameters, a standard SIR model will generally have one wave of cases where the endemic equilibrium has no infected (all of the population has either stayed in susceptible, or moved through the system of equations to recovered) ^[30]. Figure 6 reports the optimal policy with equal daily cost of USD 10 per person for all the three types of lockdown. We notice that, despite different death rates and population sizes, the optimal lockdowns are very similar. In fact, all three run until exactly day 74, with a resurgence of only locking down the elderly. This suggests that a general lockdown would make the most sense for this situation.

We then keep the daily cost constant across lockdown types, but increase it. This causes a a differentiation in optimal lockdown policies, as seen in Figure 7. The higher daily cost, now USD 70, causes a more subtle lockdown than in Figure 6. However, while there is a reduction in lockdown for school age children and the general public, the lockdown in fact increases in time for the elderly who have a much higher death rate than the younger populations.

Now, a lockdown of the general public comes with a much higher cost than one of a subgroup of the population, in part because there are more people affected, but also because when everybody is affected, the economic strain begins to take hold on small businesses, public transportation, and more. For this reason, we see a tendency like in Figure 2 to minimize the lockdown of the public compared to our subpopulations.

In Figure 8, we see an extreme case where the parameters are the same as Figure 7, but the cost of a death is $＄10$ million dollars instead of 1.5 million. Interestingly, while the total cost has a significant increase, the optimal policy is not significantly different; a prolonged public lockdown and a slightly lessened resurgence of lockdown for the elderly. This could indicate that the optimal policy is fairly robust to a difference in cost of death.

4.6.   Cost analysis by lockdown

A major goal of this work is to exemplify the power of stratifying the population into subgroups with specific sociodemographic characteristics and using these characteristics to develop targeted mitigation strategies when addressing an epidemic. In the following simulations, we will look at a situation where it takes a month to start the first mitigation strategies, the maximum lockdown is 60%, and like in the previous section, we will use a replication rate of 1.7. Now, we study the net benefit of a stratify-able population in the face of lockdown measures through the study of the cases of no lockdown, single targeted lockdowns, and when one can use the three lockdowns independently. To set the stage, we refer back to Figure 1, which provides a baseline for how the virus would progress over the given time horizon if there were no lockdown policies enacted at all. We see immediately the recognizable peak and descent of a virus wave in the infected compartment.

Researchers using a variety of models have reported conflicting results on the impact on the viruses spread of school closures. In the case of COVID-19, many estimates show that a full school closure results in less than $5\%$ of prevented deaths ^[45]. In comparison, the costs of no lockdown outweigh the costs in Figure 9 of a school lockdown by only about 15%, meaning that optimizing school closure will at most reduce the direct societal cost by 15% (this does not take into account the loss of education long-term). It is also worth noting that even when school closure is the only method used to mitigate virus risk, the optimal strategy is not to fully lock down the school students for long-term.

Comparing the results from the school-only lockdown in Figure 9 with those of an elderly-only lockdown in Figure 10, we see an immensely more important role played by the elderly. While school closure reduces cost by 15%, a full lockdown of the elderly over the time horizon (the optimal decision) results in a 42% reduction in cost. Notice the steep reduction of slope in infected at the start of the lockdown period.

Next, we look at a situation where there is only the opportunity to lock down the entire population to varying levels with no consideration of subpopulation targeting. We see from Figure 11 that a general lockdown provides even stronger protection. In reality, locking down the public results in many unforeseen costs such as strain on small business and loss of education for children missing school. These results also suggest that, if given the opportunity, a focus should be on maximizing the possible lockdown of the elderly specifically, with this lockdown being almost the entire time horizon.

Lastly, we see from Figure 12 that when given the ability to lock down various parts of the population separately, the cost is the lowest, resulting in a reduced cost of about 70%. While this is expected to be the best-case scenario, it highlights not only the power of having a flexible population lockdown but also the sheer importance of considering heterogeneity within a population.

Holding a population to a 60% lockdown across the board may be considered unrealistic, as, one could consider the ability to close schools to be 100%, general population maybe only around 60%, and the elderly somewhere in between since they can be supported in many ways by caretakers, say 80%. In Figure 13 we see that our model is able to handle such heterogeneity in order to better fit the population being modeled and the sociodemographic variables being considered.

5.   Conclusion

In this paper, we introduce a method of combining classical epidemiological models with data about sociodemographic characteristics of a population to develop an information tool that maybe crucial to public policymakers. We choose data related specifically to age and interactions between age groups, and tune the model to census data of New Jersey. We then chose subsets of these interactions to be governed by different levels of lockdown policies, namely general interactions with the elderly, interactions most likely to occur in a school setting, and the remaining interactions after those two are considered. Costs are developed and used to guide an optimal lockdown schedule for these three strategies. We saw from preliminary COVID-19 data, that our model is strong at predicting deaths over time. We found overwhelmingly that when lockdowns must be employed, there is a heightened effect to lockdown not the most active interactions (school-age children), but the most vulnerable population (the elderly). Also, it seems that early intervention is key in curbing a virus from spreading and keeping the total mortality low. In fact, even a reduced lockdown (here 60%) greatly reduces mortality if employed early enough. Lastly, a lockdown of the general public is crucial in curbing virus spread, but if given the opportunity to fit lockdown measures to specific sociodemographic groups, the cost to the population is further lessened.

Future work will focus on an advanced model able to handle multiple sociodemographic variables as well as multiple lockdown policies. Different lockdowns may have overlapping effects, thus variables and controls must be accurately defined in the development of the model itself. Other extensions will focus on including virus mutation dynamics as well as vaccination.

All script files are written in the Python programming language and all code produced here is placed on Github and made publicly available.

Author contributions

Weightman, Piccoli, and Akinode collaborated across various stages of the paper. In the conceptualization phase, Weightman and Piccoli contributed their expertise. Methodologically, Weightman, Akinode, and Piccoli worked together to develop research methodologies. During writing, Weightman and Piccoli drafted the original manuscript, while also engaging in review and editing alongside Akinode. Piccoli provided supervision throughout the project, and played a crucial role in securing funding.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors acknowledge the support of the Joseph and Loretta Lopez Chair endowment, the support of the NSF CMMI project # 2033580 "Managing pandemic by managing mobility" in collaboration with Cornell University and Vanderbilt University, Rutgers Global grant "Innovative models for pandemics with mutating viruses: mitigation strategies for vulnerable populations", and the Institute for Advanced Study in Princeton for hosting last author.

Conflict of interest

Benedetto Piccoli is an editor-in-chief for [Networks and Heterogeneous Media] and was not involved in the editorial review or the decision to publish this article. All authors declare that there are no competing interests.

Appendix A: Data, Parameters, and Contact Matrix

Each simulation reported in Section 4 uses the SEIR parameter set in Table A1 based on data collected during the Covid-19 pandemic.

To calibrate our model, we choose an age structure to drive the interactions and a lockdown as our mitigation strategy. We begin by partitioning our population into the seven age groups in Table A2.

We partition the total population into these age groups based on percentages provided by the United States census ^[9]. For each piece of information captured in a SEIR model, initial conditions and parameters must be estimated. The Johns Hopkins University (JHU) repository of Covid-19 data collected case and death numbers for each state in the US carefully from the start of the pandemic until early 2023 ^[37]. We then choose a starting date and use the JHU data to get the number of infected. From this, we build initial conditions for our compartments.

The number of exposed is set equal to a certain percent of the infected from the first 5 days of the model (December 16–20 2019). Then, the number of infected is taken to be the number of cases from the first 5 days of the model minus exposure assigned proportionally based on the number of people in each age group. We consider the number of recovered to be the sum of cases reported from the start of reporting until two weeks before the chosen start date. Lastly, we define the number of susceptibles to be the sum of the population in each group minus the sum of the other compartments.

We take the general SEIR model and couple it with an age-based interaction matrix. Our interaction matrix ${l}$ is developed starting from data from ^[41], which takes data from population-based contact diaries in eight countries and projects to various other countries using a Bayesian model. The contact matrix is actually inferred data for the United States. The age groups in the original contact matrix provided by the authors contains $16$ age groups at a higher resolution than we include here, so averages are taken where two of the original age groups make up our age group. For example, the second element of ${l}$ is $0.8003$ which makes up two of the age groups in the original dataset; the 5–9 element being $1.1013$ and the 9–14 element being $0.4994$, the average of which is $0.8003$. Notice also that this matrix is not symmetric. The data is used to develop vectors for each age group which represent their interactions at home, interactions in the workplace, interactions at school, etc. Therefore, a phenomenon occurs where, for example, the $70+$ age group, in their daily interactions, are very likely to interact with someone of the 40–59 age group, as these are often caretakers, children, or employees at consumer establishments. So, their element of interacting with 40–59 is high. However, relative to other interactions, 40–59 age group members are unlikely to interact with the $70+$ age group the most, resulting in their relative interactions with $70+$ being much lower. It is worth noting that while we used only data from a more recent source, the contact matrices found in other papers are very similar in structure to the ones we use ^[8].

We then use matrix l coupled with general information about the population of the United States of America to separate our subpopulations $l^i$.

Death rates are computed from CDC Covid-19 data, but follow the general pattern of many viruses where as age increases, as does death rate ^[39]. They are as follows in order of our age groups: 0.00016, 0.00016, 0.00006, 0.00007, 0.00229, 0.01915, 0.13527.

In the end, we would like to use optimal control to answer the question of which sub-groups lockdown would be most impactful in curbing the effect a virus as on a population. As such, the cost function must be developed carefully to reflect the population in question and the age categories being locked down. To achieve this, we define the cost per child per day of a school lockdown by averaging estimates provided by the Survey of Income and Program Participation and the Current Population Survey ^[9,28]. This number fell in the range of as low as $＄25$ per day to as high as $＄60$ per day. It is important to note that this does not include the extremely high potential reduction of education. For the cost of locking down the elderly populations, we consider the same cost as locking down another adult person, which we find to be ＄40–＄350. This number has a large range, with some reporting a cost of as much as $＄350$ per day for a strict lockdown and much less for a less severe, more flexible lockdown ^[12,38].

Lastly, choosing parameters such as costs for a cost function requires quantifying the "cost of death" in an economic sense. We follow the estimates of one analyst who claims $＄1.5$ million dollars to be the average cost of a Covid-19 death in the United States ^[10]. All codes used to generate simulations are available at ^[48].