Age-structured non-pharmaceutical interventions for optimal control of COVID-19 epidemic

In an epidemic, individuals can widely differ in the way they spread the infection depending on their age or on the number of days they have been infected for. In the absence of pharmaceutical interventions such as a vaccine or treatment, non-pharmaceutical interventions (e.g. physical or social distancing) are essential to mitigate the pandemic. We develop an original approach to identify the optimal age-stratified control strategy to implement as a function of the time since the onset of the epidemic. This is based on a model with a double continuous structure in terms of host age and time since infection. By applying optimal control theory to this model, we identify a solution that minimizes deaths and costs associated with the implementation of the control strategy itself. We also implement this strategy for three countries with contrasted age distributions (Burkina-Faso, France, and Vietnam). Overall, the optimal strategy varies throughout the epidemic, with a more intense control early on, and depending on host age, with a stronger control for the older population, except in the scenario where the cost associated with the control is low. In the latter scenario, we find strong differences across countries because the control extends to the younger population for France and Vietnam 2 to 3 months after the onset of the epidemic, but not for Burkina Faso. Finally, we show that the optimal control strategy strongly outperforms a constant uniform control exerted over the whole population or over its younger fraction. This improved understanding of the effect of age-based control interventions opens new perspectives for the field, especially for age-based contact tracing.

Introduction Following its emergence in December 2019, COVID-19 has become an international public health emergency [1]. The infection has many similarities with that caused by influenza virus regarding its clinical manifestations and transmission mechanisms [1]. However, contrary to seasonal influenza, COVID-19 has become pandemic by spreading rapidly among completely naive host populations, i.e. with no pre-existing immunity [2][3][4][5]. Initially, no pharmaceutical interventions such as vaccines or treatments were available and it usually takes several months before their deployment. For this reason, non-pharmaceutical intervention (NPI) strategies, such as social distancing, are key to controlling the pandemic [6].
When an intervention is summarized by one or few parameter values, identifying an optimal strategy according to some criterion variable can readily be done, e.g. using a gradient approach [7]. Things become more challenging when the intervention parameter value is a function of time. Optimal control theory [8], specifically addresses this issue by identifying a function of time such that some criterion is optimized. This has allowed studies to identify optimal non-pharmaceutical interventions to control infectious diseases such as influenza and COVID-19 [9][10][11][12]. However, a strong limitation of these studies is that they all ignore at least one aspect of the host population structure. First, infection parameters (contagiosity, recovery) vary with infection age, i.e. depending on the number of days since infection. Second, hosts vary in age. The latter point is particularly important because, in addition to being a function of time since the onset of the outbreak, optimal strategies involving physical distancing can also vary depending on host age [13][14][15][16][17]. Accounting for two dimensions, time and host age, make the optimization procedure much more challenging because optimal control theory is usually applied to ordinary differential equations (ODEs) -something very common-while here we are working on partial differential equations (PDEs) -which is less common, and much more challenging. Here, we address this challenge and identify interventions varying in intensity with time and host age that significantly reduce morbidity associated with COVID-19 at a minimal cost. Furthermore, we compare the situation in countries with contrasted agestructure, namely Burkina-Faso, France, and Vietnam, to show how this affects optimal strategies.
The age structure of the population is a known key determinant of acute respiratory diseases, especially when it comes to infection severity. For example, children are considered to be responsible for most of the transmission of influenza virus [18] but the related hospitalization and mortality burden is largely carried by people of ages over 65 years [19,20]. While much remains unknown about the COVID-19 epidemics, evidence to date suggests that mortality among people who have been tested positive for the coronavirus is substantially higher at older ages and near zero for young children [3,21]. Moreover, the infectiousness of an individual has been reported to vary as a function of time since infection [22], which is known to affect epidemic spread [23][24][25][26].
Our model for the disease stage-progression is structured both by a (continuous) age of the host and a (continuous) age of infection. A variety of epidemiological models allow for one or the other type of structure [27][28][29][30], starting with a seminal article from the 1920s [23]. However, models allowing for a double continuous structure are rare [30][31][32][33][34][35][36][37], even though it is particularly suited to investigate infections such as COVID-19, with strong effects of host and infection age. Indeed, in addition to taking into account the age structure of the host population, as well as the gradient of disease severity from mild to critical symptoms, the model readily captures the variation in infectiousness as a function of the time since infection. From a theoretical point of view, age-structured models have been proposed to investigate the spread of acute respiratory diseases [38][39][40][41][42], and some rare models of acute respiratory diseases consider both structures as continuous variables [30,32], although not in the context of optimal control theory.
We first introduce the mathematical model and define its parameters and outputs. Next, we characterize the optimal control strategy that minimizes the number of deaths as well as the cost due to the implementation of the control strategy itself. The main body of the results then follows. We first analyze the epidemic spread without any intervention for three countries with contrasted age distributions (Burkina-Faso, France, and Vietnam). Second, we compare the performance of optimal control in terms of deaths and hospitalizations for different costs of the control measure. Finally, we compare our optimal control strategy to two other strategies that use the same amount of resources to control the outbreak. Finally, we discuss our model in results in the context of the current pandemic and identify perspective for future work.

Materials and methods
An age-structured epidemiological model Model overview. We denote by S(t, a) the density of individuals of age a 2 [0, a max ] that are susceptible to the infection at time t 2 [0, T]. These individuals can become infected with a rate called the force of infection and denoted λ(t, a). We assume that a fraction p of these individuals are paucisymptomatic, which means that they will develop very mild to no symptoms, and enter group I p . Note that p is likely to depend on age, but, because it is currently unknown, we assume it to be constant. This class I p can also be interpreted as the fraction of the population that will not isolate themselves during their infection. Other individuals are assumed to develop more symptomatic infections, either severe I s with proportion q(a) depending on the age a, or mild I m with proportion 1 − q(a).
Each of the three infected host populations are structured in time since infection, so that I v (t, a, i), v 2 {p, s, m}, denotes the density at time t of individuals of age a that have been infected for a duration i 2 R þ . Upon infection, all exposed individuals are assumed to remain non-infectious during an average period i lat . Next, they enter an asymptomatic period during which they are infectious. Only I m and I s develop significant symptoms after an average time since infection i sympt , which can allow them to self-isolate to limit transmission. During their infection, individuals can recover at a rate h v (a, i) (v 2 {p, m, s}) that depends on the severity of the infection and the time since infection i. Severely infected individuals of age a may also die from the infection at rate γ(a, i).
The infection life cycle is shown in Fig 1.
Remark. Contrarily to classical SEAIR models, disease-stage progression in our model is not captured by discrete compartments (exposed, asymptomatic, and infected) with exponentially distributed waiting times to switch between compartments. The advantage of our formalism is that disease progression can be modeled using a continuous variable, called the time since infection (in days) denoted here by i. Every infected person then remains in the "infection compartment" from exposure until recovery (or death). Latency from exposed to asymptomatic and time of symptoms onset are not needed for this modeling approach because these are captured through the functions describing the transmission rate, the mortality rate, and the recovery rate at time i post infection. More precisely, the average latency from exposed to asymptomatic (i lat ) is simply mentioned to define the average time to infectiousness onset (i sympt ), and also to help the readers to understand the model flow diagram (Fig 1). On the other hand, i sympt is used to define infectiousness reduction factors (ξ s , ξ m ) and the mortality rate due to the infection (γ).
Age-structured transmission and severity. We use two components to model the infection process. First, we define the transmission probability β v (a, i) (v 2 {p, m, s}) for each contact between an infected of age a and a susceptible person, which depends on the time since infection i. Second, we introduce the kernel K(a, a 0 ) that represents the average number of contacts by unit of time between an individual of age a 0 and an individual of age a. The numerical values of this contact matrix are based on data from an earlier study [43]. The force of infection underwent by susceptible individuals of age a at time t is then given by lðt; a; cÞ ¼ ð1 À cðt; aÞÞ where c = c(t, a) is the percentage of contacts reduction towards people with age a, due to control measures, at time t. The total force of infection at time t in the whole population is computed as R a max 0 lðt; a; cÞda. The dynamics of newly infected individuals (i.e. i = 0) in each group is thus defined by I s ðt; a; 0Þ ¼ ð1 À pÞqðaÞlðt; a; cÞSðt; aÞ; I m ðt; a; 0Þ ¼ ð1 À pÞð1 À qðaÞÞlðt; a; cÞSðt; aÞ; I p ðt; a; 0Þ ¼ plðt; a; cÞSðt; aÞ: Susceptible hosts of age a at time t (S(t, a)) are exposed to the virus with a force of infection λ(t, a). A fraction p of exposed individuals, which are infected since time i, will never develop symptoms and enter the group of paucisymptomatic infections (I p (t, a, i)). The rest will develop symptomatic infections, either severe (I s (t, a, i)) with proportion q(a) depending on age a of individuals, or mild (I m (t, a, i)). Exposed individuals remain non-infectious for a duration i lat after infection. Next, they become asymptomatic infectious and only symptomatic infected will develop symptoms at time i sympt after infection. Infected individuals recover at rate h v (a, i). Only severely infected of age a die from the infection at rate γ(a, i). Notations are shown in Table 1. where μ nat denotes the natural mortality rate when hospitals are not saturated. We assume that this rate increases significantly as soon as the number of severe cases exceeds the healthcare capacity H sat and denote by μ add this additional death rate due to hospital saturation. We apply the same reasoning by assuming that the disease-related mortality can increase because of hospital saturation. Therefore, severely infected individuals of age a who have been infected since time i die at time t at rate γ(a, i, H(t)) defined by gða; i; HðtÞÞ ¼ ðg dir ðaÞ þ g indir ða; HðtÞÞÞ1 ½i sympt ;i s max � ðiÞ: where γ dir and γ indir are mortality rates directly and indirectly due to the COVID-19 respectively. The disease-related mortality occurs after the emergence of symptoms and before the mean final time of infection for severe cases, i.e. for i 2 ½i sympt ; i s max �: Finally, infected individuals of age a infected since time i recover at rates h s (a, i), h m (a, i), and h p (a, i) for severe, mild, and paucisymptomatic infections respectively.
The boundary conditions (3)  for each ða; iÞ 2 ½0; a max � � R þ . Numerical values for initial conditions are detailed later. Using equation system (3) and an integration over i of system (5), one may observe that the total population N defined by Eq (1) is strictly decreasing since it satisfies the following inequality: @N @t ðt; aÞ � À m nat ðaÞNðt; aÞ; 8a 2 ½0; a max �; 8t � 0: This is due to the fact that population aging and births are neglected in this model since we consider a time horizon of only one year. Basic properties of the model such as existence and positiveness of solutions is out of the primary scope of our study but these can be specifically addressed using an integrated semigroup approach and Volterra integral formulation (see e.g. [44][45][46][47] and references therein). More specifically, one can refer to [31], where the well-posedness of an epidemiological model with a double continuous structure is handled.

Epidemiological outputs, model parameters and initial conditions
In this section we briefly describe some useful epidemiological outputs, the shape of agedependent parameters considered for the simulations of model (3)-(5), and the initial conditions. All state variables and other parameters are summarized in Table 1.
Epidemiological outputs. In addition to the total number of hospitalized cases H(t) at time t defined by Eq (4), we define additional epidemiological outputs. The first one is the number of non-hospitalized cases (N H (t)) which encompasses paucisymptomatic, mildly infected, and severely infected but not yet hospitalized hosts.
For the cumulative number of deaths, we distinguish between those directly due to COVID-19 infections (D cum dir ðtÞ), and those indirectly due to the epidemic (D cum indir ðtÞ), which originate from the saturation of the health system: where D dir ðtÞ and D indir ðtÞ are the number of deaths at time t respectively defined by HðtÞÞ Every aforementioned output implicitly depends on parameter c = c(t, a), which we will omit in the notations when no confusion is possible. However, for clarity, we do explicitly write this dependence to compare public health measures. The relative performance between two strategies c 1 and c 2 , denoted by Δ(c 1 , c 2 ), is estimated by assessing the cumulative number of deaths in the whole population during the T days of control period with the strategy c 1 relatively to deaths with the strategy c 2 . Formally, we have Hence, a relative performance Δ(c 1 , c 2 ) = 0.1 implies that strategy c 1 reduces the number of deaths by 10% relatively to strategy c 2 .
Model parameters. Mortality rates. We assume that indirect mortality, i.e. not directly due to COVID-19, increases when the number of hospitalisations H(t), at time t, exceeds a healthcare capacity threshold H sat (which we approximate with the maximal intensive care capacity). The natural mortality rate then increases by μ add (a, H) for the whole population, and by γ indir (a, H) for severely infected individuals of age a. These rates are modeled by logistic functions that are arbitrarily chosen as: for each a 2 [0, a max ], meaning that the natural mortality rate is only increased by 1%, while the disease-induced mortality rate γ is doubled. Indeed, according to [48], less than 50% of patients in critical care will die in case of no saturation of hospitals. Transmission rates. The infectiousness of an individual aged a who is infected since time i, is given by β v (a, i) (v 2 {s, m, p}). Based on estimates described in [22], we assume that β v does not depend on host age a, i.e., β v (a, i) = β v (i). As discussed later, this assumption is only made for parameterization purpose and does not impact the general formulation of the model proposed here.
The transmission rate at a given day i post infection of a given type of infectious host is defined such that bðiÞ, for v 2 {s, m, p}. As detailed below, α is a scaling parameter obtained from the value of the basic reproduction number R 0 , which is the mean number of secondary infections caused by an infected host [24]. As [22], we assume that parameter � b, which strongly depends on the generation interval, follows a Weibull distribution � b � Wð3; 5:65Þ. Finally, parameters ξ v (i) are factors that capture variations in infectiousness based on the type of host. For paucisymptomatic individuals, for instance, these are assumed to be constant (ξ p (i) = ξ p ), while the reduction factor in more symptomatic infections (severe and mild) is assumed to vary after symptom onset to capture admission in a healthcare facility or self-isolation at home. More precisely, we assume that These two functions are chosen arbitrarily by assuming that individuals do not isolate before symptoms onset (i � i sympt ), and that isolation is stronger when symptoms are more severe (Fig 2A). We therefore assume that the transmission probability � b is divided by 10 (respectively 2) every day after the average time of symptoms onset for individuals severely (resp. mildly) infected.
Recovery rates. We assume that recovery rates h v (a, i), v 2 {s, m, p}, of infected individuals of age a infected since time i are independent of a and take the following form: That is, one can recover from severe (resp. mild and paucisymptomatic) infections only after a time since infection i s max (resp. i m max ) corresponding to the mean duration of infection. Initial conditions. The initial susceptible population S 0 and epidemic size I 0 are given in Table 1. Since, initially, screening is usually restricted to individuals with severe symptoms, we assume that all initial cases are severe infections. Thus, we set the initial severely infected individuals, which we assume to be uniformly distributed with respect to the time since infection i on the interval ½0; i s max �. Using estimates from [55,60] on the age distribution of hospitalised people, we derive an estimation of I s,0 (a, i) for each ða; iÞ 2 ½0; a max � � R þ . Next, following the life cycle (Fig 1), we obtain an estimation of the total initial infected population by I s;0 ða;iÞ ð1À pÞqðaÞ . From there, we deduce the initial mildly and paucisymptomatic infected populations, which can be denoted respectively by I m;0 ða; iÞ ¼ 1 À qðaÞ qðaÞ I s;0 ða; iÞ and I A;0 ða; iÞ ¼ p qðaÞð1 À pÞ I s;0 ða; iÞ:

Optimal intervention
As explained above, our goal is to find an optimal control strategy that is allowed to vary depending on the number of days since the onset of the epidemic (t) and on host age (a). In this section, following well established methodology in optimal control theory [13][14][15][16]61], we search for the optimal control effort function c � that minimizes the objective functional J : L 1 ðR þ � ½0; a max �Þ∍c7 !JðcÞ 2 R, where wherein the set U is defined by with c max � 1 a positive constant. That is to say, the function c � will minimize the cumulative number of deaths during T days, as long as the cost of the control strategy is not too large. Let (S, I s , I m , I p , R) be a given solution of equation systems (3)-(5) then let λ and H be respectively defined by Eqs (2) and (4). After some computations (S1 Text), we find that the adjoint system of (5) Finally, the Hamiltonian H of (11) is detailed in S1 Text. By solving @H @c ¼ 0, it comes that c � ðt; aÞ ¼ maxð0; minðĉðt; aÞ; 1ÞÞ; with λ 0 detailed in S1 Text. We also assume that the cost B(a) of the control measure over individuals aged a 2 [0, a max ] is proportional to their density in the initial susceptible population S 0 , i.e.
where B � 2 R þ is a variable parameter characterizing the relative cost in implementing the strategy.
Additionally, one could also factor in the age distribution of the economic cost on the shape of the function B. For example, the economic cost is likely to be more important for the working population (i.e. age group 20-60) compared to the older, mostly retired, population. However, in absence of relevant references regarding this topic, we stand with our primary assumption.
The state system (3)-(5) and the adjoint system (12) and (13) together with the control characterization (14) form the optimality system to be solved numerically. Since the state equations have initial conditions and the adjoint equations have final time conditions, we cannot solve the optimality system directly by only sweeping forward in time. Here, we use an iterative algorithm, forward-backward sweep method [8]. In other words, finding c � numerically first involves solving the state variables (3)-(5) forward in time, then solving the adjoint variables (12) and (13) backward in time, before finally plugging the solutions for the relevant state and adjoint variables into (14), subject to bounds on the control function. The proof of the existence of such control is standard and is mostly based on Ekeland's variational principle [62]. Therefore, we assume the existence of the solution to the above problem and refer to [13] for additional details.

Results
Here we consider three countries as case studies: Burkina Faso, France, and Vietnam. Their populations have quite contrasted age-structure and social or physical contacts (Fig 3). In Burkina Faso, the vast majority (96.1%) of the population is less than 60 years old, whereas this is less the case in Vietnam and in France (87.7% and 73.4% respectively) as shown in Fig 3A and  3B. It shows that a higher proportion of the population is older than 60, hence at risk for COVID-19 infection, in France 26.6%, in comparison to Vietnam 12.3%, or to Burkina Faso 3.9%. Also, contacts are more frequent among the older population in France compared to Vietnam (Fig 3D and 3E). By contrast, very few contacts are observed among older populations in Burkina Faso (Fig 3C).

Global sensitivity analysis
We study the sensitivity of infected individuals, hospitalizations, and deaths to five parameters: the proportion of paucisymptomatic infections (p), the average time of symptoms onset (i sympt ), the infectiousness reduction of paucisymptomatic infections (ξ p ), the healthcare capacity (H sat ), and the population structure (including the natural mortality, the size of the population, age-structure, and social contacts). The variation range of the above parameters is  Table 1. Sensitivity indices are estimated by fitting an ANOVA (Analysis Of Variance) linear model, including third-order interactions, to the data generated by simulation. This ANOVA linear model fitted well with 99% of variance explained.
Overall, the population structure is the main parameter highlighted by the sensitivity analysis. It explains 70% of the variance for the number infected individuals and 40% for hospitalizations and deaths (S1 Fig). The population structure is followed by ξ p , p, and i sympt which have a similar impact on the number infected individuals with a slight dominance of ξ p (S1 Fig). For hospitalizations and deaths, the population structure is followed by p with 40% and 30% of the variance explained respectively; while ξ p and i sympt have very marginal impact (S1 Fig). Finally, the importance of H sat on the three output variables is largely negligible, with, however, greater importance on deaths as compared to hospitalizations and infections (S1 Fig).

The basic reproduction number R 0
An explicit expression of the R 0 of the model defined by Eqs (3)-(5) is difficult to obtain in general. We show in S2 Text that it is possible to write R 0 ¼ a � rð � U Þ, where α is the scaling parameter introduced earlier, and rð � U Þ is the spectral radius of the next generation operator � U defined on L 1 (0, a max ) by where S 0 is the initial susceptible population, K is the contact matrix, and ω(a, i) is the infectiousness of individuals of age a infected since time i (S2 Text). It follows that In our numerical approach, we set R 0 = 3.3 [63,64] for all three countries and corresponding values for S 0 and K for each country. We then successively determine rð � U Þ and α by (15) and (16) respectively.

Uncontrolled epidemic
We first use the model (3)-(5) to describe the outbreak of the epidemics for all three countries, in absence of public health measure (i.e. c � 0), with R 0 = 3.3 and other parameters defined previously and summarized in Table 1.
The peak of the epidemics is then reached approximately at day t = 51 for hospitalised people, and day t = 46 for non-hospitalised people in the France scenario (Fig 4E). Such times to peaks for hospitalised and non-hospitalised people are 47 and 41 (resp. 50 and 45) for Burkina Faso (resp. Vietnam) scenario (Fig 4A, resp. Fig 4I). The delay between the two peaks is explained by the latency time i sympt for symptoms onset (Table 1).
In absence of control measures, the healthcare capacity is quickly exceeded, about twenty days for the 'France' scenario (Fig 4E), and the number of deaths increases sharply from then on. A similar configuration is observed for Vietnam (Fig 4I). By contrast, because of the very low healthcare capacity in Burkina Faso, the health system is exceeded within a few days ( Fig  4A). However, this overloading does not have the same consequences in terms of mortality in Burkina Faso compared to France and Vietnam. This is partially explained on the one hand by the fact that less than 4% of the population is above 60 years in Burkina Faso (Fig 3A) and on the other hand by the fact that very few contacts are observed within the older population in Burkina Faso compared to France or Vietnam (Fig 3B and 3D).

PLOS COMPUTATIONAL BIOLOGY
At the end of the simulation (t = 150 days), without any control measures, the herd immunity threshold (1 − 1/R 0 � 69.7%) is clearly reached in France and in Vietnam (Fig 5), where the average size of the epidemic (severe, mild, and paucisymptomatic infections) is close to 90%. The threshold is also reached in Burkina Faso but only with an epidemic size of 78% ( Fig  5). In all three countries, the proportion of the population less than 20 years old that has been infected is around 93%. In the  group, we find a similar percentage in France and Vietnam (94%), but only 65% in Burkina Faso. This proportion then decreases for the population older than 60, more or less quickly depending on the country: it is around 73% in France, 56% in Vietnam, and 33% in Burkina Faso. Finally, among the infected population, more than 98% are less than 60 in Burkina Faso, while this proportion is 92% in Vietnam and 76% in France. This age structure of infected populations is particularly important since most of the infections that occur in the young population do not require hospitalization (Fig 4C, 4G and 4K) while people older than 60 represent the age class with the highest cumulative number of deaths ( Fig  4B, 4F and 4J).

Optimal intervention
We now investigate the result of implementing an optimal intervention that accounts for the age structure of the population. Strategies performances are here compared in terms of the cumulative number of deaths for three costs of control measures; low (B � = 10 2 ), intermediate (B � = 10 3 ), and high (B � = 10 4 ).
The optimal control strategy varies in time and depends on host age. In general, regardless of the country (Burkina Faso, France, or Vietnam), the control is stronger early in the epidemic and for older populations (Fig 6, S2 and S3 Figs). Overall, the level of optimal control is lower in Burkina Faso compared to France and Vietnam (Fig 6, S2 and S3 Figs). If the cost of implementing the measures B � is intermediate or high, the optimal control is almost restricted to individuals above 55 and to the first third of the time interval considered, with a significant reduction in deaths (Fig 6D and 6E, S2 and S3 Figs). In France, the relative performance of the optimal control c � compared to a 'doing nothing' scenario (Δ(c � , 0)) is at least 92% (resp. 82%) when the cost is B � = 10 3 (resp. 10 4 ). For Burkina Faso, Δ(c � , 0) is at least 50% (resp. 4%) when B � = 10 3 (resp. 10 4 ). Finally, for Vietnam Δ(c � , 0) is at least 87% (resp. 62%) when B � = 10 3 (resp. 10 4 ). In the case of Burkina Faso, note that the level of the optimal control is quite low when the cost of implementation is high (S2 Fig), and, as a result, the effect of this control in reducing mortality at the population level is negligible. This is due to the relatively small number of deaths in the whole population in Burkina Faso without any control measures (Fig 4A).
If the implementation of the control measure comes at a low cost (B � = 10 2 ), the optimal control significantly extends to younger populations in all three countries (Fig 6A, S2 and S3  Figs), with a maximum intensity reached near the 4th month of the epidemics and a steady decrease until the end of the control period. Overall, the optimal control lasts less longer in Burkina Faso (S2 Fig) compared to the cases of France and Vietnam (Fig 6A and S3 Fig). At first, the control is mainly applied to people above 35 in all three countries (Fig 6A, S2 and S3  Figs). But, while the control extends to people less than 35 in France and Vietnam after 2 or 3 months (Fig 6A and S3 Fig), such an extension is very moderate (or even negligible) in Burkina Faso (S2 Fig). The resulting reduction in the number of deaths is very pronounced with a relative performance Δ(c � , 0) of at least 80% (resp. 99%, 97%) in Burkina Faso (resp. France, Vietnam).

Performance and practical implementation
To illustrate how the strategy identified using optimal control theory outperforms constant uniform control exerted over the whole population or its younger fraction, we derive optimal strategies that do not vary in time and use the same amount of 'resources' (that is the same cumulative cost). Assuming a relatively high cost B � = 10 3 , we first investigate a control strategy that targets the younger fraction of the population (Fig 7A), a second strategy that uniformly targets the whole population (Fig 7B). Both strategies have a control level c max = 0.95 and vary in duration (the total amount of resource used being constant).
In France, when targeting the population uniformly, the epidemic is under control for approximately 60 days. However, once the control resources are exhausted, the epidemic reemerges ( Fig 7C). With the (longer) control over the younger fraction of the population, the first epidemic peak is slightly delayed and the epidemic appears to be under control for a

PLOS COMPUTATIONAL BIOLOGY
Age-structured non-pharmaceutical interventions for optimal control of COVID-19 epidemic longer time (180 days). Unfortunately, resources also become exhausted and a second peak appears a few months later (Fig 7C). In both cases, i.e. targeting the whole population or its younger fraction only, the results are similar. A uniform control leads to a cumulative mortality over the time period of interest comparable to that without any control measure (Fig 7D). The performance of the optimal control relative to a uniform control is approximately 92%. In the end, only 55% of the whole population has been infected with the optimal control but at least 85% with a uniform control (Fig 7E). Similar results are obtained for the case of Vietnam (S4 Fig). By contrast, for the case of Burkina Faso, regardless of the control strategy (optimal, uniform without targeting, and uniform targeting the younger fraction) the proportion of the population that is infected is approximately the same as without control (78%). The herd immunity threshold (1 − 1/R 0 � 69.7%) is then reached for all the three control measures ( Fig  8E) and the epidemic cannot restart (Fig 8C). The cumulative mortality with a uniform control exerted over the whole population or over its younger fraction is comparable to that without any control measure (Fig 8D). However, the optimal control performs at least 50% better than PLOS COMPUTATIONAL BIOLOGY the uniform control (whether it concerns the whole population or its younger fraction). This result holds despite the nearly identical proportions of infections in the population under 60 years of age (Fig 8E).
A practical issue regarding the implementation of the optimal control strategy we identified is that it is a continuous function. One possibility to address this problem is to use step functions. In Supplementary S5 Fig, we subdivided the population into 10-year amplitude classes and updated the control every 3-weeks. Importantly, even though it is assumed to be constant during each 3-weeks period for each age-class, the control intensity directly originates from the results of the continuous optimal control strategy. This discrete implementation of the optimal strategy achieves similar efficiencies (S5 Fig), with a relative performance of 91% compared to a doing nothing scenario.

Discussion
Non-pharmaceutical public health interventions can be implemented either to mitigate the COVID-19 epidemic wave or to suppress the wave long enough to develop and implement a vaccine or a treatment. Here, we explicitly factor in the age heterogeneity of the host population in the identification of the optimal allocation of the control efforts. We focus on three countries (Burkina Faso, France, Vietnam) with contrasted population age-structures and social or physical contacts (Fig 3).
We use optimal control theory [8] to characterize an optimal strategy that significantly reduces the number of deaths while being sustainable at the population level. Our formulation assumes a quadratic cost for the control effort. Overall, the optimal control lasts less in Burkina Faso compared to France and Vietnam. With this strategy, we find that the intensity of the control is always relatively high on the older fraction of the population during at least a hundred days, before decreasing more or less rapidly depending on the cost associated with the control and the social structure of the host population. The control over the younger fraction of the population is weak and only occurs when the cost associated with the optimal control is relatively low. However, while the control applies to the younger population in France and Vietnam after 2 or 3 months, this is very moderate (or even negligible) in Burkina Faso. This late control over the younger part of the population mimics the results found in [10], which did not include host age but found that control did not peak right away. Intuitively, if control strategies come at a high cost for the population, it is best to focus on the age classes that are the most at risk. Conversely, if the control measures are more acceptable to the population, the optimal strategy is to aim wide to completely suppress the epidemic wave.
Information on the natural history of paucisymptomatic infections of COVID-19 remains relatively poorly documented [65,66]. We estimated that a proportion p of infected individuals remain asymptomatic throughout their infection, but this proportion remains largely unspecified in the literature [65,66]. We performed a sensitivity analysis of p on the optimal control strategy c � . Overall, the proportion of paucisymptomatic infections has marginal effects on the optimal control strategy (S6 Fig). The optimal control remains strong over the older population from the beginning of the epidemic, before being progressively relaxed. The control over the younger population is weaker and occurs only if the control cost itself is low. But, the level of control over the younger fraction increases when the proportion of paucisymptomatic infections decreases. Further, for high values of B � (10 3 or 10 4 ), the shape of the optimal control is qualitatively the same when the proportion p varies, except for extremely high values of p = 0.9 and B � = 10 4 , for which the control becomes low for the whole population (S6 Fig). The interpretation is that if the epidemic cannot be stopped, the best strategy is then to reduce mortality by protecting the population the most at risk (here the older population). However, with a low value of B � (10 2 ), different shapes of the optimal control yield the same result since there are enough resources to stop the epidemics.
Given the leverage represented by school and university closures, we investigated the effect of targeted control measures over individuals aged under 25. Our results show that NPIs targeting the younger fraction of the population are not very efficient in reducing cumulative mortality unless they can be implemented strongly and for a relatively long period. Indeed, the number of deaths recorded when implementing a control only over the younger population is similar to a doing nothing scenario for cases of Burkina Faso, France, and Vietnam (Figs 7 and  8, and S4 Fig). This result seems independent from the age-structure and social contacts of the population considered. However, as we further discuss below, the variation of the transmission probability with age could matter.
The formulation of the objective functional considered here aims to minimize the cumulative number of deaths. However, other objective functions could be envisaged and factor in, for example, long-term hospitalizations or long-term health consequences. Practically, this could be implemented similarly to the one presented here by considering the number of hospitalizations as the variable to be minimized and costs associated with long-term hospitalizations and long-term health consequences. This formulation may indeed be interesting to investigate in detail but would deserve a dedicated study.
Our model is an extension of models based on ordinary differential equations that tackled the issue of the optimal control of COVID-19 outbreaks [9][10][11][12]. However, the whole population is here structured by age (a) and additionally by the time since infection (i) for infectious individuals, which echoes the model developed in [64] using a discrete-time formulation of the infection. With our continuous structure, we show that the number of new cases I N (t, a) at time t in individuals of age a is given by the renewal equation Kða; a 0 Þoða 0 ; iÞI N ðt À i; a 0 Þda 0 di; where K is the contact matrix and ω(a, i) is the infectiousness of individuals aged a which are infected since time i (S2 Text). For parameterization purposes, we assume that ω(a, i) is the product between the proportion of individuals of age a in the whole population and the infectiousness � bðiÞ of individuals infected since time i. This is potentially a limitation -not in the model formulation proposed here, but rather with regards to our parameterization using the existing literature-since infectiousness � b could depend on the age a, thereby creating an additional heterogeneity to that of the time since infection i. This issue can be particularly important since some studies suggest a low risk of transmission in the young population (e.g. [67]). On the other hand, although superspreading events (of young people) have been documented, there is still much uncertainty about their relative role in the spread of the epidemic and about their origin (superspreading could be linked to environmental conditions, such as massive gatherings, rather than individual properties). Therefore, assuming independence from age appears to be the most parsimonious assumption given the current data.
Another potential limitation is the lack of gender structure and comorbidities in the model formulation. Given the observed male-biased in mortality during the COVID-19 pandemic, it has been suggested that males are more at risk of developing severe infections [68]. This heterogeneity could readily be introduced in the model if its biological importance is further demonstrated.
Contact networks have an important role in transmission dynamic models [69]. Epidemic models that determine which interventions can best control an outbreak may benefit from accounting for social structure and mixing patterns. Contacts are known to be assortative with age across a given country, but regional differences in age-specific contact patterns are noticeable [43]. The current model could be modified to explore epidemiological dynamics in a spatially structured population with non-homogeneous mixing, e.g. by using a meta-population model [70].
Another potential extension of the model would be to allow for the isolation of symptomatic cases and their contacts, following the method developed in [71] and applied recently to digital contact tracing [22]. Indeed, these measures strongly depend on the relative timing of infectiousness and the appearance of symptoms, which are both well-captured by the formulation of our model. However, this also raises technical challenges due to the double continuous structure. Being able to identify age classes to follow in priority with contact tracing could be, though, an asset in controlling epidemic spread.
Finally, a challenge of optimal control theory remains its application in the field. We show in our results that technical obstacles related to the continuous nature of the strategy identified can readily be addressed by using discrete time periods. However, as mentioned above, the choice of the objective function remains delicate because its social implications. Being able to factor in multiple feed-backs (from clinicians, economists, but also the general public) could make optimal control approaches even more impactful in epidemiology.