key: cord-0021177-sr2l8swf authors: Bolzoni, Luca; Della Marca, Rossella; Groppi, Maria title: On the optimal control of SIR model with Erlang-distributed infectious period: isolation strategies date: 2021-09-22 journal: J Math Biol DOI: 10.1007/s00285-021-01668-1 sha: 559d20fccd435a461b9c02fe315fbd730634cc5c doc_id: 21177 cord_uid: sr2l8swf Mathematical models are formal and simplified representations of the knowledge related to a phenomenon. In classical epidemic models, a major simplification consists in assuming that the infectious period is exponentially distributed, then implying that the chance of recovery is independent on the time since infection. Here, we first attempt to investigate the consequences of relaxing this assumption on the performances of time-variant disease control strategies by using optimal control theory. In the framework of a basic susceptible–infected–removed (SIR) model, an Erlang distribution of the infectious period is considered and optimal isolation strategies are searched for. The objective functional to be minimized takes into account the cost of the isolation efforts per time unit and the sanitary costs due to the incidence of the epidemic outbreak. Applying the Pontryagin’s minimum principle, we prove that the optimal control problem admits only bang–bang solutions with at most two switches. In particular, the optimal strategy could be postponing the starting intervention time with respect to the beginning of the outbreak. Finally, by means of numerical simulations, we show how the shape of the optimal solutions is affected by the different distributions of the infectious period, by the relative weight of the two cost components, and by the initial conditions. Infectious diseases are responsible for significant social and economic losses worldwide. Several recent examples, one for all the ongoing pandemic caused by SARS-CoV-2 virus, pointed out the need for suitable public health interventions for effective disease control and eradication. Also in the case of livestock and wildlife, infectious diseases lead to huge economic losses to farmers (DEFRA/DCMS 2002) and endanger wild species' survival (McCallum and Dobson 1995) , in addition to representing potential risks for spillover in human (Daszak et al. 2000) . In the last 20 years, we have witnessed to an increase of the role of mathematical models in improving the understanding and making predictions on the infection dynamics and in informing the policies for disease management (Ferguson et al. 2001; Gatto et al. 2020; Gumel et al. 2004 ). The susceptible-infected-removed (SIR) class of models represent the most commonly used framework to investigate epidemiological systems through ordinary differential equations (Anderson and May 1992) . In SIR models, the host population units are subdivided into different compartments according to their epidemiological status as: susceptibles (S), i.e. healthy units that can be infected, infected (I ), i.e. diseased units that can infect other ones, and removed (R), i.e. units who died or recovered from the infection. The units of observation may represent single individuals or groups (e.g. farms in case of livestock diseases). The classical SIR models are built using different simplified assumptions on the epidemiological processes. One of the major simplifications consists in assuming a constant rate of leaving the infected class, namely an exponential distribution of the infectious period. Effectively, this means that the chance of a unit to recover during any given time period does not depend on the duration of the time that it has already been infected. Whilst such an assumption may provide significant mathematical convenience and be reasonably realistic in some situations, most often the assumption is violated, then requiring the inclusion of more realistic distributions of infectious periods in epidemic models (Bailey 1954; Gough 1977; Simpson 1952) . The consequences of assuming non-exponential distributions of the infectious period in SIR models on (i) the qualitative (Krylova and Earn 2013; Lloyd 2001a) and quantitative (Wearing et al. 2005; Feng et al. 2016) predictions of the infection dynamics, (ii) the estimate of epidemiological parameters (Wearing et al. 2005) , and (iii) the effectiveness of time-invariant intervention policies (Feng et al. 2007; Wang et al. 2017 ) have received some attention in recent years. However, no attempts have been made, to our knowledge, to investigate the consequences of assuming nonexponential distributions of the infectious period on the performances of time-variant disease control strategies by using optimal control theory (Pontryagin et al. 1962) . Here, by investigating the problem of minimizing the epidemic burden by using isolation control, we will show that different assumptions on the shape of the infectious period distribution may lead to qualitatively and quantitatively different shapes of the optimal control function, thus critically affecting the effectiveness of disease control policies. The paper is organized as follows. In Sect. 2, we present an SIR model with Erlangdistributed infectious period and isolation control (Erlang distributions are special cases of gamma distributions (Evans and Peacock 1993) ). An optimal control problem aimed at minimizing the control costs and the epidemic size is formulated in Sect. 3, where also necessary and sufficient conditions of optimality are derived. The profiles of admissible optimal isolation strategies are analytically characterized in Sect. 4. Then, we provide numerical simulations in order to understand the sensitivity of the optimal solutions to: (i) different assumptions on the distribution of the infectious period; (ii) relative cost of the two components of the objective functional; and, (iii) initial conditions (Sect. 5). Finally, discussion of the presented results and concluding remarks are given in Sect. 6. The inclusion of non-exponential distributions of the infectious period in epidemic models means that the chance of recovery depends on the time since infection, and hence the model needs to keep track of this information. Although this can be achieved in several ways, e.g. via an integro-differential equation formulation (Hethcote and Tudor 1980; Keeling and Grenfell 1997) or a partial differential equation formulation (Anderson and May 1992) , a significant body of research has been focused on the method of stages (see, just to mention a few papers, Cunniffe et al. 2012; Lloyd 2001a, b; Sherborne et al. 2015) , first introduced-to the best of our knowledgeby Cox and Miller (1965) . It consists in assuming that the infectious period obeys a distribution of the type Erlang(n, nγ ), where n ∈ N + and 1/γ > 0 is the mean duration of infection (the Erlang variate is a gamma variate with integer shape parameter (Evans and Peacock 1993) ). The corresponding probability density function is The popularity of this assumption comes from the fact that an Erlang distribution with shape parameter n is nothing else than the distribution of a sum of n independent exponentially distributed random variables of identical mean 1/(nγ ). The case n = 1 reproduces the exponential distribution; as the parameter n increases, the infectious period distribution becomes more closely centered on its mean, with n → +∞ corresponding to all units having exactly the same duration of infection. Thus, the use of the Erlang distribution in an epidemic model is equivalent to split the infected compartment in a finite number of fictitious stages in series: newly infected units enter the first stage, pass through each successive stage, and recover as they leave the nth stage. These multiple stages of infection can be also used to represent periods of increased or decreased risk of transmitting the disease (Ma and Earn 2006) . By adopting the method of stages, we consider the following SIR model with an Erlang(n, nγ )-distributed infectious period and including a bounded control function u(t) ∈ U = [0, u max ] representing the rate of isolation of infected units at time t: S = −β S I , (1) Here, the upper dot denotes the time derivative, S is the number of susceptible units at time t and I is the number of infected units at time t, which equals the sum of the infected units, I j , in each jth stage, with j = 1, . . . , n, namely The units of observation in model (1) may represent either individuals (in the case of human diseases) or farms (in the case of livestock diseases). The parameters are positive constants: β represents the transmission rate and γ represents the removal rate. The assumption of upper boundedness for u(t) reflects the realistic scenario of practical limitations on the maximum rate at which the control strategy may be employed. Such a restriction is due to several factors, including the limited number of health officers, the time necessary for detecting the target individuals or farms, the reduced capacity of infrastructures and so on. Specifically, by setting an upper bound on u(t), we are assuming a limitation on the maximum isolation effort relative to the target population instead of on the overall isolation activities performed per time unit (which are represented by u(t)I ). For convenience of notation, let us denote μ = nγ and adopt the following change of variables: namely, Y k is the cumulative of the infected units in the first k stages of infection. In particular, Y n = I . One can easily check that system (1) becomeṡ (2) Then, the general controlled system in compact form reads: where x = (S, Y 1 , . . . , Y n ) T is the state-variables vector, and f and g are given by We denote with T the terminal time of planning (i.e. time horizon). In general, it is not clear how large T must be in order to coincide with a reasonable interpretation of the end of the epidemic. To overcome the problem, by following the approach adopted by Bolzoni et al. (2019) , Bolzoni et al. (2017) and Hansen and Day (2011) , we choose T as the disease extinction time. Formally, let ε be a positive constant such that ε < 1, then which implies that the optimal control problem will have free final time. To make the definition (4) meaningful, we choose the initial number of infected units Y n (0) = I (0) strictly greater than ε. As a consequence, T is the first time at which the state variable Y n drops to ε. Using this approach, we also avoid the undesirable possibility of subsequent peaks of infection caused by a fractional number of infected units (Hansen and Day 2011) . From a biological point of view, by setting ε < 1 as the condition to identify the extinction time means assuming the infection goes extinct when there is less than one infected unit in the host population. To model equations (2), we associate the following epidemiologically well-posed initial conditions: Note that, by simply looking at the first and the last equation in (2), one gets that the set is positively invariant, i.e. solutions of (2) starting in D remain in D for all t ≥ 0. Therefore, it is not restrictive to limit our analyses to region D. This section is devoted to the formalization of our control problem. Specifically, we set an optimal control problem involving two different and coexisting goals, namely minimization of: • costs due to the isolation efforts per time unit: T 0 A 0 u(t)dt; • costs due to the number of individuals or farms that have been infected during the whole outbreak: where A 0 , A 1 ∈ R + are the balancing weight factors, representing the costs per time unit of isolation efforts and the costs of a single new infection, respectively. Of course, A 1 must be strictly positive for the problem being meaningful (otherwise, we are searching for the control that minimizes its own cost). Hence, it is not restrictive to set A = A 0 /A 1 ≥ 0, with A representing the relative cost per time unit of isolation efforts over the cost of a single new infection. By setting the costs of isolation per time unit as A 0 u(t), we are assuming that the control costs are proportional to control effort instead of to the overall control activities, which are represented by u(t)Y n . This assumption fits the scenarios in which the procedures of isolation are forerun by surveillance and testing activities aimed at detecting the target populations, as in the cases of test-and-isolate and test-and-cull control strategies (Bolzoni et al. 2019) . In sum, we are considering the following optimal control (OC) problem, aiming at minimizing simultaneously control efforts and epidemic size: (OC) Find an optimal control function u * (t) that minimizes the objective functional with T given by (4), on the admissible set subject to (2) with initial conditions (5). The existence of the solution for the OC problem is guaranteed since the requirements of classical existence theorems are satisfied. More precisely, we refer to the formulation of Filippov theorem and its corresponding corollary given by Agrachev and Sachkov (2004) for proving the following theorem. Theorem 1 There exists an optimal control solution u * (t) to the OC problem. Proof We list the requirements of the corollary of Filippov existence theorem in Agrachev and Sachkov (2004) : (6)), are convex; 3. there exists a constant C such that: Condition 1 is trivial. Since x ∈ D and u(t) ∈ U , the state and the control variables are a priori bounded and also Condition 3 is verified. To check that Condition 2 holds, firstly we take x ∈ D and prove that the line connecting two points of A U (x) entirely lies in A U (x). Let us take y 1 , y 2 ∈ A U (x), then there exist u 1 (t), u 2 (t) ∈ U such that y j = f(x) + u j (t)g(x), j = 1, 2. One easily gets: As regards the characterization of the optimal control, according to Pontryagin's minimum principle (Chachuat 2007; Pontryagin et al. 1962; Schmitendorf 1976) , the problem of finding the time-dependent control variable u * (t) that minimizes the objective functional (7) is equivalent to the problem of minimizing the Hamiltonian: given by the sum of the integrand of (7) and the scalar product of λ and the right-hand side of the state system (2). Here λ = λ S , λ Y 1 , . . . , λ Y n T denotes as usual the adjoint variables vector at time t, which is the solution of the adjoint systeṁ with the transversality conditions The notation with superscript * is used to denote the solution to the state system (2) corresponding to the optimal control u * (t). The above transversality conditions are consequence of the equality terminal constraint on the state variable Y n given by (4). For analytical details see Chachuat (2007), Pontryagin et al. (1962) and Schmitendorf (1976) . Further, H is constantly equal to zero along the optimal solution: In sum, the optimal control u * (t) that solves the control problem is the time-dependent function that minimizes the Hamiltonian (9) on the admissible set (8): and the control horizon T is defined by (4). Since H is linear respect to u(t), the value of u * (t) is determined by the sign of the so-called switching function: In particular, If ψ(x * , λ) vanishes on a set of isolated points, the corresponding control function is a bang-bang control, as well known, and its value at these points is not defined. However, for convenience of numerical simulations, we set it equal to u max . The isolated time instants, internal to [0, T ], in which ψ vanishes are called switching times. Unless otherwise specified, we restrict the use of the superscript * only to the optimal control function u * (t), in order to simplify the notation. The main result is that the admissible optimal controls for the OC problem under investigation are bang-bang with at most two switching times. The optimal solution u * (t) to the OC problem is a bang-bang control. with τ 1 , τ 2 ∈ [0, T ], and τ 1 ≤ τ 2 . Proof Throughout the proof we omit the superscript * for the optimal quantities, in order to simplify the notation. First we compute the Hamiltonian (9) along the optimal solution: and explicitly write the adjoint system (10): The switching function with its first and second time derivatives read, respectively: Let us prove that ψ can vanish only in isolated points. By contradiction, suppose that ψ vanishes on some interval I. Then, also all its time derivatives vanish there; in particular, ψ =ψ =ψ = 0 on I yields: implying n k=1λ Y k By explicitly writing the addenda of (21), one obtains n k=2 λ Y k = 0 on I, that, along with (20), gives λ Y 1 = 0 on I. Then, alsoλ Y 1 = 0 on I and, by proceeding in series, we have (see the corresponding adjoint equations in (15)). Two cases must be considered: • A > 0. Then, equalities (22) are in contrast with (18). (18)-(22) it follows that, once ψ vanishes, it must remain null, namely T ∈ I. However, the equality (19) is in contrast with the transversality condition λ S (T ) = 0. In both cases a contradiction arises. As a consequence, ψ can vanish only in isolated points, namely u(t) is a bang-bang control: it can assume only the extreme values u(t) = 0 and u(t) = u max . Let us now prove that u(t) changes its value at most two times and, in such a case, the two switches are from u(t) = 0 to u(t) = u max and again back to u(t) = 0. As first step, we consider the value assumed by the quantities involved in the OC problem at the extinction time T . Taking into account the transversality conditions (11), we have H From the adjoint equations (15), one yieldṡ Substituting (25) into the formula of the Hamiltonian (14) provideṡ . . , n − 2. Take k = n − 2 and assume, by contradiction, that λ Y n−2 (T − ) < 0 <λ Y n−2 (T − ). By handling the corresponding adjoint equation (see (15)), one yields that is not possible because λ Y n−1 (T ) = 0 andλ Y n−1 (T ) < 0. Thus, a contradiction arises and it must be λ Y n−2 (T − ) > 0 >λ Y n−2 (T − ). By applying the same argument, one can proceed in series for 1 ≤ k < n − 2. As a further step, we show that λ Y k ≥ 0 in [0, T ], ∀k = 1, . . . , n. To this end, let us definet = inf A, with T ], ∀k = 1, . . . , n} and prove thatt = 0. From the above considerations it follows that T ∈ A = ∅ and t < T . By contradiction supposet > 0, namely there existsk ∈ {1, . . . , n} such that λ Yk (t) = 0 <λ Yk (t). Two cases must be considered: •k < n. Then, λ Yk +1 (t) = −λ Yk (t)/μ < 0 (see (15)), which is in contrast with the definition oft. •k = n. Then, sgn(λ S (t)) = sgn(λ Y n (t)) > 0. Sinceλ S (T ) < 0 as given in (24), there existst ∈ A,t 0 >ψ(τ s2 ). Sinceψ(T ) > 0 (see (11)- (17)),ψ vanishes at least two times in (0, T ), say them t M ∈ (τ s1 , τ s2 ) and t m ∈ (τ s2 , T ), that are, respectively, relative maximum and minimum points for ψ. From (17) it follows that λ S (t M ) = 1 = λ S (t m ), which, however, contradicts the condition (29). As a consequence, the presence of three or more switching times is definitively excluded. Also, if u(t) changes its value two times, it is necessarily from 0 to u max and again back to 0. This concludes the proof. In Theorem 2, the time instants τ 1 ∈ [0, T ] and τ 2 ∈ [τ 1 , T ] are introduced to extend the concept of switching time. Specifically, τ 1 [resp. τ 2 ] is the starting [resp. ending] intervention time, i.e. the first [resp. last] time instant at which the control u * (t) assumes the value u max . The knowledge of such time instants allows to immediately obtain the control time-profile. Indeed, except for the case u * (t) ≡ 0 when the best strategy is to 'do nothing', from Theorem 2 it follows that the following four possible control profiles are allowed: (i) τ 1 = 0 and τ 2 = T , i.e. u * (t) ≡ u max , (ii) τ 1 > 0 and τ 2 = T , i.e. τ 1 is the only switching time and u * (t) changes its value from 0 to u max (say, delayed control), (iii) τ 1 = 0 and τ 2 < T , i.e. τ 2 is the only switching time and u * (t) changes its value from u max to 0 (say, reactive control), (iv) τ 1 , τ 2 ∈ (0, T ) are two switching times and u * (t) changes its value from 0 to u max and again back to 0. We also prove that the control profile (i) is the only admissible if one wants to minimize just the epidemic size and not also the isolation efforts, namely A = 0 in the objective functional (7). (7), then the optimal solution u * (t) to the OC problem is constantly equal to u max . Proof Let us resume the proof of Theorem 2 and consider the case A = 0. Then, the switching function reads (23)). From (12), it follows that u(T ) = u max . By contradiction, let τ s be a switching time for implying that there exist at least two indexes k 1 , k 2 ∈ {1, . . . , n} such that λ Y k 1 (τ s ) < 0 < λ Y k 2 (τ s ). However, this in contrast with the sign of the adjoint variables (28). As a consequence, it must be u(t) ≡ u max , which is the assert. As far as the numerical simulations are concerned, the shape of admissible optimal controls allows us to adopt a simple ad hoc numerical scheme, which generalizes those implemented by Bolzoni et al. (2017 Bolzoni et al. ( , 2019 . Such a scheme is based on the idea of identifying each bang-bang solution u * (t) with two real parameters: the starting and ending intervention times. Therefore, the functional J (u) to be minimized can be seen as the function that links the pair of characteristic intervention times to the objective integral (7) of the OC problem: An optimal control u * (t) will be identified by the pair (τ 1 , τ 2 ) such that (τ 1 , τ 2 )=argminJ . The numerical solution will be computed by evaluating the function J (τ 1 , τ 2 ) over a suitable interval and looking for its minimum value. Of course, if the minimum value of J is greater than the corresponding value in absence of control (i.e. T 0 β SY n dt as predicted by system (2) with u(t) ≡ 0), then the method selects for u * (t) ≡ 0 as the optimal isolation policy. For the numerical integrations of the Cauchy problem (2)-(5) the fourth-order Runge-Kutta scheme with uniform time steps is used. For a detailed discussion of the implemented numerical method see Bolzoni et al. (2017) . Under the assumption of different shapes in the distribution of the infectious period in model (2), we numerically investigate the optimal isolation strategies focusing on highly contagious livestock diseases, such as foot and mouth disease (FMD) and highly-pathogenic avian influenza (HPAI). However, thanks to the generality of the proposed framework, the obtained results can be qualitatively applied to different contexts, like that of non-pharmaceutical interventions (such as self-isolation) in human epidemics. As in Bolzoni et al. (2019) , we assume the farm as the unit of observation, since the control measures in livestock diseases are generally implemented at farm level (Probert et al. 2016) . Moreover, our choice to work at farm level is strengthen by the finding obtained by Ferguson et al. (2001) , who estimated that the infectious period of a farm infected by FMD can be described through an Erlang distribution with n 1. The initial conditions for model (2) variables are set to a single infected unit (Y k,0 = 1, ∀k = 1, . . . , n) in a totally susceptible population (S 0 = 2000). By setting Y k,0 = 1, ∀k = 1, . . . , n, we assume that the infection is introduced in the population through a newly infected farm, i.e. I 1,0 = 1 and I k,0 = 0, k = 2, . . . , n. While this assumption seems obvious when the epidemiological units are represented by farms, other assumptions on the initial conditions can be reasonable when the epidemiological units are represented by individuals (e.g. in human diseases). Thus, we also test the effects on the optimal control of different assumptions on the initial conditions. Specifically, we derive the optimal solutions by varying the stage of infection of the first infected unit (say, k 0 ), which corresponds to assume Y k,0 = 0, 1 ≤ k < k 0 and Y k,0 = 1, k 0 ≤ k ≤ n. We set the infectious period, corresponding to the average time a unit stays infected before recovering, to 1/5 months. Thus, the recovery rate (γ = μ/n) in model (2) is set to γ = 5 month −1 . The chosen value of infectious period falls into the range of those obtained both for FMD and HPAI (Boender et al. 2007; Bouma et al. 2003; Haydon et al. 1997; Rossi et al. 2017; Sharkey et al. 2008) . We set the basic reproduction number (R 0 ), which represents the average number of secondary infections produced by a single infected in a completely susceptible population (Hethcote 2000) , to R 0 = 4 (Bouma et al. 2003; Haydon et al. 1997; Sharkey et al. 2008; Stegeman et al. 2004 ). Since R 0 in model (2) can be written as R 0 = β S 0 /γ (Hethcote 2000) , we set β = 0.01 farm −1 month −1 . As in Bolzoni et al. (2019) , we set the maximum isolation effort to u max = 1 month −1 and ε = 0.5. Hansen and Day (2011) pointed out that any value of ε smaller than one fits the purpose of the analysis, while Bolzoni et al. (2019) showed in a similar framework that the shape of the optimal control does not significantly depend on the chosen value of ε. In our simulations, we also test the sensitivity of the optimal control solutions to different combinations of the number of stages of infection n and the control parameter A (which represents the ratio between the costs per time unit of the isolation effort and the costs of a single new infection). Hethcote and Tudor (1980) observed that changing the distribution of the infectious period does not lead to significant differences in the asymptotic behavior of SIR-like epidemic models. However, it can crucially affect the quantitative predictions of the epidemic model on disease burden, epidemic duration, and effectiveness of disease control strategies, as shown by Wearing et al. (2005) . In Fig. 1 , we display the numerical simulations of model (2) infected dynamics in the absence of control (i.e. u(t) ≡ 0). Figure 1a shows the dynamics of the infected state variables Y k , k = 1, . . . , n, as a function of time (t), for n = 1 (black curve) and n = 10 (gray curves). (Thick curves represent the total number of infected, Y n .) In Fig. 1b , we display the values of the extinction time (T ), the time at which the epidemic peak occurs (t p = argmax(Y n )), and the total number of infected at the peak (max(Y n ), see the box) obtained in the simulations as functions of the number of stages of infection n. Similarly to previous results obtained with SIR-like epidemic models (Wearing et al. 2005) , our numerical analyses on model (2) in the absence of control confirm that the shape of the infectious period distribution quantitatively affects the model predictions. Specifically, an Erlang-distributed model (n > 1) predicts a faster initial increase of the infected curve, a larger peak of infection, and a faster extinction time with respect to the exponential-distributed one (n = 1). From n = 1 to n = 20, In Fig. 2 (resp. Fig. 3) , we show the numerical solutions of the OC problem for four values of the number of stages of infection: n ∈ {1, 5, 10, 20}, by setting the relative cost in the objective functional (7) equal to A = 0.1 [resp. A = 10]. Gray curves represent the temporal dynamics of the infected state variables Y k , k = 1, . . . , n, black curves those of the optimal isolation strategies u * (t). Figures 2 and 3 show that, analogously to that observed in the uncontrolled model, different assumptions on the shape of the infectious period distribution quantitatively affect the optimal control solutions. Specifically, OC models assuming exponentially distributed infectious periods tend to overestimate the costs associated to control (we estimate an increase of about 50% in the case of n = 1 with respect to n = 10, see Figs. 2, 3) . In addition to the quantitative differences described above, our numerical analyses show that different assumptions on the distribution of the infectious period can affect qualitatively the solutions of the OC problem. In Fig. 4 , we display the values of the starting (τ 1 ) and ending (τ 2 ) intervention times and the extinction time (T ) for the optimal solutions as functions of the relative cost A, by setting n = 1 (black curves) and n = 10 (gray curves). From Fig. 4 , we notice that for a wide range of A values (approximately from A = 2 · 10 −2 to A = 10 −1 ) the optimal control problem selects for a constant control u * (t) ≡ u max (i.e. τ 1 = 0 and τ 2 = T ) when n = 1, while the optimal control admits a one-switch reactive solution (i.e. τ 1 = 0 and τ 2 < T ) when n = 10. Furthermore, for more limited ranges of A, we show that the optimal control problem can: (i) admit a one-switch solution when n = 1, while it admits a twoswitches solution when n = 10 (see A ≈ 10 −1 in Fig. 4) ; (ii) admit a two-switches solution when n = 1, while it selects for u * (t) ≡ 0 when n = 10 (see the black and gray dash-dotted lines in Fig. 4) . Interestingly, for both n = 1 and n = 10, the starting and ending intervention times are monotone (increasing and decreasing, respectively) functions of A, while the disease extinction time T is increasing for a wide range of A and then starts to decrease when A approaching the value above which u * (t) ≡ 0. Finally, we also consider variations with respect to the initial conditions of the OC problem. Specifically, we set the relative cost A = 10, the number of stages of infection n = 10, and vary the stage entered by the first infected unit (k 0 ). Namely, Y k,0 = 0, 1 ≤ k < k 0 , Y k,0 = 1, k 0 ≤ k ≤ n. The values of the starting (τ 1 ) and ending (τ 2 ) intervention times and the extinction time (T ) for the optimal solutions are displayed in Fig. 5 . We notice that different assumptions on the initial class k 0 of the first infected do not substantially change the shape of the optimal solution. Both the intervention times and the extinction time slightly increase with k 0 : from k 0 = 1 to k 0 = 10 the growth is approximately 0.065 months. In this work, we investigate optimal isolation strategies in an SIR model with Erlang distribution of the infectious period. We adopt the method of stages (Cox and Miller 1965) , by splitting the infected compartment in n fictitious stages in series (n ∈ N + ); infected units at any stage are isolated with time-variant bounded rate u(t) (here, units Parameters n = 10, A = 10, Y k,0 = 0, 1 ≤ k < k 0 , Y k,0 = 1, k 0 ≤ k ≤ n; unspecified parameters as in Sect. 5.1 represent individuals or farms). The aim of the optimal control (OC) problem is to jointly minimize the costs associated to the control efforts and the epidemic size. By characterizing optimal control shapes according to Pontryagin's minimum principle (Pontryagin et al. 1962) , we prove that, for both n = 1 (exponential distribution of the infectious period) and n > 1 (non-exponential distribution of the infectious period), the OC problem admits only bang-bang solutions with at most two switches. When the optimal control u * (t) is constant for the whole time horizon [0, T ], it can be either null or equal to its maximum value u max . When only one switch is admitted, the OC problem can select for both reactive isolation (characterized by u * (0) = u max and u * (T ) = 0) and delayed isolation (characterized by u * (0) = 0 and u * (T ) = u max ). When the optimal solution admits two switches, then the OC problem selects for controls where isolation is implemented within a temporal window between times τ 1 and τ 2 (with τ 1 < τ 2 and τ 1 , τ 2 ∈ (0, T )). Numerical analyses show that increasing in the number of stages of infection n is positively associated to faster, but with higher peaks, epidemic outbreaks both in absence and in presence of control. When optimal isolation strategies are implemented, even if, a priori, the spectrum of admissible control profiles is the same for any n, the optimal solutions vary both quantitatively and qualitatively. Specifically, the OC problem in the case n = 1 can select for solutions with a different number of switches (if any) with respect to the case n > 1. In addition, we observe that the control costs associated to the epidemic significantly decrease with n, as a consequence of the reduced epidemic lengths. We also test the sensitivity of the optimal control solutions to the weight parameter A in the objective functional, that is the ratio between the costs per time unit of the isolation effort and the costs of a single new infection. From numerical simulations (see Fig. 4 ), we observe that, by increasing A, u * (t) passes from a control constantly equal to u max (that is the unique admissible solution in the special case A = 0, see Corollary 3), to a reactive control up to a two-switches control. Lastly, when A is sufficiently large, the best strategy is to 'do nothing' (i.e. u * (t) ≡ 0). While it seems unrealistic to expect that stakeholders and policy-makers can express their preferences on specific values of parameter A, the outputs of the optimization problem can identify the parameter regions where the benefits of a control campaign outweigh the costs and those where the costs outweigh the benefits, thus helping stakeholders and policymakers in their evaluations. To the best of our knowledge, our work is the first attempt to use an SIR model with non-exponential distribution of the infectious period in an optimal control framework. As is widely known, for many infectious diseases the infectious period distributions are not exponential. Hence, models with exponential distributions may lead to biased outcomes in both the understanding of behaviors of disease dynamics and the evaluations of disease control programs. The Erlang distribution proved to be a valuable compromise between realism and analytical tractability, and was largely used to estimate infectious period distributions by fitting with empirical data (see, just to mention few papers, Eichner and Dietz 2003; Ferguson et al. 2001; Wearing et al. 2005 ). This work serves as a proof-of-principle verification of the impact on optimal isolation strategies of different assumptions about the distribution of the infectious period. Our analytical and numerical results stress the conclusion that simple epidemic models, such as the ones assuming exponential distribution on the infectious period or homogeneous mixing in the host population, can not provide reliable solutions to real-world cases of disease control from the quantitative point of view. However, the theoretical findings provided by these models can be used to inform more complex simulation models developed for specific epidemiological scenarios where more realistic descriptions of the biological, ecological and epidemiological processes are included. We investigate the specific control problem where the costs of control are proportional to isolation efforts instead of to the overall number of units isolated. This assumption fits the scenarios where the costs linked to the activities of surveillance and testing aimed at detecting the units to isolate (which depend on the effort and not on the target population) are larger than the costs of isolating the infected units. We find that the optimal solutions obtained under this assumption substantially differ from those obtained for optimal control problems where the control costs are proportional to the number of infected units (Morton and Wickwire 1974; Wickwire 1975) . Roughly speaking, there are two alternative formulations for the cost of the control campaign in the objective functional: (i) control costs proportional to the overall control activities, that is the control rate multiplied by the units to whom it is addressed (e.g., in our notations, u(t)S in case of vaccination, u(t)I in case of isolation or treatment, see e.g. Behncke 2000; Morton and Wickwire 1974; Wickwire 1975) ; (ii) control costs proportional to control efforts (that is the control rate u(t) as considered here, see e.g. Behncke 2000; Bolzoni et al. 2014; Zhou et al. 2014) . The assumption (i) was largely adopted in OC applications to the SIR model under different control strategies. In such a case, the problems only support the trivial solutions of applying the maximal effort for the entire epidemic (u * (t) ≡ u max ) (Wickwire 1975) or having at most one switch from u * (t) = u max to u * (t) = 0 (reactive control) (Behncke 2000; Morton and Wickwire 1974) . These are also the only admissible solutions for optimal vaccination policies when the assumption (ii) is adopted (Zhou et al. 2014 ). Conversely, we prove here that optimal isolation policies under the assumption (ii) are not necessarily active since the initial time of the epidemic. Namely, the optimal starting intervention time could be postponed, which yields delayed controls or even two-switches controls (depending on if the ending intervention time coincides or not with the disease extinction time). The biological explanation for the optimality of delayed [resp. two-switches] control when the costs are proportional to the control efforts relies on the remark that, when the maximum effort admissible (i.e. u max ) is not large enough to significantly flatten the epidemic curve, employing resources to find infected units at the beginning [resp. the beginning and the end] of the epidemic may not be efficient, since several susceptible units would be screened/tested in the face of few infected ones. The described scenario is similar to the case where the disease control is applied only when the number of infected (or new infections) exceeds a given threshold, which can be mathematically investigated through sliding models (Bolzoni et al. 2020; Qin et al. 2016 ). Encyclopaedia of mathematical sciences Infectious diseases of humans: dynamics and control A statistical method of estimating the periods of incubation and infection of an infectious disease Optimal control of deterministic epidemics Risk maps for the spread of highly pathogenic avian influenza in poultry React or wait: which optimal culling strategy to control infectious diseases in wildlife Time-optimal control strategies in SIR epidemic models Optimal control of epidemic size and duration with limited resources Dynamics of a metapopulation epidemic model with localized culling The foot-and-mouth disease epidemic in The Netherlands in 2001 Time-dependent infectivity and flexible latent and infectious periods in compartmental models of plant disease Emerging infectious diseases of wildlife-threats to biodiversity and human health Joint Working Paper of the Department for the Environment, Food and Rural Affairs and the Department of Culture Transmission potential of smallpox: estimates based on detailed data from an outbreak Epidemiological models with non-exponentially distributed disease stages and applications to disease control Mathematical models of Ebola-consequences of underlying assumptions The foot-and-mouth epidemic in Great Britain: pattern of spread and impact of interventions Spread and dynamics of the COVID-19 epidemic in Italy: effects of emergency containment measures The estimation of latent and infectious periods Modelling strategies for controlling SARS outbreaks Optimal control of epidemics with limited resources An analysis of foot-and-mouth-disease epidemics in the UK The mathematics of infectious diseases Integral equation models for endemic infectious diseases Disease extinction and community size: modeling the persistence of measles Effects of the infectious period distribution on predicted transitions in childhood disease dynamics Destabilization of epidemic models with the inclusion of realistic distributions of infectious periods Realistic distributions of infectious periods in epidemic models: changing patterns of persistence and dynamics Generality of the final size formula for an epidemic of a newly invading infectious disease Detecting disease and parasite threats to endangered species and ecosystems On the optimal control of a deterministic epidemic The mathematical theory of optimal processes. International series of monographs in pure and applied mathematics Effects of limited medical resource on a Filippov infectious disease model induced by selection pressure Modelling farm-to-farm disease transmission through personnel movements: from visits to contacts, and back Pontryagin's principle for problems with isoperimetric constraints and for problems with inequality terminal constraints Epidemiological consequences of an incursion of highly pathogenic H5N1 avian influenza into the British poultry flock Dynamics of multi-stage infections on networks Infectiousness of communicable diseases in the household (measles, chickenpox, and mumps) Avian influenza A virus (H7N7) epidemic in The Netherlands in 2003: course of the epidemic and effectiveness of control measures Evaluations of interventions using mathematical models with exponential and non-exponential distributions for disease stages: the case of Ebola Appropriate models for the management of infectious diseases Optimal isolation policies for deterministic and stochastic epidemics Optimal vaccination policies for an SIR model with limited resources Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations We are thankful to the anonymous reviewer for the useful comments. The work of LB was supported by the Italian Ministry of Health (Grant IZSLER PRC2014003). The work of MG was supported by the Italian Ministry of Education, University and Research (Prin 2017YBKNCE). Support by INdAM-GNFM and by University of Parma is also gratefully acknowledged by RDM and MG. This research benefited from the HPC (High Performance Computing) facility of the University of Parma.Funding The work of LB was supported by the Italian Ministry of Health (Grant IZSLER PRC2014003). The work of MG was supported by the Italian Ministry of Education, University and Research (Prin 2017YBKNCE). The authors declare they have no competing interests. Code availability Not applicable.