1 Introduction

Coronavirus was treated as a non-fatal virus till 2002. Still, after its transmission in various countries like the USA, Singapore, Hong Kong, Taiwan and severe respiratory syndrome and deaths reported in 2003, it is also considered as a deadly virus [3]. COVID-19, the respiratory disease caused by the SARS-CoV-2 strain of coronavirus has prompted a global pandemic, since the initial outbreak in Wuhan, China in 2019. According to an estimate by WHO, in the early stages, the disease is less contagious than smallpox, chickenpox, and mumps but more contagious than similar conditions like seasonal flu and MERS. The estimated disease fatality rate for COVID-19 is 3.4%, which is less than that of smallpox, SARS and MERS [7]. As of now, almost every country has been affected by this disease, with 14.7 million infected people across the world with more than 6,00,000 deaths. WHO considered Europe as the active center of this pandemic on 13 March 2020 [29]. On 24 January 2020, the first COVID-19 case of Europe was confirmed in France [22]. Currently, we observe countries like the U.K., Italy, Spain, France still have high infection rates, whereas Germany, Croatia, Iceland seem to have successfully mitigated the disease.

As human to human virus transmission chain has been established, the primary route of transmission is through droplet spread [12]. As there is no well-proven treatment to cure the disease and lack of COVID-19 vaccine to build immunity against the disease, non-pharmaceutical interventions like social distancing emerge as an unprecedented measure taken by most countries followed by isolation of affected individuals by testing methods. But how does the virus run out of the people to infect? It could either be due to the lockdown and social distancing measures, which prevent the virus from spreading, or it may be because the disease has raced through a population, infecting a big chunk of it and leaving behind a small portion that does not have the virus, the case of herd immunity scenario. While restrictions and bans introduced to deal with the outbreak have been in place for months, concerns that reopening will bring back a second wave have not been universally borne out, and confirmed cases remain a fraction of the total population in most places. But the study in such areas suggests that more people may have had the disease than officially recorded [8, 25].

From the perspective of understanding the dynamics of the pandemic, several disease modeling approaches have been tried. Compartmental models are the most straightforward way out to study an epidemic. As the disease dynamics depend upon contact between people, an age-structured model can help provide critical insights, which can be instrumental in developing policies to control the outbreak. This model can be used to study the effect of physical distancing measures ranging from staying at home, leaving home for essential services, closure of schools and colleges, and reopening non-essential sectors. Some authors tried to show up the closure effects in resetting the number of susceptible people with a good old SIR model. Still, this oversimplification does not allow to understand measures taken in different fields like school, work or other activities [4]. Delay differential equations had also been used in predicting future conditions along with an examination of health care capacity in the pandemic in vaccine-free scenario [37].

As already mentioned in the case of COVID-19, the means of transmission were unknown, and a vaccine was not available. In such situations, health authorities have recommended isolating individuals who may be infected. In the absence of other means to curb the spreading of a disease, the only way to slow down its propagation is to deny possible infection pathways. Adding quarantine compartment to classical SIR model facilitates us to consider the real-world scenario where a set of people can be kept in isolation and can not spread the disease. Quarantine measures play an essential role in preventing human diseases and epidemics, such as smallpox, tuberculosis, SARS, and the current outbreak of COVID-19. Ruschel et al. [23] used the SIQ model to study the impact of quarantine to control the spread of infectious diseases. Cao et al. [2] used the SIQR model with quarantine to calculate the system’s different transmission thresholds to regulate the disease dynamic. In this paper, we have formulated an age-structured SIQR model and fitted it with case data from France, Spain and Germany separately. This paper is designed as follows. In Sect. 1, we introduce the mathematical structure of the SIQR model and discuss the theoretical perspective of the model along with justification of model assumptions and their limitations. In Sect. 2, we have performed a mathematical analysis of the model and established its soundness. In Sect. 4, we describe the methodology data sources and discuss the parameters. In Sect. 5, we analyze and interpret the predictions generated by numerical simulations of this model. Finally, in Sect. 6, we concluded with implementation, potential shortcomings, and future aspects of our model with further improvements.

2 Model formulation and its limitations

Starting from classical SIR model, we have added quarantined population (Q) compartment in our model to keep track of isolation. In this paper, we propose a age-structured model and have divided our population into 16 age-groups of five years each. Mathematical framework of our model depends upon 64 coupled ordinary differential equations (ODEs) accounting for four ODEs for each age-group. The total population N(t) is divided into the following 4 compartments.

  1. 1.

    Susceptible Population, S(t): These individuals are not yet infected but may become infected or remain susceptible depending upon how they mix up with Infected individuals.

  2. 2.

    Infected population, I(t): These individuals are infected, but either pre-symptomatic (they haven’t started exhibiting symptoms) or asymptomatic. These are untested individuals, capable of spreading the infection.

  3. 3.

    Quarantined population, Q(t): These individuals are symptomatic and tested. They are removed from contact with the population, so not capable of spreading the disease.

  4. 4.

    Recovered population, R(t): These individuals have recovered from the disease and assumed to be immune to reinfection.

Thus, \(N(t) = S(t) + I(t) + Q(t) + R(t)\). A fifth compartment, F(t), which contains all disease-induced fatalities, is introduced to keep track of deaths and to determine mortality rate.

Each compartment is further subdivided into M age classes, indexed \(S_i\), \(I_i\), \(Q_i\), \(R_i\) and \(F_i\) for \(i = 1\ldots M\). The degree of contact between two such age classes is determined by the coefficient \(C_{ij}(t)\). This is further broken down into contact contributions from home (H), work (W), schools (S), and other sources (O), so the social contact matrix \([C_{ij}]_{M\times M}\) is given by

$$\begin{aligned} C(t) \;=\; \alpha _H(t) C_H + \alpha _W(t) C_W + \alpha _S(t) C_S + \alpha _O(t) C_O. \end{aligned}$$
(1)

The coefficients \(\alpha _m\), \(m \in \{H, W, S, O\}\) vary with time as social distancing measures change the degrees of contact in these different places. For example, \(\alpha _S\) is set to zero when all educational institutes are closed. Similarly, \(\alpha _m\) are reduced during periods of partial or total lockdown and fitted to existing data. We note that contributions from work and other sources never drops to zero due to the continued operation of essential services. Lockdowns also increase the degree of contact within a household [14].

Thus, we develop the parameter

$$\begin{aligned} \lambda _i(t) \;=\; \sum _{j = 1}^M C_{ij} \frac{I_j}{N_j}, \end{aligned}$$
(2)

which regulates the degree of contact of age class i with infected individuals from all ages.

Susceptible individuals in compartment S only become infected upon contact with an already infected individual, thus moving into the infected compartment I. We assume that a fraction a of such individuals are asymptomatic. Thus, the symptomatic fraction of infected individuals, \((1-a)\) move into quarantine Q at a rate \({\bar{\delta }}\), and the asymptomatic fraction a recover and move into R at a rate \({\bar{\epsilon }}\). Quarantined individuals also recover and move into R at a rate \(\gamma \).

Note that we have assumed that symptomatic individuals are moved immediately into quarantine, while asymptomatic individuals are left behind and proceed straight to recovery. Thus, only individuals in compartment Q suffer disease induced mortality at a rate \(\mu \). There is no further natural death rate, aging, or recruitment of the population.

Fig. 1
figure 1

Schematic illustration of movement between compartments for the model (3)

Full size image

The movement of individuals between compartments has been illustrated in Fig. 1. The dynamics of this model obeys the following system of ordinary differential equations.

$$\begin{aligned} \begin{aligned} {\dot{S}}_i(t) \;&=\; -\beta \lambda _i S_i, \\ {\dot{I}}_i(t) \;&=\; \beta \lambda _i S_i \,-\, ((1 - a)\bar{\delta } \,+\, a\bar{\epsilon })I_i, \\ {\dot{Q}}_i(t) \;&=\; (1 - a)\bar{\delta }I_i \,-\, (\gamma \,+\, \mu _i)Q_i, \\ {\dot{R}}_i(t) \;&=\; a{\bar{\epsilon }} I_i \,+\, \gamma Q_i, \\ {\dot{F}}_i(t) \;&=\; \mu _i Q_i. \end{aligned} \end{aligned}$$
(3)

We can absorb the terms \(\delta = (1-a){\bar{\delta }}\) and \(\epsilon = a{\bar{\epsilon }}\). The initial conditions are set as \(S_i(0) > 0\), \(I_i(0) > 0\), \(Q_i(0) > 0\), \(R_i(0) = 0\) and \(F_i(0) = 0\). All parameters in the system are positive quantities, and have been tabulated in Table 1.

Table 1 Parameters and their values
Full size table

We have considered 50% of the cases to be asymptomatic, which have the potential of spreading the disease. This fraction has been assumed based on studies such as [9, 10, 36], but this parameter is difficult to determine precisely as it may vary significantly between different regions or countries. Moreover, many people reportedly came into the contact of an infected person had also been put under home quarantine. We have not considered any spontaneous quarantine in our model. We might note that in actual scenarios, mortality rate of asymptomatic people can be far lesser but in our model for simplicity we have considered the same mortality rate for symptomatic and asymptomatic across all age-groups. As per study, mortality rate depends upon the pre-existing medical condition and mostly on age. The older people face a high mortality rate for having lower immunity. We have obtained different value of mortality rate for each age-group by fitting the curve, though the later part falls fast suggesting medical institutions are able to gain experiences of treating the disease, giving rise to higher recovery rate.

Our model is not gender structured, which is another major limitation as it has been suggested males are more likely to develop critical infections than female due to observed gender bias in mortality [24]. As of now, we have not seen enough evidence of herd immunity [18]. In all of the three countries, disease has spend more than five months. Our model suggests second wave of the pandemic affecting all of these countries again if enough measures and precautions are not taken. By now many potential vaccines have already completed two phases of clinical trial and will be ready for production by completing efficiency trial and getting approval [5]. Our paper completely ignores these conditions as nothing can be surely predicted. These two conditions taking active part in the pandemic the second peak can be avoided. There are very rare evidences of people getting infected twice with COVID-19 so we have not assumed any loss of immunity here [6].

As, after a few days every country has declared closure of school and college, we haven’t considered any further contribution from that matrix. The rescaling factors have been obtained in order to match the reported cases’ curve also keeping in mind the dates on which restriction have been imposed. Germany developed a very fast and novel way of testing on 26 March, which helped them to detect cases quickly [31]. Their testing rate per confirmed cases and population size is much higher than other two countries and rapid testing facility developed very early there, which caused the infected curve to go down so easily [17]. Due to this fast-acting process in Germany, we have considered \(\delta \) to increase after that date signifying faster tracing and isolation.

3 Mathematical analysis

For the wellposedness of our model, we need to check positivity and boundedness of solutions. Then we will prove existence of unique solution [13].

3.1 Positivity of solutions

Theorem 1

For \(i \in \{1,2,\ldots M\}\) if initial data of the model system (3) be \(S_{i}(0)\ge 0\), \(I_{i}(0)\ge 0\), \(Q_{i}(0)\ge 0\) and \(R_{i}(0)\ge 0\), then the solutions \(S_{i}(t), I_{i}(t), Q_{i}(t), R_{i}(t)\) of the model are positive for all time \(t>0\).

Proof

$$\begin{aligned}&\frac{dS_{i}(t)}{dt}= -\beta \lambda _{i}(t)S_{i}(t), \\&\frac{dS_{i}(t)}{dt}+\beta \lambda _{i}(t)S_{i}(t)=0. \end{aligned}$$

Multiplying both side by \(\exp \left( \int _0^t \beta \lambda _{i}(u)\mathrm {d}u\right) \) we obtain

$$\begin{aligned} \exp \left( \int _0^t \beta \lambda _{i}(u)\mathrm {d}u\right) \frac{d S_{i}(t)}{dt}+\exp \left( \int _0^t \beta \lambda _{i}(u)\mathrm {d}u\right) \beta \lambda _{i}(t)S_{i}(t)=0, \end{aligned}$$

then

$$\begin{aligned} \frac{d\left( S_{i}(t)\exp \left( \int _0^t \beta \lambda _{i}(u)\mathrm {d}u\right) \right) }{dt}=0. \end{aligned}$$

Integrating this from 0 to t, we obtain that

$$\begin{aligned} S_{i}(t)=S_{i}(0)\exp \left( -\int _0^t \beta \lambda _{i}(u)\mathrm {d}u\right) \ge 0. \end{aligned}$$

So, \(S_{i}(t)\ge 0\,\, \forall \, \, i\in \{1,2, \ldots M\}\). Note that

$$\begin{aligned} \lambda _{i}(t)S_{i}(t)&= \sum _{j=1}^{M} \frac{C_{ij}(t)I_{j}(t)S_{i}(t)}{N_{j}(t)}, \nonumber \\&= \sum _{j=1}^{M} \frac{C_{ij}(t)S_{i}(t)}{N_{j}(t)} \times I_{j}(t), \nonumber \\&= \left( \frac{C_{ii}(t)S_{i}(t)}{N_{i}(t)}\right) I_{i}(t)+\sum _{j\ne i}^{M} \frac{C_{ij}(t)S_{i}(t) I_{j}(t)}{N_{j}(t)}. \nonumber \\ \frac{dI_{i}(t)}{dt}&=\beta \lambda _{i}(t)S_{i}(t)-(\delta +\epsilon )I_{i}(t), \nonumber \\ \frac{dI_{i}(t)}{dt}&\ge \left( \beta \frac{C_{ii}(t)S_{i}(t)}{N_{i}(t)}-(\delta +\epsilon )\right) I_{i}(t). \end{aligned}$$
(4)

Similarly, we find the integrating factor be

$$\begin{aligned} \exp \left( -\left( \int _{0}^{t}\beta \frac{C_{ii}(u)S_{i}(u)}{N_{i}(u)}\mathrm {d}u\right) +(\delta +\epsilon )t\right) . \end{aligned}$$

So,

$$\begin{aligned} I_{i}(t)\ge I_{i}(0)\exp \left( \left( \int _{0}^{t}\beta \frac{C_{ii}(u)S_{i}(u)}{N_{i}(u)}\mathrm {d}u\right) -(\delta +\epsilon )t\right) \ge 0. \end{aligned}$$

So, \(I_{i}(t)\ge 0 \quad \quad \forall \,\, i\in \{1,2, \ldots M\}\).

$$\begin{aligned} \frac{dQ_{i}(t)}{dt}&= \delta I_{i}(t)-(\gamma +\mu _{i})Q_{i}(t). \end{aligned}$$

Note as \(I_{i}(t)\ge 0, \, \, \delta >0\),

$$\begin{aligned} \frac{dQ_{i}(t)}{dt}\ge -(\gamma +\mu _{i})Q_{i}(t). \end{aligned}$$
(5)

So,

$$\begin{aligned} Q_{i}(t)&\ge Q_{i}(0)\exp \left( -\int _{0}^{t}(\gamma +\mu _{i})\mathrm {d}t \right) \ge 0, \\ Q_{i}(t)&\ge 0 \quad \quad \forall \,\, i\in \{1,2,\ldots M\}, \\ \frac{dR_{i}(t)}{dt}&= \epsilon I_{i}(t)+\gamma Q_{i}(t). \end{aligned}$$

As \(I_{i}(t)\) and \(Q_{i}(t)\ge 0\) and \(\epsilon , \gamma >0\), we conclude \(R_{i}(t)\) to positive by integrating from 0 to t. \(\square \)

3.2 Boundedness of the solutions

Theorem 2

The feasible region \(\Gamma =\{(S,I,Q,R,\ldots ,Q_{M},R_{M}^{N})^T\in R_{+}^{4M}: N_{i}(t)\le N_{i}(0)\}\) is positively invarient for the model system (3).

Proof

We know that

$$\begin{aligned} N_{i}(t)&= S_{i}(t)+I_{i}(t)+Q_{i}(t)+R_{i}(t). \end{aligned}$$

After adding first four equations model system (3), we obtain,

$$\begin{aligned} \frac{dN{i}(t)}{dt}&=-\mu _i Q_{i}(t). \end{aligned}$$
(6)

From the positivity of \(S_{i}(t), I_{i}(t)\) and \(R_{i}(t)\), we conclude that \(Q_{i}(t)\le N_{i}(t)\).

Equation (6) transform to,

$$\begin{aligned} \frac{dN_{i}(t)}{dt}\ge & {} -\mu _{i}(t)N_{i}. \end{aligned}$$

So, solving we obtain, \(N_{i}(t) \ge N_{i}(0)\exp \left( -\int _{0}^{t}\mu _{i}(t)\mathrm {d}t\right) \).

Thus the \(\max \) value attained by function is \(N_{i}(0)\). So, all the solutions are bounded by some upper limit. \(\square \)

3.3 Existence of solutions

Theorem 3

The model system (3) with initial condition \(S_{i}(0)\ge 0,\, I_{i}(0)\ge 0,\, Q_{i}(0)\ge 0,\, R_{i}(0)\ge 0\) for all \(i \in \{1,2, \ldots M\}\) has a unique solution.

Proof

Let

$$\begin{aligned} X = \left( \begin{array}{cccc} {X_{1}} \\ X_{2}\\ {\vdots } \\ X_{M} \end{array}\right) \qquad \mathrm{{where}}\,\, X_{i} = \left( \begin{array}{cccc} {S_{i}(t)} \\ I_{i}(t)\\ Q_{i}(t) \\ R_{i}(t) \end{array}\right) . \end{aligned}$$

So, the model equation can be rewritten in the matrix form:

$$\begin{aligned} \varphi (x)&= AX+B(X). \end{aligned}$$

where

$$\begin{aligned} \qquad A = \left( \begin{array}{cccc} A_{1} &{}\quad 0 &{}\quad 0&{}\quad 0 \\ 0 &{}\quad A_{2} &{}\quad &{}\quad \vdots \\ 0 &{}\quad &{}\quad \ddots &{}\quad \vdots \\ 0 &{}\quad \cdots &{}\quad \cdots &{}\quad A_{M} \end{array}\right) \end{aligned}$$

with

$$\begin{aligned} A_{i} = \left( \begin{array}{cccc} 0 &{}\quad 0 &{}\quad 0 &{}\quad 0 \\ 0 &{}\quad -(\delta +\epsilon ) &{}\quad 0 &{}\quad 0\\ 0 &{}\quad \delta &{}\quad -(\mu _{i}+\gamma ) &{}\quad 0 \\ 0 &{}\quad \epsilon &{}\quad \gamma &{}\quad 0 \end{array}\right) \end{aligned}$$

Also,

$$\begin{aligned} B(X) = \left( \begin{array}{cccc} {B(X_{1})} \\ \vdots \\ {\vdots } \\ B(X_{M}) \end{array}\right) , \qquad B(X_{i}) = \left( \begin{array}{cccc} -\beta S_{i}(t) \sum _{j=1}^{M} \frac{C_{ij}(t) I_{j}(t)}{N_{j}(t)}\\ \beta S_{i}(t) \sum _{j=1}^{M} \frac{C_{ij}(t) I_{j}(t)}{N_{j}(t)}\\ 0 \\ 0 \end{array}\right) . \end{aligned}$$

Now,

$$\begin{aligned} \begin{vmatrix}B(X_{i}^{1})-B(X_{i}^{2})\end{vmatrix}= 2 \beta \begin{vmatrix}S_{i}^{1}(t) \sum _{j=1}^{M} \frac{C_{ij}(t)I_{j}^{1}(t)}{N_{j}(t)} -S_{i}^{2}(t) \sum _{j=1}^{M} \frac{C_{ij}(t)I_{j}^{2}(t)}{N_{j}(t)}\end{vmatrix}. \end{aligned}$$

From the positivity, note that \(S_{i}(t)<N_{i}(t)\) and \(I_{j}(t)<N_{j}(t)\). So,

$$\begin{aligned} \begin{vmatrix}B(X_{i}^{1})-B(X_{i}^{2}\end{vmatrix}&\le 2\beta \left( \left( \sum _{j=1}^{M} C_{ij}(t)\begin{vmatrix}S_{i}^{1}(t)-S_{i}^{2}(t) \end{vmatrix}\right) +\left( \sum _{j=1}^{M} \frac{C_{ij}(t)N_{i}(t)}{N_{j}(t)} \begin{vmatrix}I_{j}^{1}(t)-I_{j}^{2}(t)\end{vmatrix}\right) \right) , \\ \begin{vmatrix}B(X_{i}^{1})-B(X_{i}^{2})\end{vmatrix}&\le V \left( \begin{vmatrix}I_{i}^{1}(t)-I_{i}^{2}(t)\end{vmatrix}+\begin{vmatrix}S_{i}^{1}(t)-S_{i}^{2}(t)\end{vmatrix}\right) , \\ \begin{vmatrix}B(X_{i}^{1})-B(X_{i}^{2})\end{vmatrix}&\le V \begin{Vmatrix} X^{1}-X^{2} \end{Vmatrix}, \end{aligned}$$

where \(V = max \left( 2 \beta \left( \sum _{j=1}^{M} C_{ij}(t)\right) , 2 \beta \left( \sum _{j=1}^{M} \frac{C_{ij}(t)N_{i}(t)}{N_{j}(t)}\right) , \Vert A_i\Vert \right) .\)

Hence, finally we have

$$\begin{aligned} \begin{vmatrix}\varphi (X_{i}^{1})-\varphi (X_{i}^{2})\end{vmatrix}&\le V \begin{Vmatrix} X^{1}-X^{2} \end{Vmatrix}. \end{aligned}$$

From this we can say function \(\varphi \) is uniformly Lipschitz continuous and we conclude the existence of unique solution of the model system (3) [1]. \(\square \)

4 The optimal control problem

Optimal control theory was developed by Pontryagin [19]. The Pontryagin’s maximum principle allow us to lower the implementation cost and improve the efficiency functional, as a consequence of which we can optimize the performance criterion by controlling the required model parameters. Control functions for such models are mostly functions of time appearing as coefficients in the model [15]. While formulating an optimal control problem, it is very important to decide where and how to introduce the control in the system of differential equations.

Now we include control \(u_{1i}\) and \(u_{2i}\) in the model (3) which represents the medical care of infected patients with COVID-19 and state efforts to encourage people to sanitize and use mask. Our main objective in this proposed control strategy is to minimize the number of infected and quarantined individuals. Thus, the controlled mathematical system is given by the following system of differential equations:

$$\begin{aligned} {\dot{S}}_i(t)= & {} -\beta \lambda _i S_i (1-u_{1i})=f_1, \nonumber \\ {\dot{I}}_i(t)= & {} \beta \lambda _i S_i (1-u_{1i})\,-\, ({\delta } \,+\, {\epsilon }+u_{2i})I_i=f_2, \nonumber \\ {\dot{Q}}_i(t)= & {} {\delta }I_i \,-\, (\gamma \,+\, \mu _i+u_{2i})Q_i=f_3, \nonumber \\ {\dot{R}}_i(t)= & {} \epsilon I_i \,+\, \gamma Q_i=f_4, \nonumber \\ {\dot{F}}_i(t)= & {} \mu _i Q_i=f_5, \end{aligned}$$
(7)

where \(S_i(0) \ge 0, I_i(0) \ge 0, Q_i(0) \ge 0, R_i(0) \ge 0, F_i(0) \ge 0,\) for all \(i \in {1,2,...M}\) are given initial states. Next, we define our objective function using bounded measurable control as

$$\begin{aligned} {\mathcal {J}}(u_{1i},u_{2i})=\int _0^T\left( I_i(t)+Q_i(t)+\dfrac{B_1}{2} u_{1i}^2+\dfrac{B_2}{2} u_{2i}^2\right) dt, \end{aligned}$$
(8)

subject to the state system (7). In order to find an optimal solution, we first find the Lagrangian and Hamiltonian for the optimal control problem (7)–(8). The Langrangian of the control problem is given by

$$\begin{aligned} L=I_i(t)+Q_i(t)+\dfrac{B_1}{2} u_{1i}^2+\dfrac{B_2}{2} u_{2i}^2, \end{aligned}$$
(9)

where the coefficients \(B_1, B_2\) represents the balancing cost factors for control strategies \(u_{1i}, u_{2i}\) respectively.

The optimal control problem involves in finding \(u_{1i}^*,u_{2i}^*\) such that the associated state trajectories \((S^*, I^*, Q^*, R^*, F^*)\) is solution of the controlled system of equation in the interval [0, T] with the initial conditions and minimizing the cost functional \({\mathcal {J}}\).

$$\begin{aligned} {\mathcal {J}}(u_{1i}, u_{2i})=min_{u_{1i},u_{2i} \in \Delta } \{{\mathcal {J}}(u_{1i},u_{2i})|~u_{1i},~u_{2i} \in \Delta \} \end{aligned}$$
(10)

where \(\Delta \) is the set of admissible controls given by

$$\begin{aligned} \Delta =\{(u_{1i},u_{2i})\in L^1[0,T]| 0\le u_{1i}\le 1, 0\le u_{2i}\le 1\}, t\in [0,T]. \end{aligned}$$
(11)

Pontryagin’s Maximum Principle converts (7), (8) and (9) into a problem of minimizing a Hamiltonian function (H), defined as

$$\begin{aligned} H(t)= I(t)+Q(t)+\sum _{k=1}^5 \delta _i(t) f_k(S^j,I^j, Q^j, R^j, F^j), \end{aligned}$$
(12)

where \(f_k\) are given in (7).

Theorem 4

Given the optimal control \(u_{1i}, u_{2i}\) and the solutions \(S_i, I_i, Q_i, R_i, F_i\) of the corresponding model system (3) there exists adjoint variables \(\delta _1, \delta _2, \delta _3, \delta _4, \delta _5 \) satisfying

$$\begin{aligned} \dot{\delta _1}(t)= & {} -(u_{1i}-1) \beta (\delta _{1}-I_i \delta _{2}) \lambda _{i}, \nonumber \\ \dot{\delta _2}(t)= & {} -1- u_{2i} \delta _{2}+ \delta (\delta _{2}-\delta _{3}) - (\delta _{2}+\delta _{4}) \epsilon +S_{i}(t) (-1+u_{1i}) \beta \delta _{2} \lambda _{i},\nonumber \\ \dot{\delta _3}(t)= & {} - 1+u_{2i} \delta _{3}+\gamma (\delta _{3}-\delta _{4}) +\delta _{3} \mu _{i}-\delta _{5} \mu _{i}, \nonumber \\ \dot{\delta _4}(t)= & {} 0, \nonumber \\ \dot{\delta _5}(t)= & {} 0, \end{aligned}$$
(13)

with the tranversality conditions at time T: \(\delta _1=0\); \(\delta _2=0\); \(\delta _3=0; \delta _4= 0; \delta _5= 0.\) Furthermore, for \(t\in [0,T]\), the optimal controls \(u_{1i}\) and \(u_{2i}\) are given by:

$$\begin{aligned} u_{1i}^*= & {} min\left( u_{1i}^{max},max \left( u_{1j}^{min},\frac{S_{i} \beta (I_i \delta _{2}-\delta _{1}) \lambda _{1}}{B_{1}}\right) \right) , \end{aligned}$$
(14)
$$\begin{aligned} u_{2i}^*= & {} min\left( u_{2i}^{max},max \left( u_{2j}^{min},\frac{Q_{i} \delta _{3}-I_i \delta _{2}}{B_{2}}\right) \right) . \end{aligned}$$
(15)

Proof

The Hamiltonian at time t is given by

$$\begin{aligned} H(t)= I(t)+Q(t)+\sum _{k=1}^5 \delta _i(t) f_k(S^j,I^j, Q^j, R^j, F^j), \end{aligned}$$

where

$$\begin{aligned} \begin{aligned} f_1 \;&=\; -\beta \lambda _i S_i (1-u_{1i}), \\ f_2 \;&=\; \beta \lambda _i S_i (1-u_{1i})\,-\, ({\delta } \,+\, {\epsilon }+u_{2i})I_i, \\ f_3 \;&=\; {\delta }I_i \,-\, (\gamma \,+\, \mu _i+u_{2i})Q_i, \\ f_4 \;&=\; \epsilon I_i \,+\, \gamma Q_i, \\ f_5 \;&=\; \mu _i Q_i, \end{aligned} \end{aligned}$$
(16)

for \(t\in [0,T]\), the adjoint equations and transversality conditions can be obtained by Pointryagin’s maximum principle, such that

$$\begin{aligned} \dot{\delta _1}= & {} -(u_{1i}-1) \beta (\delta _{1}-I_i \delta _{2}) \lambda _{i}, \nonumber \\ \dot{\delta _2}= & {} -1- u_{2i} \delta _{2}+ \delta (\delta _{2}-\delta _{3}) - (\delta _{2}+\delta _{4}) \epsilon +S_{i} (-1+u_{1i}) \beta \delta _{2} \lambda _{i}, \nonumber \\ \dot{\delta _3}= & {} - 1+u_{2i} \delta _{3}+\gamma (\delta _{3}-\delta _{4}) +\delta _{3} \mu _{i}-\delta _{5} \mu _{i}, \nonumber \\ \dot{\delta _4}= & {} 0, \nonumber \\ \dot{\delta _5}= & {} 0, \end{aligned}$$
(17)

with the tranversality conditions at time T: \(\delta _1=0\); \(\delta _2=0\); \(\delta _3=0; \delta _4= 0; \delta _5= 0\). Furthermore, for \(t\in [0,T]\), the optimal controls \(u_{1i}^*\) and \(u_{2i}^*\) can be solved from the optimality condition: \(\dfrac{dH}{du_{1i}}=0\) and \(\dfrac{dH}{du_{2i}}=0\), that is

$$\begin{aligned} \dfrac{dH}{du_{1i}}= & {} B_1 u_{1i}+ S_{i} \beta (\delta _{1}-I_i \delta _2) \lambda _{i}=0,\\ \dfrac{dH}{du_{2i}}= & {} B_{2} u_{2i}+I_i \delta _{2}-Q_{i} \delta _{3}=0, \end{aligned}$$

so we have \(u_{1i}^*\), \(u_{2i}^*\) as given in (14)–(15). This completes the proof. \(\square \)

5 Data and methodology

For each country, 5M differential equations were set up and numerically integrated using odeint from the ‘numpy’ package in python. The case and fatality data was obtained from Worldometers [34], population data from Population Pyramid [20], and contact matrix data from Prem et al. [21]. Data-set for each country is structured as day number, date, deaths, active cases from left to right column-wise. The website presents various type of data (active cases, total confirmed cases, deaths, daily new confirmed cases, daily new deaths) in graphs for the countries collected from the respective governmental sites and other authenticated sources.

The values of the parameters \({\bar{\delta }}\) (rate of infection), \({\bar{\epsilon }}\) (rate of recovery) and \(\gamma \) (appearance of symptoms) have been set to pre-determined constant [16, 28]. Though several other study suggests slightly different value that are well in the confidence range. For the purposes of simulation, we assume \(a = 1/2\), which is consistent with studies such as [9, 36]. Initial conditions set as \(I_i(0) = Q_i(0)\) on the assumption that positive cases are under reported such that only half of all infected individuals have been identified and quarantined. Depending upon the restrictions imposed, we initially started with a reasonable value of \(\alpha \)’s. Then, we optimized the curve with respect to coefficient \(\beta \) first using the scipy optimizer as this parameter is generally responsible for the trend of the graph. After that, we manually tried to modify \(\alpha \)’s by small amounts to closely match the curve. By fitting numerically to existing data the parameters \(\mu \) and \(\beta \) are obtained, as seen in Fig. 2.

Fig. 2
figure 2

Simulation data fitted to case data in (a) France, (b) Spain and (c) Germany. Note that the ‘infected’, ‘quarantined’, and ‘total infected’ curves are simulated, while the ‘actual cases’ data-points have been obtained from Worldometers [34]

Full size image

The choice of transmission coefficient \(\beta \) being same for all age groups is justified due to computational complexity of fitting 16 parameters as well as the absence of proper estimation of age-group specific transmission probabilities. It should be noted that usage of single \(\beta \) does not affect much the result as the main disease dynamics is encoded in contact matrix. For different age group, transmission factor is different which is multiplication of transmission probability (\(\beta \)) and contact matrix (\(\ C_{ij}\)). But, disease induced mortality rates across all age groups has been fitted to be the different values, \(\mu \) in Table 4. While these rates may indeed differ from time to time in reality due to change in the nature of virus and availability of specified medical facility in later stage, which justifies this assumption with this model closely imitating the trends in reality without additional parameters and complexity being introduced. Furthermore, due to lack of age dependent immunity data and its implementation in our model, and the coefficients \(\lambda _i\) denotes normalized contacts between age groups. Our initial assumption of half of the cases being reported, can be better set by getting proper estimation of the fraction of cases reported.

6 Results and discussion

As a response to the pandemic, every country has implemented physical distancing measures, along with travel bans and nationwide or partial lockdowns. On 13 March, 2020, the Prime Minister of France ordered the closure of all non-essential public places. The French President issued a mandatory notice of staying at home for the citizens on 16 March, 2020 [30]. With an increasing number of cases in Spain, the Prime Minister declared “a state of alarm”, closing all non-essential activities and putting citizens under a ‘stay at home’ order, only exempting essential work [32]. On 26 February, 2020, Germany started the closure of schools, libraries and swimming pools, and canceled important events throughout the country. State authorities declared further closures day by day. Bans on gathering and travel, and the mandatory wearing of masks started to come in the picture from 8 March onwards, though not many of these restrictions were imposed by the Federal Government [31]. Spain and France have faced criticism for their late reactions to the disease. On the other hand, Germany was quick enough to put down restrictions and contain the disease very well. Some countries in Europe like Iceland, Croatia have successfully tackled the disease, while other countries like UK, Spain allowed sports activity, International Women’s Day celebration to go on allowing gatherings, which proved to be spreading centres in later stages [33]. Thus, taking early measures would have helped these countries to make their way out of the global pandemic with fewer cases.

Fig. 3
figure 3

Predicted infected and quarantined populations in (a) France, (b) Spain and (c) Germany. The dashed curves in (a) and (b) represent the consequences of a delayed relaxation of restrictions by a month (For France May 22 to June 21 and for Spain May 4 to June 3). The vertical green line denotes the date 25 July, 2020

Full size image

The infection curves generated by this model are displayed in Fig. 3. It’s important to note that a fast rising infection curve has been predicted, but which takes much longer time to go down for France and Spain. The peak of reported cases(interpreted as quarantined individuals in our model) trails behind the peak of the total infected curve. In the case of France, we observed a long enough flat top. Spain and Germany have already crossed the peak. The disease will be soon eradicated in Germany, but as people are going back to activities, there is also a little bit chance of a second wave there if they drop further restrictions. Spain and France will continue to see upsurge in number of cases in the absence of any herd immunity or vaccines paving the way for second peak, which has been caused by their implementing restart policies at wrong time.

As previously discussed, it is important to impose lockdown and restriction at proper time so is its withdrawal. Ideally a country should put complete lockdown until its active cases go down to zero, but as lockdown comes with its devastating effects on people’s income and country’s economy in a realistic scenario it is not possible. As countries started to observe lowering trend of new affected patients they started to figure out de-escalation plans comprising of several stages and zone-wise restriction lifting. In Spain sports activities with no spectators, market and worship places with less capacity came into action from 11 May, 2020 in some areas [32]. In France several lockdown lifting phases came with different restrictions uplifted unlike Spain. They also started their first phase on 11 May, 2020 allowing meeting, gatherings, activities and travelling with basic social distancing rules [30]. In case of Germany state governments started to lift all restrictions from 9 June. On June 3rd, they allowed people to travel to EU countries. From the dates, we can understand measures taken in Germany were for a longer period of time than other two countries, which is the main reason of its success [31]. France and Spain tried to make life normal for people when there were still around sixty thousands active cases, where Germany initiated lifting of restrictions when there were only around eight thousands active cases.

Re-start in France and Spain with so much active cases made us curious to further simulate the effect what if their re-opening policies are delayed by a month (For France May 22 to June 21 and for Spain May 4 to June 3), as it would be the same time when Germany lifted its restrictions. We observe by this delayed relaxation in France and Spain number of cases would be significantly less by now, which can be clearly seen from dotted line in Fig. 3. Though it would not help in completely mitigating the disease like Germany as these countries had acted late during the first phase. In the hypothetical scenario, infected curves for both the countries are looking to rise with lowered peak and significant delay than the current scenario, which can helpful in developing herd immunity or availability of vaccine.

Fig. 4
figure 4

Relative decrease in susceptible populations in every age-group in (a) France, (b) Spain and (c) Germany

Full size image

In our model, due to usage of single \(\beta \) value for different age-groups, we need to verify the results from the actual scenario. As infection is less in children, we have not considered 0–15 age-groups in the category of mostly infected age-groups [35]. There is a common trend among the 30–35 age-group being the most affected age-group in Spain and France where 25–30 is most affected for Germany. Moreover, we have seen 15–40 age-groups are the mostly infected age-group for all. We estimate this by measuring the relative drop in \(S_i\) from it’s initial value \(S_{i,0}\) across every age groups over time. This has been illustrated in Fig. 4.

Fig. 5
figure 5

The variation of the effective reproductive number \({\mathfrak {R}}_0^\text {eff}\) with time in (a) France, (b) Spain and (c) Germany

Full size image
Table 2 \({\mathfrak {R}}_0^\text {eff}\) at different time periods
Full size table

The basic reproductive ratio \({\mathfrak {R}}_0\) can be thought as the new cases directly caused by a single infected person. Growth or dying out of a disease essentially depends upon this number. The basic reproductive ratio \({\mathfrak {R}}_0\) is calculated numerically using the next generation method [27]. The detailed calculation gave us that the basic reproductive number \({\mathfrak {R}}_0= \rho ({\mathcal {F}}{\mathcal {V}}^{-1})\) is simply calculated as the spectral radius of the matrix \({\mathcal {F}}{\mathcal {V}}^{-1}\) [26]. Note that the matrix \([C_{ij}N_i/N_j]_{M\times M}\) has strictly positive entries, and thus has a unique, positive, real largest eigenvalue r. Thus, we have \({\mathfrak {R}}_0= \beta r/(\delta + \epsilon )\).

\({\mathfrak {R}}_0^\text {eff}(t)\) (effective reproduction number) is the number of people in a population who an individual can infect at any specific time. It changes as the population becomes increasingly immunized, either by individual immunity following infection or as people die. Factors affecting \({\mathfrak {R}}_0^\text {eff}(t)\) include the number of people with infection, the number of susceptible people infected, and people’s behavior, such as social distancing. Now, we develop \({\mathfrak {R}}_0^\text {eff}(t)\) as the spectral radius of \(\bar{{\mathcal {F}}} \bar{{\mathcal {V}}}^{-1}\). \(\bar{{\mathcal {F}}}\) and \(\bar{{\mathcal {V}}}\) denotes effective \({\mathcal {F}}\) and \({\mathcal {V}}\) matrices. This ratio denotes the instantaneous, effective growth of the infected compartments at a given time. In the disease-free equilibrium, stability of the system depends upon value of \({\mathfrak {R}}_0\). When \({\mathfrak {R}}_0\;<\; 1,\) the system is asymptotically stable, but unstable if \({\mathfrak {R}}_0\;>\; 1\) [27]. We have presented the results in Fig. 5. Note that \({\mathfrak {R}}_0\) depends on (i) transmission coefficient \(\beta \), (ii) social contact coefficients \(C_{ij}\), (iii) the fraction of asymptomatic people a, (iv) rate of appearance of symptoms \(\delta \) and (v) rate of recovery of asymptomatic individuals \(\epsilon \). In the initial stages, we observed values between 3 and 4 for the countries. We noted the highest value of \({\mathfrak {R}}_0\) for Spain, followed by France and lowest in Germany among the three. In the midway, basic reproductive ratio of Spain, France and Germany had gone below 1, suggesting the dying out of the disease. But restart has empowered the ratio to go just little above 1 in France and Spain, which is a matter of concern for these countries as we are observing rise in infected people again. As of now the \({\mathfrak {R}}_0\) value in France and Spain, which is higher than 1, suggests spreading of the disease and the infected curve to rise. The values of \({\mathfrak {R}}_0^\text {eff}\) at these times of interest have been presented in Table 2.

On 13 March, 2020, WHO formally declared the outbreak as a pandemic, which acted as an alarm for most of the countries. As the disease started earlier in some parts of Europe, taking immediate measures could not help them control the disease. Countries like Iceland though having less number of cases put ban on public gatherings and closed school, colleges from 16 March, which proved to be fruitful for them. On the other hand countries like Croatia got enough time to act after WHO’s declaration and as of then there were less number of cases [29]. If Spain and France had acted early, they might be able to delay, but no reduction, which suggests they should increase the measures. Though total confirmed cases of Germany seems to be higher than France but looking at their population size (c.f. Table 3) it is crystal clear that Germany had done a far better job compared to France and Spain.

7 Conclusions

There is a huge chance of missing data affecting model parameters during an ongoing pandemic, leading to wrong predictions. We have also made some assumptions like no aging in population; there is no birth and normal death (except death due to disease). We have focused on people less than eighty years old as the remaining population is very small compared to the below 80 residents. In addition, this model does not consider the 75-80 year age group in Spain due to the unavailability of contact data. Keeping in mind these assumptions, our model is good enough for short-term predictions, but there may be some errors in the forecast as we have tried to simulate for nearly ten months. Moreover, partial lockdown of some hotspot (in terms of more infected) areas is not too easy to model because the disease may not spread from these areas, but in simulation, these effects can’t be shown as we have taken the country as a whole. In our model, we have continued the measure, as it is happening now. Though there is specific news of relaxing the measures in the future, we identify that maintaining these measures in highly infected areas and lifting some measures in other areas can affect in the same way as maintaining as a whole.

Some scientists/researchers tried to survey for collecting the contact pattern of age-groups [11]. But as they have rightly mentioned about the recall and selection bias, it can lead to false prediction. We have used rescaling factors for the POLYMOD contact patterns data, which have been considered as the effect of physical distancing. The disease dynamics are essentially encrypted in the contact matrix data. So collecting the information about the first observed case and restriction policies implemented in that country, this model can be used further in different parts of the world to interpret measures taken. Our model gives the scope of designing age-group-specific policies to battle the pandemic situation. To understand the importance of age-targeted interventions of epidemic contention, we fit our model with the data across three countries (i) Spain, (ii) France, (iii) Germany. Our work found that the mortality rate is maximum for the age group 65–75 for all the countries (c.f. Table 4). To minimize the death in older age groups, we need to focus on the old first and provide vaccination to them on a priority basis. However, the spread of the disease is mainly propagated by the people of the 30–35 age-group, so vaccinating them can essentially stop the spread of the disease. So, upon the availability of vaccines, it’s a tradeoff between the number of affected individuals and the number of deaths registered. With a proper abundance of data (age-group wise active cases), this model can also be used to determine age-structured transmission factors that will help to understand which age-groups are more prone to develop and transmit the disease. Policies like restrictions for particular age-groups or age-targeted rapid testing and isolating them to identify super-spreaders can be adopted with the help of this model.

Table 3 Age structured population data, used as initial values for \(S_i\)
Full size table
Table 4 Age structured mortality rate \(\mu _i\)
Full size table

We note that implementing lockdown has a tremendous effect on the country’s economy, so in the ideal case scenario, keeping the full restrictions until active cases go to zero is impossible in real condition. Lifting restrictions with proper understating of the future dynamics of the disease is very important in which part our model can help design the unlock plans. As all models try to approximate real scenarios, this age-structured model can help develop policy with the social, medical, and economic perspectives.

Future work with the age-structured model can be done by changing it to a continuum model and inclusion of birth and natural death with aging. This model can incorporate additional features for a greater realistic scenario with a true abundance of data. It will be nice to add non-local diffusion to our model to examine the spatio-temporal pattern of this pandemic.