# Deciphering dynamics of recent COVID-19 outbreak in India: An age-structured
modeling

Vijay Pal Bajiya1, Jai Prakash Tripathi1,Ranjit Kumar Upadhyay2

1Department of Mathematics, Central University of Rajasthan,
Kishangarh-305817, Ajmer, India
2 Department of Mathematics & Computing, Indian Institute of Technology (ISM),
Dhanbad-826004, India

###### Abstract

The transmission dynamics of an infectious disease are most sensitive to the
social contact patterns in a population of a particular community and to
analyze the precautions people use to reduce the transmission of the disease.
The social contact pattern depends on the age distribution of the specific
community. Therefore, knowing the age-specific prevalence and incidence of the
infectious disease is essential for modeling the future burden of the disease
and the effectiveness of interventions such as vaccination. In the present
study, we consider an SEIR age-structured multi-group epidemic model to
understand the impact of social contact patterns in controlling the disease.
To observe how fluctuations in social mixing have affected the spreading of
the emerging infectious disease, we used synthetic location-specific social
contact matrices in the community. For mathematical analysis, we computed the
basic reproduction number $(R_{0})$ for the system and also illustrated the
global behavior of the system in terms of basic reproduction number. Further,
the existence of optimal control for the associated problem has been
established and computed mathematically. The transmission rate for the
proposed model using the real data of COVID-19 for India from September 1,
2020, to December 31, 2020, has also been estimated. We simulated lifting of
different non-pharmaceutical interventions by permitting the people to go back
for their works in a phased-manners and investigated the effects of returning
to work at different stages, accordingly. Our results suggest that awareness
of symptomatic infected individuals of age groups $20-49$ years is beneficial
to reduce the number of infected individuals when all schools are closed.
However, awareness of symptomatic infected individuals of school children age
groups also plays a significant role in reducing disease cases when some
schools are partially opened. The simulation results also recommend that the
number of cases could be reduced in large numbers by controlling the contacts
at school and other gathering places. Interestingly, it has been investigated
that the time-dependent transmission rate is more realistic rather than the
constant spread rate to COVID-19 for India via estimating transmission rate
using the least square method. Our study suggests that the early and sudden
lifting of control measures could lead to other peaks and a high COVID-19
burden, which could be flattened and reduced by relaxing the interventions
gradually. We hope that our results would help health policymakers in deciding
appropriate and timely age-based vaccination distribution strategies and,
therefore, control the disease.

_Keywords_ : Social contact matrix, Age-structured epidemic model, Basic
reproduction number, Optimal control, COVID-19, Parameter estimation.

##  1 Introduction

Mathematical modeling of infectious disease spreading has become an essential
tool for understanding disease dynamics and outbreak patterns. It also
provides crucial insights to the policymakers in making timely decisions to
control and reduce the burden of the diseases when limited empirical data are
available. The accomplishment of mathematical modeling in advising critical
decisions to human shield has been confirmed for various diseases, including
pandemic Influenza, COVID-19, Flu, etc. [1, 2, 3, 4]. Since the infection
ability of an infected individual to transmission rapidly and disturb many
people in a community. Therefore, infectious diseases ( for example,
Influenza, COVID-19, etc.) directly communicated from individual to individual
by the respiratory system have been of particular interest for mathematical
modeling. The importance of epidemic models and the usefulness of policies
based on these models are dependent on the robustness of the model parameters,
which capture different features of the disease [1, 5]. The key parameter in
epidemic modeling is the probability of successful contact between an
infectious and a susceptible individual to spread the infection. The accurate
contact structures in the epidemic modeling improve the accuracy and
efficiency of the prediction and allow us to investigate the effect of control
measures targeted at specific locations, such as schools, workplaces, or
homes. Therefore, social contact structures of people are the important
factors in epidemic modeling to understand the transmission dynamics and
investigate the control measures of infectious diseases [6, 7].

Several assumptions are mandatory to make the range of human relations
straightforward into controllable mathematical models of communicable diseases
transmitted from one person to another person. The essential assumption of a
homogeneous mixing population (in which each individual has an equal
probability of contact to infected person) has been transformed by the
different realistic frameworks, in which the likelihood of contact fluctuates
between different groups, which may be defined by the age of individuals. The
extent to which individuals specially mix with people of the same age
(assortativeness mixing) is a crucial heterogeneity that is now routinely
included in epidemic modeling. The heterogeneity has also been made to further
represent the fundamental structure of social contact patterns by breakdown
the population according to their different locations such as households,
schools, workplaces, and other places [6, 8]. For the directly transmitted
respiratory virus (COVID-19, Measles, and Influenza), social mixing patterns
influence the risk of individual-level transmission of diseases [10, 11] and
population-level infection dynamics [12], as well as the effectiveness of age-
specific control measures [9]. The government of India declared a state of
emergency and announced the nationwide lockdown and shut down schools. Sports
facilities, non-food shops, restaurants, shopping malls, and traveling
prohibited from March 25, 2020 (Phase-1 lockdown) to reduce social contacts.
The authority has also instructed to follow the strict distancing measures,
avoid unnecessary social interactions, and the local people to stay at home
except for essential works. Social gathering was also prohibited, and
therefore, social contacts were also reduced significantly. During the
lockdown lifting, social activities such as market opening and traveling with
social distancing and face-masks were allowed step by step. The social
gathering was also permitted with the limited peoples (for example, marriage
functions with fifty people gathering are allowed). Thus, social contact
patterns essentially play a significant role in understanding the disease
progression.

The changes in social mixing patterns among people of various age groups and
age distribution affect the prevalence of each age group. The child age-group
individuals are subject to make more unnecessary social contacts than adults
due to their unawareness [15] and hence, children may pay more contribution to
disease transmission than adults [16, 17]. In the early stage of COVID-19 in
India, $75\%$ of confirmed cases belong to the actively working population
between the age group of 21 to 60 years [18]. The primary reason for this is
that the maximum number of persons in those age groups travel internationally
for their jobs and businesses. They were also working in essential services
during the lockdown period in India. Therefore, the number of cases or
prevalence of the disease depends strongly on the role of different age groups
and the mixing patterns. Other age distributions could show significantly
different epidemic shapes and the overall impact on the transmission dynamics
of infectious diseases. School closures are considered an essential
intervention for an epidemic due to the higher rate of unnecessary social
contact in children. Thus, the impact of school closure depends on children’s
contribution to disease transmission. The early cases of SARS-CoV-2 in Wuhan,
China, were concentrated in adults over 40 years of age due to skewed age
distribution in China. The assortative mixing patterns between adults could
have condensed infection transmission to children in China’s very early stages
of the COVID-19 outbreak. In many countries (outside China), COVID-19
outbreaks may have been initiated by working-age travelers entering the
country [14] producing a similar excess of adults in the early phases of local
epidemics. In both cases, the school closures that happened potentially
further declined the transmission among children, but to what degree is
unclear.

The heterogeneity of the infected individuals and cruelty according to their
age groups, especially for children and older people, encouraged the interest
of several researchers [19, 20, 21, 22]. Some studies have exposed that the
cruelty of infectious diseases increases with the age and morbidity of
hospitalized patients [23, 24, 25]. Wu et al. [22] have estimated that the
global symptomatic case fatality risk of Coronavirus in Wuhan was $1.4\%$,
which is significantly lower than the corresponding crude confirmed case
fatality risk of $4.5\%.$ Their findings also suggest that the risk of showing
symptoms rises by $4\%$ per year in adults who have age groups between $30$
and $60$ years. Davies et al. [26] considered an age-structured mathematical
model to examine the epidemic data of COVID-19 from China, Italy, Japan,
Singapore, Canada, and South Korea. Their study recommends that there is a
robust connection between the age of an individual and the probability of
showing symptoms. Since asymptomatic individuals can also be infectious, a
higher possibility of showing symptoms does not indeed have greater
infectiousness. The study of Zou et al. [27] has shown that the viral load in
the asymptomatic cases was very similar to that in the symptomatic cases in
many situations. Jones et al. [28] established that viral loads in infected
children with age group $0-15$ years do not fluctuate significantly from those
of older age people. However, older age people are more likely to develop
symptoms due to their immunity.

The existing works of the literature suggest that examination of the
transmission dynamics between two generations is essential to a better
understanding of COVID-19 transmission and most fundamental to studying the
mitigation interventions of COVID-19 efficiently. Therefore, social contact
patterns among age groups in the community (like work, school, home, and other
locations) are a fundamental factor to consider while modeling the
transmission of the COVID-19 pandemic. To incorporate these factors, various
mathematical models have been established [20, 21, 26, 29, 30]. Chikina and
Pegden [30] studied an age-structured model to explore age-targeted mitigation
interventions. They use the numerical values of the parameter from the
literature and make discussion using the age-structured temporal series to fit
their considered model. However, Davies et al. [26] also illustrate age-
related effects in controlling the COVID-19 pandemic. The authors also discuss
the efficiency of different control measures by using statistical inference to
fit an age-structured SIR model output to empirical data.

In the present study, we use synthetic location-specific contact patterns
(contact matrices) for India to observe how these fluctuations of social
mixing among different age groups have affected the disease outbreak
progression. We modify these in the case of school and workplace closures,
reduction in mixing of the population in the general community, and some other
intervention scenarios. Using these matrices and the estimated epidemiological
parameters for the COVID-19 outbreak in India [26, 51, 53], we estimate the
infection rate $\beta$ according to lockdown lifting scenarios and simulates
the trajectory of COVID-19 outbreak in India through an age-structured (a
multi-group) $SEI^{a}I^{s}R$ epidemic model for various types of intervention
scenarios. The main findings of present work are (i) awareness of the
symptomatic infected individuals of age groups $20-49$ years has been an
important control measure to COVID-19 in India (ii) the number of cases could
be reduced in large number by controlling the contacts at school and other
gathering places (iii) the time-dependent transmission rate may be more
realistic rather than the constant rate to COVID-19 for India.

The rest of the paper is formulated as follows. In Section 2, we formulate an
$SEI^{a}I^{s}R$ multi-group epidemic model incorporating the social mixing
pattern to quantify the control measures of the disease. The computation of
basic reproduction number $(R_{0})$, local and global stability of equilibria
in terms of $R_{0}$ are discussed in Section 3. We also formulated and
computed the optimal control problem in Section 4 and Existence of an optimal
control. In Section 5, we presented the numerical simulations to the case
study of COVID-19 in India, including the model parameter estimation. We also
discussed the effects of awareness of the symptomatic infected individuals and
various scenarios of social mixing patterns. Finally, we discuss the findings
and conclude our work in Section 6.

##  2 Model Formulation

In the case of a homogeneous mixing population, we can divide the total
population $(N)$ into five main epidemiological classes, susceptible $(S)$,
exposed $(E)$, asymptomatic infected $(I^{a})$, symptomatic infected $(I^{s})$
and recovered $(R)$. Thus, we can achieve the following system of ordinary
differential equations to describe the infectious disease dynamics.

$\displaystyle\frac{dS(t)}{dt}=$ $\displaystyle\Gamma-\beta
S(t)\lambda_{1}\left(\frac{I^{a}(t)+\eta I^{s}(t)}{N(t)}\right)-\mu S,$ (1)
$\displaystyle\frac{dE(t)}{dt}=$ $\displaystyle\beta
S(t)\lambda_{1}\left(\frac{I^{a}(t)+\eta
I^{s}(t)}{N(t)}\right)-(\alpha+\mu)E(t),$ $\displaystyle\frac{dI^{a}(t)}{dt}=$
$\displaystyle(1-\sigma)\alpha
E(t)-\left(\gamma^{a}+\delta^{a}+\mu\right)I^{a}(t),$
$\displaystyle\frac{dI^{s}(t)}{dt}=$ $\displaystyle\sigma\alpha
E(t)-\left(\gamma^{s}+\delta^{s}+\mu\right)I^{s}(t),$
$\displaystyle\frac{dR(t)}{dt}=$
$\displaystyle\gamma^{a}I^{a}(t)+\gamma^{s}I^{s}(t)-\mu R(t),$

where $N(t)=S(t)+E(t)+I^{a}(t)+I^{s}(t)+R(t).$ $\Gamma$ represents the new
recruitment in the susceptible population and $\mu$ be the natural death rate
of individuals. $\alpha$ is the rate at which the exposed are fetching
infectious individuals, and $\sigma$ represents the proportion of symptomatic
infectious in total infectious individuals. $\delta^{a}$ and $\delta^{s}$ are
the mortality rates due to the disease of asymptomatic and symptomatic
infectious individuals, respectively. $\gamma^{a}$ and $\gamma^{s}$ are the
recovery rates of asymptomatic and symptomatic infectious individuals,
respectively. Assuming that the symptomatic infectious individuals avoid
unnecessary contact with others (by self-isolation and other precautions), the
transmission of disease by the symptomatic infectious individuals is reduced
than asymptomatic infectious individuals. Therefore, $\eta$ represents the
proportion of contacts avoided by symptomatic infectious individuals due to
their awareness. Thus, $\eta$ also represents the awareness of symptomatic
infectious individuals. $\lambda_{1}$ be the average number of contact between
the susceptible and infectious individuals, and represents the probability of
infection on the contact between infectious and susceptible individuals.
During the initial phases of the infectious disease, disease spreading between
individuals is statistically independent, which means the probability of
making contact between an infectious individual and someone no longer
susceptible is very low. in, Epidemiology, the basic reproduction number
$R_{0}$ is used to forecast the curve of an epidemic such that disease will
eliminate when $R_{0}<1$ and persist when $R_{0}>1$. $R_{0}$ represents the
number of people that an infectious is expected to infect and can be
calculated by next-generation matrix approach [31]. The basic reproduction
number $(R_{0})$ for the model system (1) is given by.

$R_{0}=\frac{\beta\lambda_{1}\alpha(1-\sigma)}{(\alpha+\sigma)(\gamma^{a}+\delta^{a}+\mu)}+\frac{\beta\lambda_{1}\alpha\eta\sigma}{(\alpha+\mu)(\gamma^{s}+\delta^{s}+\mu)}$
(2)

where
$\frac{\beta\lambda_{1}\alpha(1-\sigma)}{(\alpha+\sigma)(\gamma^{a}+\delta^{a}+\mu)}$
and
$\frac{\beta\lambda_{1}\alpha\eta\sigma}{(\alpha+\mu)(\gamma^{s}+\delta^{s}+\mu)}$
represent the contribution of asymptomatic and symptomatic infectious
individuals in the disease transmission, respectively.
In general, it must be noted that the value of the basic reproduction number
is not a biological constant. However, it depends on various epidemiological
aspects like contact patterns among individuals, applied interventions. Also,
the number $R_{0}$ of the disease usually depends on newly susceptible
populations. In particular, a value of $R_{0}>1$ specifies that the infectious
disease will initiate to spread in the population if there does not exist any
intervention. However, a higher value of $R_{0}$ represents the faster
exponential growth of infection in a community. For example, measles is known
to be one of the most contagious diseases, with $12\leq R_{0}\leq 18$ [32].
The CDC determined that COVID-19 has an $R_{0}$ approximately value of $5.7$
for the United States [33]. It is near to that of Polio and Rubella [34]. For
many infectious diseases such as COVID-19, the mortality rate due to disease
does not distribute homogeneously. However, it varies according to the age
groups of individuals [26].
For many diseases, such as COVID-19, the impact of the different age groups on
the disease transmission dynamics varies drastically. It may happen due to
different contact patterns among age groups and different mixing patterns.
Therefore, we formulate the transmission dynamics model for an outbreak of the
disease in the heterogeneously mixing population to quantify the impact of the
contact patterns and age of individuals. We consider an age-structured
epidemic model (multi-group epidemic model) in which each group of a
compartment represents the dynamics of individuals of the age groups like
$0-04,05-09,10-14,15-19,\cdots,75-79.$ That means
$S_{i},E_{i},I^{a}_{i},I^{s}_{i}$ , and $R_{i}$ represent the number of
susceptible, exposed, asymptomatic, symptomatic, and recovered individuals of
$i^{th}$ age group. In this age-structured epidemic model, we also incorporate
the contact patterns according to their age groups. It is described by the
social contact matrix, which contains the average number of contact between
each pair of age groups. Therefore, the force of infection for age-structured
epidemic model takes the following form:

$\text{force of infection}=\beta
S_{i}(t)\sum_{j=1}^{n}\lambda_{ij}\left(\frac{I^{a}_{j}(t)+\eta_{j}I^{s}_{j}(t)}{N_{j}(t)}\right)$

where $N_{j}=S_{j}+E_{j}+I^{a}_{j}+I^{s}_{j}+R_{j}$ be the total population of
$j^{th}$ group, $\lambda_{ij}$ be the $ij^{th}$ element of the social contact
matrix which represents the number of contact between $i^{th}$ and $j^{th}$
age groups, and $\beta$ is the probability of infection of that contact. Here,
it is assumed that symptomatic infectious individuals reduce their contacts
compared to asymptomatic individuals, i.e., $\eta_{j}$ is the proportion of
contacts avoided by these self-isolating individuals of $j^{th}$ age group
(also allowing for compliance rates). It shows the awareness of symptomatic
infectious individuals. An increase in the value of $\eta_{j}$ shows decrease
the awareness of individuals of $j^{th}$ age group. Thus, by incorporating the
social contact matrix (contact patterns), the model system (1) takes the form
of following age-structured epidemic model (multi-group epidemic model) for
the $i^{th}$ age group.

$\displaystyle\frac{dS_{i}(t)}{dt}=$ $\displaystyle\Gamma_{i}-\beta
S_{i}(t)\sum_{j=1}^{n}\lambda_{ij}\left(\frac{I^{a}_{j}(t)+\eta_{j}I^{s}_{j}(t)}{N_{j}(t)}\right)-\mu_{i}S_{i},$
(3) $\displaystyle\frac{dE_{i}(t)}{dt}=$ $\displaystyle\beta
S_{i}(t)\sum_{j=1}^{n}\lambda_{ij}\left(\frac{I^{a}_{j}(t)+\eta_{j}I^{s}_{j}(t)}{N_{j}(t)}\right)-(\alpha_{i}+\mu_{i})E_{i}(t),$
$\displaystyle\frac{dI^{a}_{i}(t)}{dt}=$
$\displaystyle(1-\sigma)\alpha_{i}E_{i}(t)-\left(\gamma^{a}_{i}+\delta^{a}_{i}+\mu_{i}\right)I^{a}_{i}(t),$
$\displaystyle\frac{dI^{s}_{i}(t)}{dt}=$
$\displaystyle\sigma\alpha_{i}E_{i}(t)-\left(\gamma^{s}_{i}+\delta^{s}_{i}+\mu_{i}\right)I^{s}_{i}(t),$
$\displaystyle\frac{dR_{i}(t)}{dt}=$
$\displaystyle\gamma^{a}_{i}I^{a}_{i}(t)+\gamma^{s}_{i}I^{s}_{i}(t)-\mu_{i}R_{i}(t),$

where $i=1,2,\cdots,n$ and $n$ represents the number of age-groups in the
proposed model system. The biological meaning of parameters
$\Gamma_{i},\mu_{i},\alpha_{i},\sigma,\gamma^{a}_{i},\gamma^{s}_{i},\delta^{a}_{i}$
and $\delta^{s}_{i}$ for $i^{th}$ age-group are same as described for model
(1).

Figure 1: The schematic diagram for proposed model system (3) which represents
the transmission and transition between different compartments and groups.

######  Remark 2.1.

We can ignore the demographics of the population (i.e., birth and natural
deaths equal to zero) for the disease for which the infectious period is very
shorter compared with the lifespan of an individual of the considered
population. For example, in the case of the COVID-19 pandemic, we can ignore
the demographics as the infectious period ($7-14$ days) is shorter than the
lifespan of an individual (approximately $70$ years).

###  2.1 Characterization of the Social Contact Matrix

To determine the effect of different interventions, we divide the age-
structured contact matrix $\lambda=(\lambda_{ij})_{n\times n}$ into the
contributions from home location, schools (all educational institutions like
schools, colleges, universities, different types of coaching, etc.),
workplaces (govt. and private offices, other workplaces), and other locations
(markets, cinema halls, restaurants shopping malls, etc.), given by
$\lambda^{h},\lambda^{s},\lambda^{w}$ and $\lambda^{o}$ respectively. The
weight of each contact matrix is given by the coefficients
$\alpha_{h},\alpha_{s},\alpha_{w}$, and $\alpha_{o}$, we can change the
numerical value of these coefficients between zero and one overtime to reflect
the effect of different scenarios of interventions. If any of these locations
do not contribute to the disease transmission, we set the weight coefficient
for the specific contact matrix equal to zero. For a partial contribution, we
set the weight coefficient between 0 and 1, accordingly. Therefore, we can
divide the contact matrix $\lambda_{ij}$ in four different social contact
matrices with corresponding weight coefficients such that

$\lambda_{ij}=\alpha_{h}\lambda^{h}_{ij}+\alpha_{s}\lambda^{s}_{ij}+\alpha_{w}\lambda^{w}_{ij}+\alpha_{o}\lambda^{o}_{ij}$

For example, schools have been shut down during the strict lockdown in India,
so that, we can set $\alpha_{s}=0$ for this duration. It must also be
distinguished that the contributions to the work locations and other areas
like markets are never zero ( $\alpha_{w}\neq 0$ and $\alpha_{o}\neq 0$ ) even
during a strict lockdown because people were involved in essential services
and markets were opened functioning to a reduced degree. The contributions of
contacts at the home location are never zero. Moreover, it could also be
noticed that the different types of lockdown can encourage an increase in
contact at the home location because people are staying at home more. The
contributions of contacts at the home location can not be zero

##  3 Mathematical Analysis

In this section, we emphasize the existence, positivity, and boundedness of
solutions of model system (3). We also investigate the global dynamics around
the equilibria of the system.

###  3.1 Nonnegativity and Boundedness

Due to the biological feasibility of solutions of system (3), we are only
interested in nonnegative and bounded solutions of the system. Therefore, we
will show that all the solutions of the model system (3) with nonnegative
initial conditions are nonnegative i.e., $S_{i}(t)\geq 0,E_{i}(t)\geq
0,I^{a}_{i}(t)\geq 0,I^{s}_{i}(t)\geq 0,R_{i}(t)\geq 0$ for all $t\geq 0$ and
$i=1,2,\cdots,n.$

######  Theorem 3.1.

The solutions of model system (3) with the positive initial condition are
nonnegative for all $t\geq 0$ and uniformly bounded in the feasible region
$\Omega.$

###### Proof.

For the proof of Theorem 3.1, refer the Appendix A. ∎

###  3.2 Disease Free Equilibrium and Basic Reproduction Number

The disease free equilibrium (DFE) for model system (3) is

$E^{0}=\left(S^{0}_{1},0,0,0,0,S^{0}_{2},0,0,0,0,\cdots,S^{0}_{n},0,0,0,0\right)\in\Omega$

where $S^{0}_{i}=\frac{\Gamma_{i}}{\mu_{i}}$ and always exists. If the value
of $R_{0}$ is less than $1$, then the likelihood of producing new infection
cases by infectious people is insufficient for an outbreak to be sustained. If
$R_{0}$ is greater than $1$, the number of secondary cases increases the
infection and an epidemic arises until the proportion of susceptible
individuals declines. The next-generation matrix method in [31] is applied to
determine the basic reproduction number $(R_{0})$.
The basic reproduction number $(R_{0})$ is given by the following equation;

$R_{0}=\rho(FV^{-1})=\rho\left(F_{1}V_{21}V^{-1}_{11}V^{-1}_{22}+F_{2}V_{31}V^{-1}_{11}V^{-1}_{33}\right)$
(4)

Here matrix $F_{1}V_{21}V^{-1}_{11}V^{-1}_{22}$ gives the contribution from
asymptomatic infected individuals and matrix
$F_{2}V_{31}V^{-1}_{11}V^{-1}_{33}$ gives the contribution from symptomatic
infected individuals. For more detail calculations and
$F_{1},\,F_{2},V_{11},\,V_{22},\,V_{33},\,V_{21},\,V_{31}$ refer the Appendix
A.

###  3.3 Global Dynamics when $R_{0}\leq 1$

This section shows that the disease could be eliminated from a community when
the basic reproduction number $R_{0}\leq 1$ and whatever size of the initial
outbreak and that the infection persists otherwise.

######  Theorem 3.2.

Assume contact matrix $\lambda_{ij}$ is irreducible. If $R_{0}\leq 1$, then
the DFE $(E^{0})$ is globally asymptotically stable in $\Omega.$ If $R_{0}>1$,
then $E^{0}$ is unstable and system (3) is uniformly persistent and there
exists at least one endemic equilibrium (EE).

###### Proof.

See Appendix A. ∎

###  3.4 Global Dynamics when $R_{0}>1$

This section illustrates that the endemic equilibrium is unique and globally
asymptotically stable (GAS) in the interior of the feasible region $\Omega$.
Biologically, we say that the disease always remains in the population and
persists at a unique endemic level, whether the initial outbreak’s size is
small or larger. By Theorem 3.2, we have that an endemic equilibrium

$E^{*}=\left(\bar{S_{1}},\bar{S_{2}},\cdots,\bar{S_{n}},\bar{E_{1}},\bar{E_{2}},\cdots,\bar{E_{n}},\bar{I^{a}_{1}},\bar{I^{a}_{2}},\cdots,\bar{I^{a}_{n}},\bar{I^{s}_{1}},\bar{I^{s}_{2}},\cdots,\bar{I^{s}_{n}},\bar{R_{1}},\bar{R_{2}},\cdots,\bar{R_{n}}\right)$

exists and satisfy the following system of equations:

$\displaystyle\Gamma_{i}-\beta
S_{i}(t)\sum_{j=1}^{n}\lambda_{ij}\left(\frac{I^{a}_{j}(t)+\eta_{j}I^{s}_{j}(t)}{N_{j}(t)}\right)-\mu_{i}S_{i}=$
$\displaystyle 0,$ (5) $\displaystyle\beta
S_{i}(t)\sum_{j=1}^{n}\lambda_{ij}\left(\frac{I^{a}_{j}(t)+\eta_{j}I^{s}_{j}(t)}{N_{j}(t)}\right)-(\alpha_{i}+\mu_{i})E_{i}(t)=$
$\displaystyle 0,$
$\displaystyle(1-\sigma)\alpha_{i}E_{i}(t)-\left(\gamma^{a}_{i}+\delta^{a}_{i}+\mu_{i}\right)I^{a}_{i}(t)=$
$\displaystyle 0,$
$\displaystyle\sigma\alpha_{i}E_{i}(t)-\left(\gamma^{s}_{i}+\delta^{s}_{i}+\mu_{i}\right)I^{s}_{i}(t)=$
$\displaystyle 0,$
$\displaystyle\gamma^{a}_{i}I^{a}_{i}(t)+\gamma^{s}_{i}I^{s}_{i}(t)-\mu_{i}R_{i}(t)=$
$\displaystyle 0.$

######  Theorem 3.3.

Assume that the contact matrix $\lambda_{ij}$ is irreducible. The system (3)
has a unique and globally asymptotically stable (GAS) endemic equilibrium
$(E^{*})$ in $\Omega$ when $R_{0}>1$

###### Proof.

For the proof of Theorem 3.3, one may refer the Appendix A. ∎

##  4 The Optimal Control Problem Formulation

There were no treatments and vaccines for COVID-19 available in markets until
the end of 2020. As a result, many scientists emphasize mainly three control
strategies for challenging this pandemic and reducing the transmission risk
with Coronavirus. First, avoid exposure to COVID-19 by the following
prevention: reducing or avoiding unnecessary contact with infectious
individuals, using face masks, sensitizers, washing hands with water and soap
frequently, cleaning and disinfecting contagion surfaces using proper
chemicals. Second, stroking the people in quarantine areas who are suffering
from immunodeficiency then making special protection measure for them. Third,
putting infected individuals in isolation or making hospitalization to
decrease the infection risk of COVID-19.

Our main objective of this optimal control problem is to minimize the total
number of exposed people $(E_{i})$, the asymptomatic infectious people
$I^{a}_{i}$ and the symptomatic infectious people $I^{s}_{i}$. For this, we
include the three types of control in the controlled system (6). The first
control $\omega_{i}$ characterizes the effort of awareness programs in
$i^{th}$ age group that target presenting the control measures to prevent
individuals from being infected (washing hands regularly, using face masks,
avoiding unnecessary contacts with infected surfaces, and individuals).
Therefore, the term $(1-\omega_{i})$ is used to reduce the force of infection
by a proportion. The second control parameter $\theta_{i}$ represents the
effort to encourage the exposed individuals of $i^{th}$ age group to join
quarantine centers. Finally, control parameter $\vartheta$ measures the effort
to aware the infected individuals to stay home in isolation or join the
institutional isolation ward. So, the mathematical system with control
parameters is given by the following system of differential equations.

$\displaystyle\frac{dS_{i}(t)}{dt}=$
$\displaystyle\Gamma_{i}-\beta(1-\omega_{i}(t))S_{i}(t)\sum_{j=1}^{n}\lambda_{ij}\left(\frac{I^{a}_{j}(t)+\eta_{j}I^{s}_{j}(t)}{N_{j}(t)}\right)-\mu_{i}S_{i},$
(6) $\displaystyle\frac{dE_{i}(t)}{dt}=$ $\displaystyle\beta
S_{i}(t)(1-\omega_{i}(t))\sum_{j=1}^{n}\lambda_{ij}\left(\frac{I^{a}_{j}(t)+\eta_{j}I^{s}_{j}(t)}{N_{j}(t)}\right)-(\alpha_{i}+\mu_{i})E_{i}(t)-\theta_{i}(t)E_{i}(t),$
$\displaystyle\frac{dI^{a}_{i}(t)}{dt}=$
$\displaystyle(1-\sigma)\alpha_{i}E_{i}(t)-\left(\gamma^{a}_{i}+\delta^{a}_{i}+\mu_{i}\right)I^{a}_{i}(t)-\vartheta(t)I^{a}_{i}(t),$
$\displaystyle\frac{dI^{s}_{i}(t)}{dt}=$
$\displaystyle\sigma\alpha_{i}E_{i}(t)-\left(\gamma^{s}_{i}+\delta^{s}_{i}+\mu_{i}\right)I^{s}_{i}(t)-\vartheta(t)I^{s}_{i}(t),$
$\displaystyle\frac{dRi(t)}{dt}=$
$\displaystyle\gamma^{a}_{i}I^{a}_{i}(t)+\gamma^{s}_{i}I^{s}_{i}(t)-\mu_{i}R_{i}(t),$

where $S_{i}(0)\geq 0,E_{i}(0)\geq 0,I^{a}_{i}(0)\geq 0,I^{s}_{i}(0)\geq 0,$
and $R_{i}(0)\geq 0$ for all $i=1,2,\cdots,n$ are given initial sizes of
variables. The main objective of the optimal control problem is to compare the
costs of these intervention strategies and their effectiveness in the fight
against the disease. To do this, we need to investigate the optimal level of
effort that would be required to control the disease or to minimize the total
number of infected individuals. For this, we use the following objective
function :

$\displaystyle J(\omega_{i},\theta_{i},\vartheta_{i})=$ (7) $\displaystyle
E_{i}(T)+I^{a}_{i}(T)$
$\displaystyle+I^{s}(T)+\int_{0}^{T}\left(E_{i}(t)+I^{a}_{i}(t)+I^{s}(t)+\frac{A^{1}_{i}}{2}\left(\omega_{i}(t)\right)^{2}+\frac{A^{2}_{i}}{2}\left(\theta_{i}(t)\right)^{2}+\frac{A^{3}_{i}}{2}\left(\vartheta_{i}(t)\right)^{2}\right)dt,$

where parameters $A^{1}_{i}>0,A^{2}_{i}>0,$ and $A^{3}_{i}>0,$ for all
$i=1,2,\cdots,n$ are the cost coefficients at time $t$ associated with applied
controls , $T$ is the final time. On the other hands, we try to find the
optimal controls $\omega^{*}_{i},\,\theta^{*}_{i},$, and $\vartheta^{*}_{i}$
such that

$J(\omega^{*}_{i},\theta^{*}_{i},\vartheta^{*}_{i})=\min_{(\omega_{i},\theta_{i},\vartheta_{i})\in
C_{ad}}J(\omega_{i},\theta_{i},\vartheta_{i})$

where $C_{ad}$ is the set of acceptable control defined by

$C_{ad}=\left\\{(\omega_{i},\theta_{i},\vartheta_{i}):0\leq\omega_{i}(t)\leq
1;0\leq\theta_{i}(t)\leq 1;0\leq\vartheta_{i}(t)\leq 1,t\in[0,T]\right\\}$

### Characterization of the optimal control

Here, we compute the necessary condition for the optimal control by using the
pontryagin’s maximum principle [42]. The idea of Pontryagin’s maximum
principle makes known to the adjoint function to the associated system to the
objective functional resulting in the formulation of a function called the
Hamiltonian. The principle interprets into a problem of minimizing Hamiltonian
$H_{i}(t)$ at time $t$ defined by

$H_{i}(t)=E_{i}(t)+I^{a}_{i}(t)+I^{s}(t)+\frac{A^{1}_{i}}{2}\left(\omega_{i}(t)\right)^{2}+\frac{A^{2}_{i}}{2}\left(\theta_{i}(t)\right)^{2}+\frac{A^{3}_{i}}{2}\left(\vartheta_{i}(t)\right)^{2}+\sum_{k=1}^{5}\psi^{k}_{i}(t)f_{k}\left(S_{i},E_{i},I^{a}_{i},I^{s}_{i},R_{i}\right),$
(8)

where $\psi^{k}_{i}$ is the $k^{th}$ adjoint variable at time $t$ and $f_{k}$
is the right-hand side function of the control system (6) corresponding to the
$k^{th}$ state at time $t.$

######  Theorem 4.1.

Given the optimal controls $\omega^{*}_{i},\theta^{*}_{i},\vartheta^{*}_{i}$
and the solutions $S_{i},E_{i},I^{a}_{i},I^{s}_{i}$ and $R_{i}$ associated to
control system (6), then there exist adjoint variables $\psi^{k}_{i}$ for
$k=1,2,3,4,5$ satisfying

$\displaystyle(\psi^{1}_{i})^{\prime}==$
$\displaystyle\psi^{1}_{i}\mu_{i}+(\psi^{1}_{i}-\psi^{2}_{i})\beta\left(1-\omega_{i}(t)\right)\sum_{j=1}^{n}\lambda_{ij}\left(\frac{I^{a}_{j}(t)+\eta
I^{s}_{j}(t)}{N_{j}}\right),$ (9) $\displaystyle(\psi^{2}_{i})^{\prime}==$
$\displaystyle\left(\alpha_{i}+\mu_{i}+\theta_{i}(t)\right)\psi^{2}_{i}-(1-\sigma)\alpha_{i}\psi^{3}_{i}-\sigma\alpha_{i}\psi^{4}_{i},$
$\displaystyle(\psi^{3}_{i})^{\prime}==$
$\displaystyle(\psi^{1}_{i}-\psi^{2}_{i})\beta(1-\omega_{i}(t))S_{i}(t)\sum_{j=1}^{n}\left(\frac{\lambda_{ij}}{N_{j}}\right)+\left(\gamma^{a}_{i}+\delta^{a}_{i}+\mu_{i}+\vartheta(t)\right)\psi^{3}_{i}+\gamma^{a}_{i}\psi^{5}_{i},$
$\displaystyle(\psi^{4}_{i})^{\prime}==$
$\displaystyle(\psi^{1}_{i}-\psi^{2}_{i})\beta(1-\omega_{i}(t))S_{i}(t)\eta_{j}\sum_{j=1}^{n}\left(\frac{\lambda_{ij}}{N_{j}}\right)+\left(\gamma^{s}_{i}+\delta^{s}_{i}+\mu_{i}+\vartheta(t)\right)\psi^{4}_{i}+\gamma^{s}_{i}\psi^{5}_{i},$
$\displaystyle(\psi^{5}_{i})^{\prime}==$ $\displaystyle\mu_{i}\psi^{5}_{i},$

with the tranversality conditions at time
$T:\psi^{1}_{i}(T)=0,\psi^{2}_{i}(T)=1,\psi^{3}_{i}(T)=1,\psi^{4}_{i}(T)=1,\psi^{5}_{i}(T)=0.$
Then the optimal controls $\omega^{*}_{i},\theta^{*}_{i},$ and
$\vartheta^{*}_{i}$ are given by

$\omega^{*}_{i}=\max\left(0,\min\left(1,\frac{\left(\psi^{2}_{i}-\psi^{1}_{i}\right)\beta
S_{i}(t)}{A^{1}_{i}}\sum_{j=1}^{n}\lambda_{ij}\left(\frac{I^{a}_{j}(t)+\eta_{j}I^{s}_{j}(t)}{N_{j}}\right)\right)\right),$
$\theta^{*}_{i}=\max\left(0,\min\left(1,\frac{\psi^{2}_{i}E_{i}(t)}{A^{2}_{i}}\right)\right),\,\,\text{and}\,\,\vartheta^{*}_{i}=\max\left(0,\min\left(1,\frac{\psi^{3}_{i}I^{a}_{i}(t)+\psi^{4}_{i}I^{s}_{i}(t)}{A^{3}_{i}}\right)\right).$

###### Proof.

For the proof of Theorem 4.1, readers may refer the Appendix B ∎

##  5 Numerical Simulations

Here, we illustrate the numerical results to our proposed multi-group model
(3).

###  5.1 Comparison of contact pattern for four countries (India, Italy,
Brazil and USA) at different locations

On August 20, 2020, the Ministry of Home Affairs (MHA), India, released the
rules for activities allowed in Unlock $4.0$ [43]. It states that lockdown
will continue in the high-risk zones till $30^{th}$ September 2020, the same
as the previous month (August). Outside the containment region, some events
were allowed in a specific manner. The metro rail was permitted to be
restarted in the categorized mode from September 7, 2020. The gathering of
limited peoples was allowed; fifty people were allowed in marriage functions,
twenty people were legalized in funereal/last rites ceremonies, and up to one
hundred people were acceptable in religious, political, sports, and academic
functions. Face masks wearing was mandatory in markets, workplaces, and during
traveling. On September 30, 2020, the MHA delivered the guidelines for
activities that were allowed in Unlock $5.0$ [44] and announced that lockdown
should continue in force stringently in the containment region still September
30, 2020 [48]. States and Union Territories will be able to take some
necessary actions from $15^{th}$ October to reopen the schools in a group
manner. Cinema halls could be reopened from $10^{th}$ October 2020 with the
half seating of their capacity, and Swimming pools were allowed to be used for
training of sportsperson [44]. On $27^{th}$ October, 2020 the MHA announced
the instructions for Unlock $6.0$ [45, 46, 47, 48] and said the guidelines of
Unlock 5.0 would continue to be employed in November too [13]. Out of
containment zones, some states have allowed more activities and partial
reopening of schools [49]. Thus, the Ministry of Home Affairs released
separate guidelines to each Unlock every month. The movements of people vary
according to the guidelines of every Unlock. Therefore, the contribution of
social contacts in the disease transmissions depends on the Unlock.

First of all, we compare the social contact structure of people of four
countries (India, Italy, Brazil, and USA). The data for the age distributions
of people for these countries are collected from the Population Pyramid
website [50] as depicted in Fig. 9. In these age distributions, after
selecting for the different age groups, the age of each individual is
distributed uniformly from the associated range. The social contact matrices
for different locations (home, work, school, others) are found in Prem et. al
[51] that are obtained from surveys and Bayesian imputation and represented by
Fig. 10.
This comparison aims to highlight the differences between their contact
patterns and emphasize their effects on the transmission of infectious
diseases. Fig. 9 shows the percentage of the population (separated by gender)
in the five-year age group terminating at the age of eighty. The Taj Mahal
dome shape of age distributions for India and Brazil is typical of a
demographic transition. The narrower base of Italian populations is
characteristic of aging people at or near sub-replacement fertility. The
rectangle shape of the age distribution of the USA population represents the
stationary distribution.
The second row of Fig. 10 shows the contact between different age groups in
the household setting, represented by matrices $\lambda^{h}_{ij}$ where darker
magenta colored squares indicate a more significant number of contacts. As
discussed in [51], the structures common to India, Italy, and Brazil are the
diagonal dominance, reflecting communication between same age groups (i.e.,
siblings and partners) and the prominent off diagonals, separated by the mean
inter-generation gap, reflecting connections between different age groups
(i.e., children and parents). However, we have observed some other types of
behavior in contact patterns at the home location in the case of USA. The
principal difference in the contact matrix for India is the presence of a
third dominant diagonal, again separated by the mean inter-generation gap,
reflecting the incidence of three-generation households. It quantifies the
significant contact between children and grandparents and the possibility of
substantial transmission of contagion from third to the first generation. Such
connections are smaller in Brazil and negligible in Italy.
The third row of Fig. 10 shows the contact patterns $\lambda^{s}_{ij}$ between
age groups at the school location in which we see that there is main diagonal
dominant till the school ages (strongly assortative, with primary contacts
within the school-going children), reflecting contact within school-age
groups, and other smaller contacts reflecting the connection between student
and teacher. As observed, this structure is typical in all four countries. The
fourth row of Fig. 10 shows the social contact matrices $\lambda^{w}_{ij}$
between age groups in the workplace. In disparity to home and school location,
the contact pattern at the workplace is more homogeneous across age groups in
all four countries, indicating that the workplace contributes to the
transmission of contagion between age groups that are, otherwise, largely
separated from each other in the household. The boundaries of these age groups
are more prominent in India and USA than in Italy and Brazil. The last row of
Fig. 10 shows the matrices $\lambda^{o}_{ij}$ for contact patterns at other
locations. It shows that the contacts are strongly assortative (forcefully
contact within age groups) for India, Brazil, and USA, reflecting the
preferential social connection within age groups in this sphere, but otherwise
do not show systematic patterns. However, the contact pattern matrix for Italy
reflects the non-preferential social connection.
In summary, in the case of India, the home location provides the main channel
to disease transmission between three generations, and the workplace provides
the main channel to (largely homogeneous) disease transmission between
working-age groups. The school is the main channel of infection within
children and between children and adult teachers to a smaller extent. Other
spheres of contact, due to the assortative mixing, contribute to transmission
within age groups.

Social mixing patterns vary according to locations, including homes,
workplaces, schools, and other locations (markets, worship places, etc.).
Therefore, we used the method set out in Prem et al. [51] which takes accounts
for these differences and obtain the location-specific social contact matrices
for various countries. In a standard setting, the average number of contacts
made at all of these locations contributes to the overall mixing pattern in
the population, i.e., the social contact matrices of the first row of Fig. 10
can be used for standard setting (before the outbreak). In an outbreak
duration, different control measures (non-pharmaceutical interventions) are
applied at reduced social mixing in different situations to lower the overall
infection spreading in the community. We used the weight coefficients between
$0$ and $1$ to incorporate each intervention scenario from these building
block matrices to simulate the effects of various control measures applied at
reduced social mixing.

###  5.2 Model fitting and parameter estimation: A case study of COVID-19 in
India

Since, from September 1, 2020, to December 31, 2020, the nationwide lockdown
was lifted in a monthly phased manner. Every month, the strictness of lockdown
was condensed, and the government of India permitted some more activities for
people. On August 20, 2020, the Ministry of Home Affairs (MHA), India,
released the rules for activities allowed in Unlock $4.0$ [43]. It states that
lockdown will continue in the high-risk zones till $30^{th}$ September 2020,
the same as the previous month (August). Outside the containment region, some
events were allowed in a specific manner. The metro rail was permitted to be
restarted in the categorized mode from September 7, 2020. The gathering of
limited peoples was allowed; fifty people were allowed in marriage functions,
twenty people were legalized in funereal/last rites ceremonies, and up to one
hundred people were acceptable in religious, political, sports, and academic
functions. Face masks wearing was mandatory in markets, workplaces, and during
traveling. On September 30, 2020, the MHA delivered the guidelines for
activities that were allowed in Unlock $5.0$ [44] and announced that lockdown
should continue in force stringently in the containment region still September
30, 2020 [48]. States and Union Territories will be able to take some
necessary actions from $15^{th}$ October to reopen the schools in a group
manner. Many activities were accepted in a particularly phased manner [44]. On
$27_{th}$ October 2020, the MHA announced the instructions for Unlock $6.0$
[45, 46, 47, 48] and said, that the guidelines of Unlock 5.0 would continue to
be employed in November too. Out of containment zones, some of the states have
allowed more activities and partial reopening of schools [49]. Thus, the MHA
released separate guidelines to each Unlock every month.
The movements of people vary according to the guidelines of every Unlock.
Therefore, the contribution of social contacts in the disease transmissions
depends on the Unlock. Hence, the social contact pattern was changing
continuously, and the transmission rate differed accordingly. Based on these
observations, we consider the following four lockdown lifting phases: (i)
September 2020 (ii) October 2020 (iii) November 2020 (iv) December 2020 and
estimate the value of the transmission rate of COVID-19 infection for each
phase. Therefore, we choose the empirical data of COVID1-19 for India for the
period from September 1, 2020, to December 31, 2020 [54], to estimate the
transmission rate of COVID-19 for the considered model. Our model fits the
cumulative cases data for the period mentioned earlier, and the value of
transmission rate $\beta$ has been estimated. It is noteworthy that the
strictness of the intervention strategies can essentially change the disease’s
transmission dynamics.

Due to the short study duration of COVID-19 epidemics compared to the human
lifespan, we ignore the new birth rate in susceptible individuals and natural
death rates of all individuals (demographics) i.e., we consider $\Gamma=0$ and
$\mu=0$. We assume that the $40\%$ individuals of age groups $0-19$ years and
the $60\%$ individuals of age groups $20-79$ years, are showing symptoms after
the incubation period [52]. So, we set $\sigma=0.4$ for $i=1,2,3,4$ and
$\sigma=0.6$ for $i=5,6,...,16.$ In terms of parametrization for incubation
period of COVID-19, a reliable estimate of $\alpha_{i}$ for COVID-19 is
difficult to obtain from observed data of COVID-19 infection. Therefore, we
assume $\alpha_{i}$ to be a constant among all age groups, i.e.,
$\alpha_{i}=\alpha$ for $i^{th}$ age group. Thus, we set
$\alpha_{i}=\alpha=0.196$ (i.e. $1/\alpha(\text{latent period})=5.1$ [53]).
For similar reason, further, we aslo assume $\gamma^{a}_{i}$ and
$\gamma^{s}_{i}$ to be a constant $\gamma^{a}_{i}=\gamma^{s}_{i}=\gamma_{i}$
among all age groups, i.e., $\gamma_{i}=\gamma=\frac{1}{\text{Infectious
period}}$ for $i^{th}$ age group. Therefore, the following system of
differential equations COVID-19 dynamic.

$\displaystyle\frac{dS_{i}(t)}{dt}=$ $\displaystyle-\beta
S_{i}(t)\sum_{j=1}^{n}\lambda_{ij}\left(\frac{I^{a}_{j}(t)+\eta_{j}I^{s}_{j}(t)}{N_{j}(t)}\right),$
(10) $\displaystyle\frac{dE_{i}(t)}{dt}=$ $\displaystyle\beta
S_{i}(t)\sum_{j=1}^{n}\lambda_{ij}\left(\frac{I^{a}_{j}(t)+\eta_{j}I^{s}_{j}(t)}{N_{j}(t)}\right)-\alpha
E_{i}(t),$ $\displaystyle\frac{dI^{a}_{i}(t)}{dt}=$
$\displaystyle(1-\sigma_{i})\alpha
E_{i}(t)-\left(\gamma+\delta^{a}_{i}\right)I^{a}_{i}(t),$
$\displaystyle\frac{dI^{s}_{i}(t)}{dt}=$ $\displaystyle\sigma_{i}\alpha
E_{i}(t)-\left(\gamma+\delta^{s}_{i}\right)I^{s}_{i}(t),$
$\displaystyle\frac{dR_{i}(t)}{dt}=$
$\displaystyle\gamma\left(I^{a}_{i}(t)+I^{s}_{i}(t)\right),$

where $i\in\\{1,2,\cdots,16\\}$ and $n=16$ represents the number of age groups
for simulations purpose of COVID-19 for India and the social contact matrix
$\lambda_{ij}$ is parameterized as

$\lambda_{ij}=\alpha_{h}\lambda^{h}_{ij}+\alpha_{s}\lambda^{s}_{ij}+\alpha_{w}\lambda^{w}_{ij}+\alpha_{o}\lambda^{o}_{ij}.$

Further, we have that the limited activities are allowed to people in
September 2020. Thus, we assume that $30\%$ social contacts at workplaces and
$20\%$ contacts at other places of the normal setting are possible and
generate the synthetic social contact matrices for workplaces and other places
using weight coefficients $\alpha_{w}=0.3$ and $\alpha_{o}=0.2$ for September
2020, respectively. In similar manner, we set weight coefficients
$\alpha_{w}=0.5$ and $\alpha_{o}=0.4$ for October; $\alpha_{w}=0.6$ and
$\alpha_{o}=0.65$ for November; $\alpha_{w}=0.65$ and $\alpha_{o}=0.75$ for
December, 2020 according to guidelines of lockdown lifting. We collected the
real data of COVID-19 for India from September 1, 2020, to December 31, 2020,
from the website [54]. Further, we perform our model fitting by using the in-
built function lsqnonlin in Matlab 2018 to minimize the sum of square function
and use the following steps;
Algorithm for Least Square Method for Data Fitting:
Step 1. We obtain the cumulative cases of model output at time $t$ by
summation of cumulative cases of each age group at time $t$. It is given by

$I_{cum}(t,\Theta)=\sum_{i=1}^{n}\int_{0}^{t}\alpha_{i}E_{i}(\tau)d\tau,\,\,\,\,\text{where}\,\theta=\\{\beta\\},$

where $I_{cum}(t,\Theta)$ represents the cumulative confirmed cases at time
$t$ and $n$ be the number of age groups. For the initial values of the
population of each age group, we assume that the total initial population of a
class could be divided into age-specific compartments according to the age
distribution pyramid discussed in Fig. 9.
Step 2. We compute the sum of squares of errors at each time point, given by

$SS_{k}(\Theta)=\sum_{j=1}^{m}\left(I_{cum}^{a}(t_{j})-I_{cum}(t_{j},\Theta)\right)^{2},$

where, $I_{cum}^{a}(t_{j})$ is the actual data at $t_{j}^{th}$ day and $m$ is
the number of data points. $k$ represents the number of iterations performed
such that $SS_{k}$ is the sum of squares of errors in the $k^{th}$ iteration.
Step 3. We compute the value of parameter $\Theta$ such that

$SS=\min\left\\{SS(\Theta)\right\\}.$

Thus, we obtain the value of the transmission rate $(\beta)$ for the best
fitting of model output to actual data. The fitting of model output to real
data of cumulative confirmed cases is shown in Fig. 2. The estimated values of
the transmission rate $\beta$ are displayed in Fig. 3.

Figure 2: The best fitting of model output to cumulative confirmed cases. The
month-wise fitting incorporates different interventions scenarios applied by
MHA in India. Figure 3: The figure illustrates the numerical values of model
parameter $\beta$ for different age groups to the best fitting of model output
according to previous months (September, October, November, and December) of
2020 (different-different scenarios of unlocking in India).

###  5.3 Age profile of Asymptomatic and Symptomatic Infected Individuals

Figure 4: The figure depicts the long term dynamics of infected individuals of
different age groups.

The age profile of infected individuals (asymptomatic and symptomatic ) are
shown in Fig. 4, in which we observe that asymptomatic infected individuals of
age groups 05-24 years are gradually increasing. The possible reason behind
this type of dynamics of infected individuals is the unawareness of children.
The asymptomatic infected individuals of age groups higher than 60 years
increase and then decrease to stable levels due to their limited contacts
(fewer contacts) to other persons. Moreover, the asymptomatic infected
individuals of middle age groups (25-59 years) increase to their saturation
levels. The individuals aged 25-59 years could not have important contacts due
to their jobs, but they use protection measures like face masks, social
distancing, etc. Therefore, the dynamics of infected individuals of age groups
25-59 years is very different from those of 10-24 years. The same reasoning is
also applicable for the dynamics of symptomatic infected individuals (please
refer to the right panel of Fig. 4).

###  5.4 Impact of Awareness of Symptomatic Individuals

From the expression (2) (which represents the basic reproduction number
$(R_{0})$ of the single group $SEI_{a}I_{s}R$ epidemic model), we noticed that
$R_{0}$ would decrease by decreasing the values of $\eta$, i.e., the increase
in the awareness in symptomatic individuals implies an increase in the
avoiding unnecessary contacts. Consequently, the disease burden reduces to a
significant level. Since the expression of the basic reproduction of system
(3) has a complicated form. Therefore, it is not easy to discuss the impact of
awareness in Symptomatic individuals on basic reproduction number $R_{0}$
(expressed in equation (4)) of age-structured epidemic model (3) explicitly.
Therefore, we observe the effects of awareness of symptomatic individuals by
time series analysis of total infected individuals with various levels of
understanding of people in Fig. 5 when schools are closed, i.e.,
$\alpha_{s}=0$ and Fig. 6 when schools are partially opened, i.e.,
$\alpha_{s}>0$. In Fig. 5a, we observe the effects of symptomatic infected
individuals of age groups $0-19$ years and notice that the number of total
infected individuals is reduced by $20\%$ (approximately) at the end of the
year due to increases in the awareness of symptomatic infected individuals of
age groups $0-19$ years by $40\%$. From Fig. 5b, we obtain that the number of
total infected individuals is reduced by $33\%$ (approximately) at the end of
the year due to increases in the awareness of symptomatic infected individuals
of age groups $20-49$ years by $50\%$. In Fig. 5c, we observe that the changes
in the number of total infected individuals are significantly less
(approximately $10\%$ reduction) at the end of the year due to increases in
the awareness of symptomatic infected individuals of age groups $50-79$ years
by $50\%$. Approximately $40\%$ reduction in the number of the total infected
population is noticed due to a $50\%$ increment in awareness of symptomatic
infected individuals of age groups $20-79$ years (please refer the Fig. 5d).

(a)

(b)

(c)

(d)

Figure 5: The figure illustrates the impacts of awareness of symptomatic
infected individuals on the long term dynamics of infected individuals when
schools are closed $(\alpha_{s}=0).$

(a)

(b)

(c)

(d)

Figure 6: The figure illustrates the impacts of awareness of symptomatic
infected individuals on the long term dynamics of infected individuals when
schools are partially opened $(\alpha_{s}>0).$

From Fig. 5a and 6a, we have that the awareness of children is the most potent
control measure to decrease the number of infected people when schools are
partially opened. Therefore, the awareness programs for children should be
exhibited by the authority to reduce the endemic level to a significant level
when schools have opened partially.

###  5.5 Impact of Social mixing pattern

(a)

(b)

(c)

(d)

Figure 7: The figure illustrates the impacts of social mixing pattern of
people on the long term dynamics of infected individuals.

In Fig. 7, we plotted long-term dynamics of total infected individuals for
various social mixing scenarios. In the first three figures of 7 (Fig. 7a, 7b
and 7c), we observe that the total number of infected individuals for four
different levels of social contacts of a particular location and the weight
coefficient of other three contacts matrices are fixed at a positive value. In
Fig. 7a, we observe that the effects of various levels of school reopening
during the lockdown lifting when other three social contacts matrices are
contributing to disease transmission, i.e., $\alpha_{h}>0,\alpha_{w}>0$ and
$\alpha_{o}>0$ are taken in the figure. Here, we notice that if $60\%$ schools
are opened, the number of infected individuals is raised by approximately
$96\%$ from that of obtained in case of totally schools closure at the end of
the year. About $35\%$ increment in the number of infected individuals is
reported at the end of the year when $60\%$ contacts of normal setting at
working places are considered.

##  6  Discussion

COVID-19, a contact-contagious infectious disease, is assumed to spread
through a population via direct contacts between peoples [2, 57]. Outbreak
control measures aimed at reducing the amount of mixing in the population can
reduce the final size of the epidemic. Mathematical models help us to
understand that how COVID-19 could spread across the population and inform
control measures that might mitigate future transmission[56, 57]. We simulated
the shape of the ongoing outbreak of COVID-19 in India using an age-structured
multi-group SEIR model. We accounted for various types of heterogeneities in
contact networks in our model system (3) to determine the effect of location-
specific physical distancing measures such as extended school closures and
interventions in workplaces and other gathering places on the disease
progression. We simulated the COVID-19 outbreak and modeled the interventions
by scaling down the appropriate component of the contact mixing matrices for
India. As individuals mixing patterns are non-random, they influence the
disease transmission dynamics [20, 26], models that evaluate the effectiveness
of physical distancing interventions, such as school closure, need to account
for social structures and heterogeneities in the mixing of individuals [58,
59, 60].

In our proposed model, we combined the changes to age-specific and location-
specific social mixing patterns to estimate the effects of location-specific
physical distancing interventions in controlling the transmission of the
outbreak. The control measures put in place to reduce contacts at schools,
workplaces, and other gathering places help us prevent the epidemic by
providing the healthcare system sufficient time and opportunity to respond.
Subsequently, if these restrictions are lifted early, there are still adequate
susceptible people to keep the $R_{0}>1,$ once contacts increase and the
number of infections grows. However, in realistic scenarios, the interventions
are lifted slowly in a specific manner, partially, as an effort to avoid a
sharp growth in infection due to logistical and practical reasons. Therefore,
we simulated variously lifting the interventions such as school reopening in
phase manners, market reopening, and many daily activities under the
particular guidelines. There is little evidence of the effects of different
physical distancing measures on containing the outbreak, is known about the
behavioral changes of individuals over time during the outbreak. Therefore, we
supposed the effect that certain types of physical distancing have on age-
specific and location-specific contact rates to model the impact of the
physical distancing measures executed in India.

To mathematical analysis, we have shown the existence of solutions for the
proposed model. The nonnegativity and boundedness of solutions also have been
discussed to demonstrate the feasibility in the context of population
dynamics. The basic reproduction number $(R_{0})$ has been computed using the
next-generation matrix approach. We have studied the global dynamics of the
system (3) in terms of the basic reproduction number to check the behavior of
the system around its equilibria. We found that the disease eliminates from
all age groups when $R_{0}<1$ and the disease remains in the population when
$R_{0}>1$. Hence, we obtained that the basic reproduction number plays the
role of a threshold quantity to disease elimination from all age groups, and
$R_{0}<1$ is a sufficient condition to eliminate COVID-19. The optimal control
problem formulation has also been discussed by considering three non-
pharmaceutical controls, namely: the first is the effort to awareness
campaigns program regarding taking care of their health to protect from
COVID-19, and second is the effort to encourage the exposed individuals to
join the quarantine centers to break the transmission of COVID-19. The last
control is the effort to make aware the infected individuals stay home in
isolation or join the institutional isolation ward. The existence of optimal
control has been discussed for the proposed optimal control problem (19). The
optimal control to minimize the total number of infected population also has
been computed mathematically.

Our numerical results exhibited that control measures aimed at reducing social
mixing in the population can effectively reduce the disease burden of the
COVID-19 outbreaks. For different scenarios of lifting the control measures
among individuals, the schools and other places reopening have had a
significant impact on the progression of the epidemic. We present the age
profile of the number of asymptomatic and symptomatic infected individuals for
one year (365 days) in Fig. 4. Our simulation result suggests that the rate of
change in the number of asymptomatic and symptomatic infected individuals vary
in different categories; the number of asymptomatic infected individuals of
age groups $5-24$ year and the number of symptomatic infected individuals of
age groups $20-44$ year fastly increase. The solution trajectories of other
age groups are going to an almost saturation level. Further, we discussed the
effects of awareness of symptomatic infected individuals (refer Fig.5) when
all schools are fully closed, i.e., $\alpha_{s}=0.$ Here, we conclude that the
awareness of symptomatic infected individuals with age groups $20-49$ year is
beneficial to reduce the number of infected individuals. It should be noticed
that the large number of people involved in the essential services during the
COVID-19 pandemic belongs to age groups $20-49$ year. The contacts of these
people were contributing to the transmission of COVID-19 on a large scale.
Therefore, the disease burden may have been reduced significantly by their
awareness. It becomes essential that healthcare workers and other service
providers are trained for awareness differently. However, the reduction in
contacts of school-age children also reduces infection in the community when
all or some schools are partially opened, i.e., $0<\alpha_{s}<1.$ ( refer the
Fig. 6). The impact of different scenarios of social mixing patterns of people
on the dynamics of infected individuals has been shown in Fig 7.

The main findings of this study are that the reduction of contacts at schools
and other locations (market, shopping malls, and other gathering places) is an
effective tool to decrease the number of infected individuals to a significant
level. The school closure is also suggested as an intervention to reduce the
COVID-19 in Wuhan, China [20]. The school reopening with limited attendance,
such as older students sitting exams or the partial reopening to younger year
groups, do not appear to affect the community spreading of COVID-19
significantly. However, fully school reopening can contribute to statistically
significant increases in the transmission rate for COVID-19 [61]. Face mask-
wearing also should be necessary for the public and aware to avoid unnecessary
contact with other people. In the case of schools closure, the awareness of
symptomatic infected individuals with age groups $20-49$ years is helpful to
significantly decrease the number of infected individuals due to their
involvement in the essential services during the COVID-19 pandemic. However,
In the case of schools partially opening, children’s awareness is also
required to reduce the number of infected people and the awareness of
symptomatic infected individuals with age groups $20-49$ years. Similar
results were also discussed for United States in [62].

In conclusion, the non-pharmaceutical interventions based on sustained
physical distancing have a strong potential to reduce the size of the epidemic
burden of COVID-19 and lead to a smaller number of overall cases. The pressure
on the healthcare system could be reduced by lowering and flattening the
epidemic curve. The reduction in the disease prevalence also helps to reduce
the load on the healthcare system. The premature and sudden lifting of control
measures could lead to other peaks, and the endemic curve could be flattened
by gradually relaxing the interventions. Therefore, gradually relaxing the
interventions is mandatory to reduce the disease burden in the community.

## References

  * [1] Grassly, Nicholas C., and Christophe Fraser. ”Mathematical models of infectious disease transmission.” Nature Reviews Microbiology 6.6 (2008): 477-487.
  * [2] Gog, Julia R. ”How you can help with COVID-19 modelling.” Nature Reviews Physics 2.6 (2020): 274-275.
  * [3] Vespignani, Alessandro, Huaiyu Tian, Christopher Dye, James O. Lloyd-Smith, Rosalind M. Eggo, Munik Shrestha, Samuel V. Scarpino et al. ”Modelling covid-19.” Nature Reviews Physics 2, no. 6 (2020): 279-281.
  * [4] Thompson, Robin N., T. Déirdre Hollingsworth, Valerie Isham, Daniel Arribas-Bel, Ben Ashby, Tom Britton, Peter Challenor et al. ”Key questions for modelling COVID-19 exit strategies.” Proceedings of the Royal Society B 287, no. 1932 (2020): 20201405.
  * [5] Duerr, H-P., M. Schwehm, C. C. Leary, S. J. De Vlas, and M. Eichner. ”The impact of contact structure on infectious disease control: influenza and antiviral agents.” Epidemiology & Infection 135, no. 7 (2007): 1124-1132.
  * [6] Pellis, Lorenzo, Neil M. Ferguson, and Christophe Fraser. ”Threshold parameters for a model of epidemic spread among households and workplaces.” Journal of the Royal Society Interface 6, no. 40 (2009): 979-987.
  * [7] Becker, Niels G., and Klaus Dietz. ”The effect of household distribution on transmission and control of highly infectious diseases.” Mathematical Biosciences 127, no. 2 (1995): 207-219.
  * [8] Kiesha, P., Y. Liu, and T. Russell. ”The effect of controlling strategies to reduce social mixing on outcomes of the COVID-19 epidemic in Wuhan.” China: A modeling study, LANCET PUBLIC Health, (20) (2020): 30073-6.
  * [9] Baguelin, Marc, Stefan Flasche, Anton Camacho, Nikolaos Demiris, Elizabeth Miller, and W. John Edmunds. ”Assessing optimal target populations for influenza vaccination programmes: an evidence synthesis and modelling study.” PLoS Med 10, no. 10 (2013): e1001527.
  * [10] de Waroux, Olivier le Polain, et al. ”Identifying human encounters that shape the transmission of Streptococcus pneumoniae and other acute respiratory infections.” Epidemics 25 (2018): 72-79.
  * [11] Latsuzbaia, Ardashel, Malte Herold, Jean-Paul Bertemes, and Joël Mossong. ”Evolving social contact patterns during the COVID-19 crisis in Luxembourg.” PloS one 15, no. 8 (2020): e0237128.
  * [12] A. J. Kucharski, K. O. Kwok, V. W. I. Wei, B. J. Cowling, J. M. Read, J. Lessler, D. A. Cummings, and S. Riley. The Contribution of Social Behaviour to the Transmission of Influenza A in a Human Population. PLoS Pathog, 10(6):e1004206, 2014.
  * [13] R. Mastrandrea, J. Fournet, and A. Barrat. Contact patterns in a high school: A comparison between data collected using wearable sensors, contact diaries and friendship surveys. PLoS ONE, 2015.
  * [14] M. Salathé, M. Kazandjieva, J. W. Lee, P. Levis, M. W. Feldman, and J. H. Jones. A high resolution human contact network for infectious disease transmission. Proc Natl Acad Sci U S A, 107(51):22020–5, Dec 2010.
  * [15] Mossong, Joël, Niel Hens, Mark Jit, Philippe Beutels, Kari Auranen, Rafael Mikolajczyk, Marco Massari et al. ”Social contacts and mixing patterns relevant to the spread of infectious diseases.” PLoS Med 5, no. 3 (2008): e74.
  * [16] Cauchemez, Simon, Alain-Jacques Valleron, Pierre-Yves Boelle, Antoine Flahault, and Neil M. Ferguson. ”Estimating the impact of school closure on influenza transmission from Sentinel data.” Nature 452, no. 7188 (2008): 750-754.
  * [17] Eames, Ken TD, Natasha L. Tilston, Ellen Brooks-Pollock, and W. John Edmunds. ”Measured dynamic social contact patterns explain the spread of H1N1v influenza.” PLoS Comput Biol 8, no. 3 (2012): e1002425.
  * [18] https://www.mohfw.gov.in/covid_vaccination/vaccination/important-information.html
  * [19] Cao, Qing, Yi-Ching Chen, Chyi-Liang Chen, and Cheng-Hsun Chiu. ”SARS-CoV-2 infection in children: Transmission dynamics and clinical characteristics.” Journal of the Formosan Medical Association 119, no. 3 (2020): 670.
  * [20] Prem, Kiesha, Yang Liu, Timothy W. Russell, Adam J. Kucharski, Rosalind M. Eggo, Nicholas Davies, Stefan Flasche et al. ”The effect of control strategies to reduce social mixing on outcomes of the COVID-19 epidemic in Wuhan, China: a modelling study.” The Lancet Public Health 5, no. 5 (2020): e261-e270.
  * [21] Singh, Rajesh, and Ronojoy Adhikari. ”Age-structured impact of social distancing on the COVID-19 epidemic in India.” arXiv preprint arXiv:2003.12055 (2020).
  * [22] Wu, Joseph T., Kathy Leung, Mary Bushman, Nishant Kishore, Rene Niehus, Pablo M. de Salazar, Benjamin J. Cowling, Marc Lipsitch, and Gabriel M. Leung. ”Estimating clinical severity of COVID-19 from the transmission dynamics in Wuhan, China.” Nature medicine 26, no. 4 (2020): 506-510.
  * [23] To, Kelvin Kai-Wang, Owen Tak-Yin Tsang, Wai-Shing Leung, Anthony Raymond Tam, Tak-Chiu Wu, David Christopher Lung, Cyril Chik-Yan Yip et al. ”Temporal profiles of viral load in posterior oropharyngeal saliva samples and serum antibody responses during infection by SARS-CoV-2: an observational cohort study.” The Lancet Infectious Diseases 20, no. 5 (2020): 565-574.
  * [24] Zhou, Fei, Ting Yu, Ronghui Du, Guohui Fan, Ying Liu, Zhibo Liu, Jie Xiang et al. ”Clinical course and risk factors for mortality of adult inpatients with COVID-19 in Wuhan, China: a retrospective cohort study.” The lancet 395, no. 10229 (2020): 1054-1062.
  * [25] Verity, Robert, Lucy C. Okell, Ilaria Dorigatti, Peter Winskill, Charles Whittaker, Natsuko Imai, Gina Cuomo-Dannenburg et al. ”Estimates of the severity of coronavirus disease 2019: a model-based analysis.” The Lancet infectious diseases 20, no. 6 (2020): 669-677.
  * [26] Davies, Nicholas G., Petra Klepac, Yang Liu, Kiesha Prem, Mark Jit, and Rosalind M. Eggo. ”Age-dependent effects in the transmission and control of COVID-19 epidemics.” Nature medicine 26, no. 8 (2020): 1205-1211.
  * [27] Zou, Lirong, Feng Ruan, Mingxing Huang, Lijun Liang, Huitao Huang, Zhongsi Hong, Jianxiang Yu et al. ”SARS-CoV-2 viral load in upper respiratory specimens of infected patients.” New England Journal of Medicine 382, no. 12 (2020): 1177-1179.
  * [28] Jones, Terry C., Barbara Mühlemann, Talitha Veith, Guido Biele, Marta Zuchowski, Jörg Hoffmann, Angela Stein, Anke Edelmann, Victor Max Corman, and Christian Drosten. ”An analysis of SARS-CoV-2 viral load by patient age.” MedRxiv (2020).
  * [29] Ayoub, Houssein H., Hiam Chemaitelly, Shaheen Seedat, Ghina R. Mumtaz, Monia Makhoul, and Laith J. Abu-Raddad. ”Age could be driving variable SARS-CoV-2 epidemic trajectories worldwide.” PLoS One 15, no. 8 (2020): e0237959.
  * [30] Chikina, Maria, and Wesley Pegden. ”Modeling strict age-targeted mitigation strategies for COVID-19.” PloS one 15, no. 7 (2020): e0236237.
  * [31] Van den Driessche, Pauline, and James Watmough. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Mathematical biosciences 180.1-2 (2002): 29-48.
  * [32] Fiona M Guerra, Shelly Bolotin, Gillian Lim, Jane Heffernan, Shelley L Deeks, Ye Li, and Natasha S Crowcroft. The basic reproduction number $(R_{0})$ of measles: a systematic review. The Lancet Infectious Diseases, 17(12):e420- e428, 2017.
  * [33] Steven Sanche and Yen Ting Lin. High contagiousness and rapid spread of severe acute respiratory syndrome coronavirus. EID Journal, 26(7), 2020.
  * [34] https://en.wikipedia.org/wiki/Basic_reproduction_number#cite_note-cdc-who-2001-38
  * [35] Shuai, Zhisheng, and Pauline van den Driessche. ”Global stability of infectious disease models using Lyapunov functions.” SIAM Journal on Applied Mathematics 73.4 (2013): 1513-1532.
  * [36] LaSalle, J. P. The stability of dynamical systems, Regional Conf. Ser. Appl. Math. SIAM, Philadephia (1976).
  * [37] Freedman, Herb I., Shigui Ruan, and Moxun Tang. Uniform persistence and flows near a closed positively invariant set. Journal of Dynamics and Differential Equations 6.4 (1994): 583-600.
  * [38] Li, M. Y., Graef, J. R., Wang, L., and Karsai, J. Global dynamics of a SEIR model with varying total population size. Mathematical biosciences, 160(2)(1999), 191-213.
  * [39] Diprima BWE, ,Elementary RC. Differential equations and boundary value problems. New York: John Wiley Sons; 2009.
  * [40] Birkhoff, Garrett, and Gian-Carlo Rota. Ordinary differential equations. John Wiley Sons; 1989.
  * [41] Fleming WH, Rishel RW. Deterministic and stochastic optimal control. New York, NY, USA: Springer; 1975.
  * [42] Pontryagin LS, Boltyanskii VG, Gamkrelidze RV, and Mishchenko EF. The Mathematical Theory of Optimal Processes. New York, NY, USA: Wiley; 1962.
  * [43] Ministry of Home Affairs Unlock $4.0$ official guidelines on their official website (PDF). Ministry of Home Affairs.
  * [44] Ministry of Home Affairs Unlock $5.0$ official guidelines on their official website (PDF). Ministry of Home Affairs.
  * [45] “MHA Order Dt, 27.10.2020 reg extension of guidelines of reopening upto 30.11.2020” (PDF). Ministry of Home Affairs, GoI. 27 October 2020. Retrieved 12 November 2020.
  * [46] MHA DO Dt. 27.10.2020 to all Chief Secretaries reg extension of reopening upto 30.11.2020” (PDF). Ministry of Home Affairs, GoI. 27 October 2020. Retrieved 12 November2020.
  * [47] MHA DO Dt. 27.10.2020 to all Administrators reg extension of reopening upto 30.11.2020” (PDF). Ministry of Home Affairs, GoI. 27 October 2020. Retrieved 12 November 2020.
  * [48] Unlock 6.0 latest guidelines: What do they mean for schools, colleges, other educational institutions. DNA India. 29 October 2020. Retrieved 11 November 2020.
  * [49] Unlock 6.0 guidelines: Which states have allowed more relaxations in Nov?. The Indian Express. 6 November 2020. Retrieved 11 November 2020.
  * [50] https://www.populationpyramid.net
  * [51] Prem, Kiesha, Alex R. Cook, and Mark Jit. Projecting social contact matrices in 152 countries using contact surveys and demographic data. PLoS computational biology 13.9 (2017): e1005697.
  * [52] Wang, Xia, Hulin Wu, and Sanyi Tang. Assessing Age-Specific Vaccination Strategies and Post-Vaccination Reopening Policies for COVID-19 Control Using SEIR Modeling Approach. medRxiv (2021).
  * [53] Lauer, Stephen A., Kyra H. Grantz, Qifang Bi, Forrest K. Jones, Qulu Zheng, Hannah R. Meredith, Andrew S. Azman, Nicholas G. Reich, and Justin Lessler. The incubation period of coronavirus disease 2019 (COVID-19) from publicly reported confirmed cases: estimation and application. Annals of internal medicine 172, no. 9 (2020): 577-582.
  * [54] https://ourworldindata.org/covid-cases?country=IND~USA~GBR~CAN~DEU~FRA
  * [55] Chan, Jasper Fuk-Woo, Shuofeng Yuan, Kin-Hang Kok, Kelvin Kai-Wang To, Hin Chu, Jin Yang, Fanfan Xing et al. A familial cluster of pneumonia associated with the 2019 novel coronavirus indicating person-to-person transmission: a study of a family cluster. The lancet 395, no. 10223 (2020): 514-523.
  * [56] Lancet, The. India under COVID-19 lockdown. Lancet (London, England) 395, no. 10233 (2020): 1315.
  * [57] Chatterjee, Patralekha. Is India missing COVID-19 deaths?. The Lancet 396, no. 10252 (2020): 657.
  * [58] Wallinga, Jacco, Peter Teunis, and Mirjam Kretzschmar. Using data on social contacts to estimate age-specific transmission parameters for respiratory-spread infectious agents. American journal of epidemiology 164, no. 10 (2006): 936-944.
  * [59] Read, Jonathan M., and Matt J. Keeling. Disease evolution on networks: the role of contact structure. Proceedings of the Royal Society of London. Series B: Biological Sciences 270, no. 1516 (2003): 699-708.
  * [60] Edmunds, W. John, C. J. O’callaghan, and D. J. Nokes. Who mixes with whom? A method to determine the contact patterns of adults that may lead to the spread of airborne infections. Proceedings of the Royal Society of London. Series B: Biological Sciences 264, no. 1384 (1997): 949-957.
  * [61] Stage, Helena B., Joseph Shingleton, Sanmitra Ghosh, Francesca Scarabel, Lorenzo Pellis, and Thomas Finnie. Shut and re-open: the role of schools in the spread of COVID-19 in Europe. Philosophical Transactions of the Royal Society B 376, no. 1829 (2021): 20200277.
  * [62] Monod, Mélodie, Alexandra Blenkinsop, Xiaoyue Xi, Daniel Hebert, Sivan Bershan, Simon Tietze, Marc Baguelin et al. Age groups that sustain resurging COVID-19 epidemics in the United States. Science 371, no. 6536 (2021): eabe8372.

## Appendix A

Proof of Theorem 3.1

###### Proof.

First, we show that $S_{i}(t)>0$ for all $t>0.$ For this, we assume that there
exists a $t^{s}_{i}>0$ such that, $S_{i}(t^{s}_{i})=0$ and $S_{i}(t)>0$ for
$0<t<t^{s}_{i}.$ Thus, the first equation of system (3) becomes

$\frac{dS_{i}(t^{s}_{i})}{dt}=\Gamma_{i}-\beta
S_{i}(t^{s}_{i})\sum_{j=1}^{n}\lambda_{ij}\left(\frac{I^{a}_{j}(t^{s}_{i})+\eta_{j}I^{s}_{j}(t^{s}_{i})}{N_{j}}\right)-\mu_{i}S_{i}(t^{s}_{i})=\Gamma_{i}>0.$

Therefore, $S_{i}(t)<0$ for
$t\in\left(t^{s}_{i}-\epsilon^{s},t^{s}_{i}\right)$ and sufficiently small
positive $\epsilon^{s}.$ This is a contradiction to $S_{i}(t)>0$ for
$0<t<t^{s}_{i}.$ Hence, $S_{i}(t)>0$ for all $t>0,$ where $i=1,2,\cdots,n.$
Next, we prove that $E_{i}(t)>0$ for $i=1,2,\cdots,n$ and all $t>0.$ By theory
of differential equations, we obtain that the solution of second equation of
system (3) is given by

$E_{i}(t)\geq E_{i}(0)e^{-(\alpha_{i}+\mu_{i})t}$

Therefore, we have that $E_{i}(t)\geq 0$ for $t\geq 0.$ In the similiar
manner, it can easily be proved that $I^{a}_{i}(t)\geq 0$, $I^{s}_{i}(t)\geq
0$ and $R_{i}(t)\geq 0$ for $t\geq 0.$
Further, we show the boundedness of solutions of model system (3). From model
system (3), we have

$\frac{d}{dt}(S_{i}+E_{i}+I_{i}^{a}+I_{i}^{s}+R_{i})=\lambda_{i}-\mu_{i}(S_{i}+E_{i}+I_{i}^{a}+I_{i}^{s}+R_{i})-\delta^{a}_{i}I_{i}^{a}-\delta^{s}_{i}I_{i}^{s}$
$\frac{dN_{i}}{dt}\leq\Lambda_{i}-\mu_{i}N_{i}\,\,\,\,\text{for}\,\,i^{th}group.$

Sinec, we obtain
$\lim_{t\to\infty}supN_{i}(t)\leq\frac{\Lambda_{i}}{\mu_{i}}$. Hence, we have

$\Omega=\left\\{\left(S_{i},\,E_{i},\,I^{a}_{i},\,I^{s}_{i},\,R_{i}\right)\in\mathbb{R}_{+}^{5}|0<S_{i}+E_{i}+I^{a}_{i}+I^{s}_{i}+R_{i}\leq\frac{\Lambda_{i}}{\mu_{i}}\right\\}$

Thus, the solutions of model system (3) is bounded. By the proof procedure of
Theorem 3.1, we know that all solutions of system (3) ultimately come in and
continue in the region $\Omega.$ Hence $\Omega$ is a bounded absorbing set for
system (3). Thus, this completes the proof of Theorem 3.1. ∎

### Existence of Solutions

######  Theorem 6.1.

The model system (3) has a unique solution when initial conditions for all the
variables are nonnegative i.e., $S_{i}(0)\geq 0,E_{i}(0)\geq
0,I^{a}_{i}(0)\geq 0,I^{s}_{i}(0)\geq 0,R_{i}(0)\geq 0$ for all
$i=1,2,\cdots,n.$

###### Proof.

Let $X=\left(\begin{array}[]{c}X_{1}\\\ X_{2}\\\ \vdots\\\
X_{n}\end{array}\right)$ where $X_{i}=\left(\begin{array}[]{c}S_{i}(t)\\\
E_{i}(t)\\\ I^{a}_{i}(t)\\\ I^{s}_{i}(t)\\\ R_{i}(t)\end{array}\right)$ for
all $i=1,2,\cdots,n.$ So, the model system (3) can be written in the following
form $\Phi(X)=AX+B(X),$ where

$A=\left(\begin{array}[]{cccc}A_{1}&0&\cdots&0\\\ 0&A_{2}&\cdots&0\\\
\vdots&\vdots&\vdots&\vdots\\\ 0&0&\cdots&A_{n}\end{array}\right)$

with

$A_{i}=\left(\begin{array}[]{ccccc}-\mu_{i}&0&0&0&0\\\
0&-(\alpha_{i}+\mu_{i})&0&0&0\\\
0&0&-(\gamma^{a}_{i}+\delta^{a}_{i}+\mu_{i})&0&0\\\
0&0&0&-(\gamma^{s}_{i}+\delta^{s}_{i}+\mu_{i})&0\\\
0&0&0&0&-\mu_{i}\end{array}\right)$

and $B(X)=\left(\begin{array}[]{c}B(X_{1})\\\ B(X_{2})\\\ \vdots\\\
B(X_{n})\end{array}\right)$ with
$B(X_{i})=\left(\begin{array}[]{c}\Gamma_{i}-\beta
S_{i}(t)\sum_{n}^{j=1}\lambda_{ij}\frac{I^{a}_{j}+\eta_{j}I^{s}_{j}}{N_{j}}\\\
\beta
S_{i}(t)\sum_{n}^{j=1}\lambda_{ij}\frac{I^{a}_{j}+\eta_{j}I^{s}_{j}}{N_{j}}\\\
0\\\ 0\\\ 0\end{array}\right).$ The functions $B(X_{i})$ satisfies

$\displaystyle\left\lvert B(X^{1}_{i})-B(X^{2}_{i})\right\rvert=$
$\displaystyle\left\lvert\beta
S^{1}_{i}(t)\sum_{n}^{j=1}\lambda_{ij}\frac{I^{a1}_{j}+\eta_{j}I^{s1}_{j}}{N_{j}}-\beta
S^{2}_{i}(t)\sum_{n}^{j=1}\lambda_{ij}\frac{I^{a2}_{j}+\eta_{j}I^{s2}_{j}}{N_{j}}\right\rvert$
$\displaystyle\leq$ $\displaystyle\left\lvert\beta
S^{1}_{i}(t)\sum_{n}^{j=1}\lambda_{ij}\frac{I^{a1}_{j}+\eta
I^{s1}_{j}}{N_{j}}-\beta
S^{1}_{i}(t)\sum_{n}^{j=1}\lambda_{ij}\frac{I^{a2}_{j}+\eta_{j}I^{s2}_{j}}{N_{j}}\right\rvert$
$\displaystyle+\left\lvert\beta
S^{1}_{i}(t)\sum_{n}^{j=1}\lambda_{ij}\frac{I^{a2}_{j}+\eta_{j}I^{s2}_{j}}{N_{j}}-\beta
S^{2}_{i}(t)\sum_{n}^{j=1}\lambda_{ij}\frac{I^{a2}_{j}+\eta_{j}I^{s2}_{j}}{N_{j}}\right\rvert$
$\displaystyle=$ $\displaystyle\left\lvert
S^{1}_{i}(t)\right\rvert\left\lvert\beta\sum_{n}^{j=1}\frac{\lambda_{ij}}{N_{j}}\left(\left(I^{a1}_{j}-I^{a2}_{j}\right)+\eta\left(I^{s1}_{j}-I^{s2}_{j}\right)\right)\right\rvert$
$\displaystyle+\left\lvert\beta\sum_{n}^{j=1}\lambda_{ij}\frac{I^{a2}_{j}+\eta_{j}I^{s2}_{j}}{N_{j}}\left(S^{1}_{i}(t)-S^{2}_{i}(t)\right)\right\rvert.$
$\displaystyle\left\lvert B(X^{1}_{i})-B(X^{2}_{i})\right\rvert\leq$
$\displaystyle\frac{\Gamma_{i}}{\mu_{i}}\left(\beta\sum_{n}^{j=1}\frac{\lambda_{ij}}{N_{j}}\left(\left\lvert
I^{a1}_{j}-I^{a2}_{j}\right\rvert+\eta_{j}\left\lvert
I^{s1}_{j}-I^{s2}_{j}\right\rvert\right)\right)+\beta(1+\eta_{j})\sum_{j=1}^{n}\frac{\lambda_{ij}}{N_{j}}\frac{\Gamma_{j}}{\mu_{j}}\left\lvert
S^{1}_{i}(t)-S^{2}_{i}(t)\right\rvert.$

we also have that

$\left\lvert B(X^{1})-B(X^{2})\right\rvert\leq\sum_{i=1}^{n}\left\lvert
B(X^{1}_{i})-B(X^{2}_{i})\right\rvert.$

Therefore, we obtain

$\left\lvert B(X^{1})-B(X^{2})\right\rvert\leq
M\sum_{i=1}^{n}\left(\left\lvert I^{a1}_{i}-I^{a2}_{i}\right\rvert+\left\lvert
I^{s1}_{i}-I^{s2}_{i}\right\rvert+\left\lvert
S^{1}_{i}(t)-S^{2}_{i}(t)\right\rvert\right)$ $\implies\left\lvert
B(X^{1})-B(X^{2})\right\rvert\leq M\lvert\left\lvert
X_{1}-X_{2}\right\rvert\rvert,$

where

$M=\beta(1+\eta_{j})\sum_{i=1}^{n}\sum_{j=1}^{n}\frac{\Gamma_{i}}{\mu_{i}}\frac{\lambda_{ij}}{N_{j}}.$

Therefore, we have that

$\lvert\lvert\Phi(X_{1})-\Phi(X_{2})\rvert\rvert\leq M^{0}\lvert\left\lvert
X_{1}-X_{2}\right\rvert\rvert,$

where $M^{0}=\min\left(M,\left\lvert\left\lvert
A\right\rvert\right\rvert\right).$ Therefore, it follows that the function
$\Phi$ is uniformly Lipschitz continuous, and the condition on $S_{i}(0)\geq
0,E_{i}(0)\geq 0,I^{a}_{i}(0)\geq 0,I^{s}_{i}(0)\geq 0,R_{i}(0)\geq 0$ for all
$i=1,2,\cdots,n.$ Hence, using results stated in [40] (section 10 of the first
chapter), we conclude that a solution of the system (3) exists. ∎

Basic Reproduction Number $(R_{0}):$ We rewrite the middle $3n$ (second, third
and fourth) equations of system (3) as follows:
$x^{\prime}=\mathcal{F}-\mathcal{V}$ and
$x=\left(E_{1},E_{2},\cdots,E_{n},I^{a}_{1},I^{a}_{2},\cdots,I^{a}_{n},I^{s}_{1},I^{s}_{2},\cdots,I^{s}_{n}\right)\in\mathbb{R}^{3n},$
where

$\mathcal{F}=\left(\begin{array}[]{c}\beta
S_{1}(t)\sum_{j=1}^{n}\lambda_{1j}\left(\frac{I^{a}_{j}(t)+\eta_{j}I^{s}_{j}(t)}{N_{j}(t)}\right)\\\
\beta
S_{2}(t)\sum_{j=1}^{n}\lambda_{2j}\left(\frac{I^{a}_{j}(t)+\eta_{j}I^{s}_{j}(t)}{N_{j}(t)}\right)\\\
\vdots\\\ \beta
S_{n}(t)\sum_{j=1}^{n}\lambda_{nj}\left(\frac{I^{a}_{j}(t)+\eta_{j}I^{s}_{j}(t)}{N_{j}(t)}\right)\\\
0\\\ 0\\\ \vdots\\\ 0\\\ 0\\\ 0\\\ \vdots\\\
0\end{array}\right)\,\,\text{and}\,\mathcal{V}=\left(\begin{array}[]{c}(\alpha_{1}+\mu_{1})E_{1}(t)\\\
(\alpha_{2}+\mu_{2})E_{2}(t)\\\ \vdots\\\ (\alpha_{2}+\mu_{2})E_{2}(t)\\\
(1-\sigma)\alpha_{1}E_{1}(t)-\left(\gamma^{a}_{1}+\delta^{a}_{1}+\mu_{1}\right)I^{a}_{1}(t)\\\
(1-\sigma)\alpha_{2}E_{2}(t)-\left(\gamma^{a}_{2}+\delta^{a}_{2}+\mu_{2}\right)I^{a}_{2}(t)\\\
\vdots\\\
(1-\sigma)\alpha_{n}E_{n}(t)-\left(\gamma^{a}_{n}+\delta^{a}_{n}+\mu_{n}\right)I^{a}_{i}(t)\\\
\sigma\alpha_{1}E_{1}(t)-\left(\gamma^{s}_{1}+\delta^{s}_{1}+\mu_{1}\right)I^{s}_{1}(t)\\\
\sigma\alpha_{2}E_{2}(t)-\left(\gamma^{s}_{2}+\delta^{s}_{2}+\mu_{2}\right)I^{s}_{2}(t)\\\
\vdots\\\
\sigma\alpha_{n}E_{n}(t)-\left(\gamma^{s}_{n}+\delta^{s}_{n}+\mu_{n}\right)I^{s}_{n}(t)\end{array}\right).$
(11)

Further, by calculating the Jacobian matrices $F$ and $V$ at the DFE
$(E^{0}),$ we obtain

$F=\left(\begin{array}[]{ccc}\bf{O}&F_{1}&F_{2}\\\ \bf{O}&\bf{O}&\bf{O}\\\
\bf{O}&\bf{O}&\bf{O}\end{array}\right)_{3n\times
3n}\,\,\text{and}\,V=\left(\begin{array}[]{ccc}V_{11}&\bf{O}&\bf{O}\\\
-V_{21}&V_{22}&\bf{O}\\\ -V_{31}&\bf{O}&V_{33}\end{array}\right)_{3n\times
3n},$ (12)

where $\bf{O}$ is the zero matrix of order $n\times n$ having all entries
equal to zero and

$F_{1}=\left(\begin{array}[]{cccc}\beta\lambda_{11}&\beta\lambda_{12}&\cdots&\beta\lambda_{1n}\\\
\beta\lambda_{21}&\beta\lambda_{22}&\cdots&\beta\lambda_{nn}\\\
\vdots&\vdots&&\vdots\\\
\beta\lambda_{n1}&\beta\lambda_{n2}&\cdots&\beta\lambda_{nn}\end{array}\right)_{n\times
n},\,\,V_{22}=\left(\begin{array}[]{cccc}\gamma^{a}_{1}+\delta^{a}_{1}+\mu_{1}&0&\cdots&0\\\
0&\gamma^{a}_{2}+\delta^{a}_{2}+\mu_{2}&\cdots&0\\\ \vdots&\vdots&&\vdots\\\
0&0&\cdots&\gamma^{a}_{n}+\delta^{a}_{n}+\mu_{n}\end{array}\right)_{n\times
n},$ $V_{11}=\left(\begin{array}[]{cccc}\alpha_{1}+\mu_{1}&0&\cdots&0\\\
0&\alpha_{2}+\mu_{2}&\cdots&0\\\ \vdots&\vdots&&\vdots\\\
0&0&\cdots&\alpha_{n}+\mu_{n}\end{array}\right)_{n\times
n},\,V_{21}=\left(\begin{array}[]{cccc}(1-\sigma)\alpha_{1}&0&\cdots&0\\\
0&(1-\sigma)\alpha_{2}&\cdots&0\\\ \vdots&\vdots&&\vdots\\\
0&0&\cdots&(1-\sigma)\alpha_{n}\end{array}\right)_{n\times n},$
$V_{31}=\left(\begin{array}[]{cccc}\sigma\alpha_{1}&0&\cdots&0\\\
0&(\sigma\alpha_{2}&\cdots&0\\\ \vdots&\vdots&&\vdots\\\
0&0&\cdots&\sigma\alpha_{n}\end{array}\right)_{n\times
n},\,\text{and}\,V_{33}=\left(\begin{array}[]{cccc}\gamma^{s}_{1}+\delta^{s}_{1}+\mu_{1}&0&\cdots&0\\\
0&\gamma^{s}_{2}+\delta^{s}_{2}+\mu_{2}&\cdots&0\\\ \vdots&\vdots&&\vdots\\\
0&0&\cdots&\gamma^{s}_{n}+\delta^{s}_{n}+\mu_{n}\end{array}\right)_{n\times
n},$
$F_{2}=\left(\begin{array}[]{cccc}\beta\eta_{1}\lambda_{11}&\beta\eta_{2}\lambda_{12}&\cdots&\beta\eta_{n}\lambda_{1n}\\\
\beta\eta_{1}\lambda_{21}&\beta\eta_{2}\lambda_{22}&\cdots&\beta\eta_{n}\lambda_{2n}\\\
\vdots&\vdots&&\vdots\\\
\beta\eta_{1}\lambda_{n1}&\beta\eta_{2}\lambda_{n2}&\cdots&\beta\eta_{n}\lambda_{nn}\end{array}\right)_{n\times
n},$

Proof for Theorem 3.2

###### Proof.

Let
$x=\left(E_{1},E_{2},\cdots,E_{n},I^{a}_{1},I^{a}_{2},\cdots,I^{a}_{n},I^{s}_{1},I^{s}_{2},\cdots,I^{s}_{n}\right)\in\mathbb{R}^{3n}$
and
$y=(S_{1},S_{2},\cdots,S_{n},R_{1},R_{2},\cdots,R_{n})\in\mathbb{R}^{2n}$ be
the disease compartment and disease-free compartment vector, respectively. We
set $f(x,y)=(F-V)x-\mathcal{F}(x,y)+\mathcal{V}(x,y)$ where vector function
$F,\,V$ are defined in Eq. (12) and vector function
$\mathcal{F}(x,y),\,\mathcal{V}(x,y)$ are defined in Eq. (11). Further, for
the disease compartments and disease free compartments, the model system (3)
can be written as

$x^{\prime}=(F-V)x-f(x,y)\,\,\text{and}\,y^{\prime}=g(x,y),$ (13)

where

$f(x,y)=\left(\begin{array}[]{c}\beta\sum_{j=1}^{n}\lambda_{1j}\left(1-\frac{S_{1}}{N_{j}}\right)\left(I^{a}_{j}+\eta_{j}I^{s}_{j}\right)\\\
\beta\sum_{j=1}^{n}\lambda_{2j}\left(1-\frac{S_{2}}{N_{j}}\right)\left(I^{a}_{j}+\eta_{j}I^{s}_{j}\right)\\\
\vdots\\\
\beta\sum_{j=1}^{n}\lambda_{nj}\left(1-\frac{S_{n}}{N_{j}}\right)\left(I^{a}_{j}+\eta_{j}I^{s}_{j}\right)\\\
0\\\ 0\\\ \vdots\\\ 0\\\ 0\\\ 0\\\ \vdots\\\
0\end{array}\right)\,\,\,\text{with}\,\,f(0,y)=0.$

It can be easily observed that $f(x,y)\geq 0$ in
$\Omega\subset\mathbb{R}^{3n+2n}_{+}\,,F$ and $V$ are nonnegative matrices.
$P_{0}=(0,y^{0})$ be the disease free equilibrium for the system (13) that is
equivalent to DFE $E^{0}$ of system (3). Since $y^{\prime}=g(0,y)$ has a
unique positive equilibrium
$y^{0}=\left(S_{1}^{0},S^{0}_{2},\cdots,S^{0}_{n},0,0,\cdots,0\right)\in\mathbb{R}^{2n},$
then it is globally asymptotically stable in $\mathbb{R}^{2n}.$ Since
$\lambda_{ij}$ is irreducible, therefore $V^{-1}F$ is also irreducible and
nonnegative. It follows by Perron-Frobenius theory that $c_{T}$ be the
nonnegative left eigenvector of the matrix $V^{-1}F$ corresponding to the
eigenvalue $R_{0}=\rho(V^{-1}F).$ Further, we assume $L_{DFE}=c^{T}V^{-1}x$ is
a Lyapunov function for the system (13) on the region $\Omega.$ By
differentiating $L_{DFE}$ along the solution of the system (13), we obtain

$\displaystyle L^{\prime}_{DFE}=$ $\displaystyle
c^{T}V^{-1}x^{\prime}=c^{T}V^{-1}\left((F-V)x-f(x,y)\right)=c^{T}V^{-1}(F-V)x-c^{T}V^{-1}f(x,y)$
(14) $\displaystyle=$ $\displaystyle(R_{0}-1)c^{T}x-c^{T}V^{-1}f(x,y).$

Thus, if $R_{0}\leq 0,$ then $L^{\prime}_{DFE}\leq 0$ in $\Omega$ and
$L^{\prime}_{DFE}=0$ implies that $c^{T}x=0$ hence $x=0.$ Using the global
stability for disease free system $y^{\prime}=g(0,y)$ and $f(0,y)=0,$
singleton $\\{P_{0}\\}$ is the only invariant set in $\mathbb{R}^{5n}_{+}$
where $L^{\prime}_{DFE}=0.$ When $R_{0}=1,$ singleton $\\{P_{0}\\}$ is largest
invariant set where $L^{\prime}_{DFE}=c^{T}V^{-1}f(x,y)=0.$ Therefore, by
LaSalle’s invariance principle [36], $P_{0}$ is globally asymptotically stable
in $\Omega$ when $R_{0}\leq 1.$
If $R_{0}>1$ then $L^{\prime}_{DFE}>0$ on condition that $x>0$ and $y=y^{0}.$
By continuity of $L^{\prime}_{DFE}>0$ in the interval
$(P_{0}-\epsilon,P_{0}+\epsilon)$ where $\epsilon$ is a small number, one can
say that all the solutions of system (13) in the positive cone and near by
$P_{0}$ diverges from $P_{0}$, except those on the invariant y-axis. This
implies that $P_{0}$ is unstable. Using the uniform persistence result form
[37] and an argument given in the proof of Proposition 3.3 of [38], it can be
shown that when $R_{0}>1,$ the uniform persistence of system (3) is assured by
instability of $P_{0}.$ Further, the existence of at least one EE can be
confirmed by using the uniform persistence and the concept of positive
invariance of compact set $\Omega.$ ∎

Proof for Theorem 3.3

###### Proof.

To determine the global stability of endemic equilibrium $E^{*}$ of system
(3), let
$D_{1i}=S_{i}-\bar{S_{i}}-\bar{S_{i}}\ln\frac{S_{i}}{\bar{S_{i}}}+E_{i}-\bar{E_{i}}-\bar{E_{i}}\ln\frac{E_{i}}{\bar{E_{i}}},$
$D_{2i}=I^{a}_{i}-\bar{I^{a}_{i}}-\bar{I^{a}_{i}}\ln\frac{I^{a}_{i}}{\bar{I^{a}_{i}}}+I^{s}_{i}-\bar{I^{s}_{i}}-\bar{I^{s}_{i}}\ln\frac{I^{s}_{i}}{\bar{I^{s}_{i}}}$
and $D_{3i}=\bar{R_{i}}-\bar{R_{i}}\ln\frac{R_{i}}{\bar{R_{i}}}.$ By
differentiating $D_{1i}$ along the solutions of system (3), we obtain

$\displaystyle D^{\prime}_{1i}=$
$\displaystyle\left(1-\frac{\bar{S_{i}}}{S_{i}}\right)S^{\prime}_{i}+\left(1-\frac{\bar{E_{i}}}{E_{i}}\right)E^{\prime}_{i}$
$\displaystyle=$
$\displaystyle\left(1-\frac{\bar{S_{i}}}{S_{i}}\right)\left(\Gamma_{i}-\beta
S_{i}\sum_{j=1}^{n}\lambda_{ij}\left(\frac{I^{a}_{j}+\eta_{j}I^{s}_{j}}{N_{j}}\right)-\mu_{i}S_{i}\right)$
$\displaystyle+\left(1-\frac{\bar{E_{i}}}{E_{i}}\right)\left(\beta
S_{i}\sum_{j=1}^{n}\lambda_{ij}\left(\frac{I^{a}_{j}+\eta_{j}I^{s}_{j}}{N_{j}}\right)-(\alpha_{i}+\mu_{i})E_{i}\right).$
By using the system of nonlinear equations (5), we obtain $\displaystyle=$
$\displaystyle\left(1-\frac{\bar{S_{i}}}{S_{i}}\right)\left(\beta\bar{S_{i}}\sum_{j=1}^{n}\lambda_{ij}\left(\frac{\bar{I^{a}_{j}}+\eta_{j}\bar{I^{s}_{j}}}{N_{j}}\right)+\mu_{i}\bar{S_{i}}-\beta
S_{i}\sum_{j=1}^{n}\lambda_{ij}\left(\frac{I^{a}_{j}+\eta_{j}I^{s}_{j}}{N_{j}}\right)-\mu_{i}S_{i}\right)$
$\displaystyle+\left(1-\frac{\bar{E_{i}}}{E_{i}}\right)\left(\beta
S_{i}\sum_{j=1}^{n}\lambda_{ij}\left(\frac{I^{a}_{j}+\eta_{j}I^{s}_{j}}{N_{j}}\right)-\beta\bar{S_{i}}\sum_{j=1}^{n}\lambda_{ij}\left(\frac{\bar{I^{a}_{j}}+\eta_{j}\bar{I^{s}_{j}}}{N_{j}}\right)\right)$
$\displaystyle=$
$\displaystyle\mu_{i}\bar{S_{i}}\left(2-\frac{\bar{S_{i}}}{S_{i}}-\frac{S_{i}}{\bar{S_{i}}}\right)+2\beta\bar{S_{i}}\sum_{j=1}^{n}\lambda_{ij}\left(\frac{\bar{I^{a}_{j}}+\eta_{j}\bar{I^{s}_{j}}}{\bar{N_{j}}}\right)-\beta\bar{S_{i}}\sum_{j=1}^{n}\lambda_{ij}\left(\frac{\bar{I^{a}_{j}}+\eta_{j}\bar{I^{s}_{j}}}{\bar{N_{j}}}\right)\frac{\bar{S_{i}}}{S_{i}}$
$\displaystyle+\beta\bar{S_{i}}\sum_{j=1}^{n}\lambda_{ij}\left(\frac{I^{a}_{j}+\eta_{j}I^{s}_{j}}{\bar{N_{j}}}\right)-\beta\bar{S_{i}}\sum_{j=1}^{n}\lambda_{ij}\left(\frac{\bar{I^{a}_{j}}+\eta_{j}\bar{I^{s}_{j}}}{\bar{N_{j}}}\right)\frac{E_{i}}{\bar{E_{i}}}-\beta
S_{i}\sum_{j=1}^{n}\lambda_{ij}\left(\frac{I^{a}_{j}+\eta_{j}I^{s}_{j}}{\bar{N_{j}}}\right)\frac{\bar{E_{i}}}{E_{i}}$
$\displaystyle\leq$ $\displaystyle
2\beta\bar{S_{i}}\sum_{j=1}^{n}\lambda_{ij}\left(\bar{I^{a}_{j}}+\eta_{j}\bar{I^{s}_{j}}\right)-\beta\bar{S_{i}}\sum_{j=1}^{n}\lambda_{ij}\left(I^{a}_{j}+\eta_{j}I^{s}_{j}\right)\frac{\bar{S_{i}}}{S_{i}}+\beta\bar{S_{i}}\sum_{j=1}^{n}\lambda_{ij}\left(I^{a}_{j}+\eta_{j}I^{s}_{j}\right)$
$\displaystyle-\beta\bar{S_{i}}\sum_{j=1}^{n}\lambda_{ij}\left(I^{a}_{j}+\eta_{j}I^{s}_{j}\right)\frac{E_{i}}{\bar{E_{i}}}-\beta
S_{i}\sum_{j=1}^{n}\lambda_{ij}\left(I^{a}_{j}+\eta_{j}I^{s}_{j}\right)\frac{\bar{E_{i}}}{E_{i}}$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{n}\beta\lambda_{ij}\bar{S_{i}}\bar{I^{a}_{j}}\left(2-\frac{\bar{S_{i}}}{S_{i}}-\frac{E_{i}}{\bar{E_{i}}}+\frac{I^{a}_{j}}{\bar{I^{a}_{j}}}+\frac{S_{i}I^{a}_{j}\bar{E_{i}}}{\bar{S_{i}}\bar{I^{a}_{j}}E_{i}}\right)+\sum_{j=1}^{n}\beta\eta_{j}\lambda_{ij}\bar{S_{i}}\bar{I^{s}_{j}}\left(2-\frac{\bar{S_{i}}}{S_{i}}-\frac{E_{i}}{\bar{E_{i}}}+\frac{I^{s}_{j}}{\bar{I^{s}_{j}}}+\frac{S_{i}I^{s}_{j}\bar{E_{i}}}{\bar{S_{i}}\bar{I^{s}_{j}}E_{i}}\right).$
$\displaystyle\text{Using the inequality}\,1-x\leq\log x,\text{we obtain}$
$\displaystyle\leq$
$\displaystyle\sum_{j=1}^{n}\beta\lambda_{ij}\bar{S_{i}}\bar{I^{a}_{j}}\left(\frac{I^{a}_{j}}{\bar{I^{a}_{j}}}-\frac{E_{i}}{\bar{E_{i}}}-\log\frac{I^{a}_{j}}{\bar{I^{a}_{j}}}+\log\frac{E_{i}}{\bar{E_{i}}}\right)+\sum_{j=1}^{n}\beta\eta_{j}\lambda_{ij}\bar{S_{i}}\bar{I^{s}_{j}}\left(\frac{I^{s}_{j}}{\bar{I^{s}_{j}}}-\frac{E_{i}}{\bar{E_{i}}}-\log\frac{I^{s}_{j}}{\bar{I^{s}_{j}}}+\log\frac{E_{i}}{\bar{E_{i}}}\right)$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{2n}a_{ij}G_{ij}\,\,\,\text{with}\,\,\,\,a_{ij}=\left\\{\begin{split}&\beta\bar{S_{i}}\bar{I^{a}_{j}}\qquad
1\leq j\leq n,\\\ &\beta\eta_{j}\bar{S_{i}}\bar{I^{s}_{j}}\qquad n+1\leq j\leq
2n.\end{split}\right.$

By differentiating $D_{2i}$ along the solutions of system (3), we obtain

$\displaystyle D_{n+i}^{\prime}=$
$\displaystyle\left(1-\frac{\bar{I^{a}_{i}}}{I^{a}_{i}}\right)(I^{a}_{i})^{\prime}+\left(1-\frac{\bar{I^{s}_{i}}}{I^{s}_{i}}\right)(I^{s}_{i})^{\prime}$
$\displaystyle=$
$\displaystyle\left(1-\frac{\bar{I^{a}_{i}}}{I^{a}_{i}}\right)\left((1-\sigma)\alpha_{i}E_{i}-\left(\gamma^{a}_{i}+\delta^{a}_{i}+\mu_{i}\right)I^{a}_{i}\right)+\left(1-\frac{\bar{I^{s}_{i}}}{I^{s}_{i}}\right)\left(\sigma\alpha_{i}E_{i}-\left(\gamma^{s}_{i}+\delta^{s}_{i}+\mu_{i}\right)I^{s}_{i}\right)$
$\displaystyle=$
$\displaystyle\left(1-\frac{\bar{I^{a}_{i}}}{I^{a}_{i}}\right)\left((1-\sigma)\alpha_{i}E_{i}-(1-\sigma)\alpha_{i}\bar{E_{i}}\frac{I^{a}_{i}}{\bar{I^{a}_{i}}}\right)+\left(1-\frac{\bar{I^{s}_{i}}}{I^{s}_{i}}\right)\left(\sigma\alpha_{i}E_{i}-\sigma\alpha_{i}\bar{E_{i}}\frac{I^{s}_{i}}{\bar{I^{s}_{i}}}\right)$
$\displaystyle=$
$\displaystyle(1-\sigma)\alpha_{i}\bar{E_{i}}\left(1-\frac{\bar{I^{a}_{i}}}{I^{a}_{i}}\right)\left(\frac{E_{i}}{\bar{E_{i}}}-\frac{I^{a}_{i}}{\bar{I^{a}_{i}}}\right)+\sigma\alpha_{i}\bar{E_{i}}\left(1-\frac{\bar{I^{s}_{i}}}{I^{s}_{i}}\right)\left(\frac{E_{i}}{\bar{E_{i}}}-\frac{I^{s}_{i}}{\bar{I^{s}_{i}}}\right)$
$\displaystyle\leq$
$\displaystyle(1-\sigma)\alpha_{i}\bar{E_{i}}\left(\frac{E_{i}}{\bar{E_{i}}}-\frac{I^{a}_{i}}{\bar{I^{a}_{i}}}-\log\frac{E_{i}}{\bar{E_{i}}}+\log\frac{I^{a}_{i}}{\bar{I^{a}_{i}}}\right)+\sigma\alpha_{i}\bar{E_{i}}\left(\frac{E_{i}}{\bar{E_{i}}}-\frac{I^{s}_{i}}{\bar{I^{s}_{i}}}-\log\frac{E_{i}}{\bar{E_{i}}}+\log\frac{I^{s}_{i}}{\bar{I^{s}_{i}}}\right)$
$\displaystyle=$ $\displaystyle
a_{n+i,i}G_{n+i,i}\,\,\,\text{with}\,\,\,\,a_{n+i,i}=\left\\{\begin{split}&(1-\sigma)\alpha_{i}\bar{E_{i}}\qquad&1\leq
i\leq n,\\\ &\sigma\alpha_{i}\bar{E_{i}}\qquad&n+1\leq i\leq
2n.\end{split}\right.$

$n+i$$i$$2n+i$$j$$n+j$$2n+j$$a_{i,2n+i}$$a_{2n+i,i}$$a_{j,2n+j}$$a_{2n+j,j}$$a_{n+i,i}$$a_{i,n+i}$$a_{n+j,j}$$a_{j,n+j}$$a_{i,n+j}$$a_{j,n+i}$$a_{i,2n+j}$$a_{j,2n+i}$
Figure 8: The weighted digraph $(\mathcal{G},\mathcal{A})$ constructed for the
model system (3) with two group

Hence assumption (1) of Theorem 3.5 stated in [35] holds. To verify assumption
(2) of the theorem 3.5 stated in [35], we define the weighted digraph
$(\mathcal{G},\mathcal{A})$associated with the weight matrix
$A=\left[a_{ij}\right]$ with $a_{ij}>0$ as defined above and all other
$a_{ij}=0.$ In digraph $(\mathcal{G},\mathcal{A})$, there are two kinds of
cycles involving direct transmission and cycles involving indirect
transmission. For each cycle, assumption (2) of Theorem 3.5 can be verified.
Therefore, by Theorem 3.5 discussed in [35], we have that
$D=\sum_{i=1}^{n}c_{i}D_{i}$ is a Lyapunov function for model system (3).
Since, $d^{-}(n+i)=d^{-}(2n+i)=1$ for each $i$ by Theorem 3.3 and 3.4 stated
in [35], we have $c_{n+i}=\sum_{j=1}^{n}c_{j}\frac{a_{j,n+i}}{a_{n+i,i}}.$
Thus, we obtain

$D=\sum_{i=1}^{n}c_{i}D_{i}+\sum_{i=1}^{n}c_{n+i}D_{n+i}=\sum_{i=1}^{n}c_{i}D_{i}+\sum_{i=1}^{n}\sum_{j=1}^{n}c_{j}\frac{a_{j,n+i}}{a_{n+i,i}}D_{n+i}.$

$D^{\prime}\leq\sum_{i=1}^{n}\sum_{j=1}^{n}c_{i}a_{ij}G_{ij}+\sum_{i=1}^{n}\sum_{j=1}^{n}c_{i}a_{i,n+j}\left(G_{i,n+j}+G_{n+j,j}\right)+\sum_{i=1}^{n}\sum_{j=1}^{n}c_{j}\frac{a_{j,n+i}}{a_{n+i,i}}a_{2n+i,n+i}G_{2n+i,n+i}.$

$D^{\prime}\leq\sum_{i=1}^{n}\sum_{j=1}^{n}c_{i}a_{ij}G_{ij}+\sum_{i=1}^{n}\sum_{j=1}^{n}c_{i}a_{i,n+j}\left(G_{i,n+j}+G_{n+j,j}\right)+\sum_{i=1}^{n}\sum_{j=1}^{n}c_{i}\frac{a_{i,n+j}}{a_{n+j,j}}a_{2n+j,n+j}G_{2n+j,n+j}.$

Let
$a^{\prime}_{ij}=\max\left(a_{ij},\frac{a_{i,n+j}}{a_{n+j,j}}(a_{2n+j,n+j})\right)$
and we obtain

$D^{\prime}\leq\sum_{i=1}^{n}\sum_{j=1}^{n}c_{i}a^{\prime}_{ij}\left(G_{ij}+G_{2n+j,n+j}\right)+\sum_{i=1}^{n}\sum_{j=1}^{n}c_{i}a_{i,n+j}\left(G_{i,n+j}+G_{n+j,j}\right).$

Let $\tilde{a}_{ij}=\max\left(a^{\prime}_{ij},a_{i,n+j}\right).$ Therefore, we
have

$D^{\prime}\leq\sum_{i=1}^{n}\sum_{j=1}^{n}c_{i}\tilde{a}_{ij}\left(G_{ij}+G_{2n+j,n+j}+G_{i,n+j}+G_{n+j,j}\right).$

Since,
$G_{ij}+G_{2n+j,n+j}+G_{i,n+j}+G_{n+j,j}=2\left(\frac{E_{j}}{\bar{E_{j}}}-\frac{E_{i}}{\bar{E_{i}}}+\log\frac{E_{i}}{\bar{E_{i}}}-\log\frac{E_{j}}{\bar{E_{j}}}\right),$
it follows that

$D^{\prime}\leq
2\sum_{i=1}^{n}\sum_{j=1}^{n}c_{i}\tilde{a}_{ij}\left(\frac{E_{j}}{\bar{E_{j}}}-\frac{E_{i}}{\bar{E_{i}}}+\log\frac{E_{i}}{\bar{E_{i}}}-\log\frac{E_{j}}{\bar{E_{j}}}\right)=0.$
(15)

Here, the equality of equation (15) follows from Theorem 3.3 stated in [35].
We can easily verify that $\\{E^{*}\\}$ is the largest invariant set for which
$D^{\prime}=0$. Using this Lyapunov function and LaSalle’s invariance
principle [36], it follows that $E^{*}$ is globally asymptotically stable in
$\Omega.$

∎

## Appendix B

## Optimal Control existence and characterization

Firstly, we are concerned about the existence of the solutions for the
controlled system (6). Further, we will prove the existence of the optimal
control [41] and then will find the optimal controls.

### Existence of an optimal control

######  Theorem 6.2.

Subject to the control system (6) with initial conditions, there exists
optimal controls $\omega^{*}_{i},\,\theta^{*}_{i},$, and $\vartheta^{*}_{i}$
such that

$J(\omega^{*}_{i},\theta^{*}_{i},\vartheta^{*}_{i})=\min_{(\omega_{i},\theta_{i},\vartheta_{i})\in
C_{ad}}J(\omega_{i},\theta_{i},\vartheta_{i}),$

if the following conditions are satisfied:

  1. (1)

The set of controls and corresponding state variables are nonempty.

  2. (2)

The control set $C_{ad}$ is compact (convex and closed).

  3. (3)

The right-hand side of the state system (6) is bounded by a linear function in
the state and control variables.

  4. (4)

The integrand
$L\left(S_{i},E_{i},I^{a}_{i},I^{s}_{i},R_{i}\right)=E_{i}(t)+I^{a}_{i}(t)+I^{s}(t)+\frac{A^{1}_{i}}{2}\left(\omega_{i}(t)\right)^{2}+\frac{A^{2}_{i}}{2}\left(\theta_{i}(t)\right)^{2}+\frac{A^{3}_{i}}{2}\left(\vartheta_{i}(t)\right)^{2}$
of the objective functional is convex on $C_{ad}$ and there exist constants
$c_{1}$ and $c_{2}$ such that

$L\left(S_{i},E_{i},I^{a}_{i},I^{s}_{i},R_{i}\right)\geq-3c_{1}+\frac{c_{2}}{2}\left(\left\lvert\omega_{i}\right\rvert^{2}+\left\lvert\theta_{i}\right\rvert^{2}+\left\lvert\vartheta_{i}\right\rvert^{2}\right)$

###### Proof.

By results discussed in [41], the existence of the optimal control can easily
be obtained and its characteristic can be discussed. First, we shall prove
that the set of controls and corresponding state variables are nonempty by
using the simplified form of the existence result in [39].
Let $Y_{1}=S_{i},Y_{2}=E_{i},Y_{3}=I^{a}_{i},Y_{4}=I^{s}_{i}$ and
$Y_{5}=R_{i}$ then from system (6), we have that
$Y^{\prime}_{k}=F_{Y_{k}}(t,Y_{1},Y_{2},Y_{3},Y_{4},Y_{5})$ with $k=1,2,3,4,5$
where $F_{Y_{k}}(t,Y_{1},Y_{2},Y_{3},Y_{4},Y_{5})$ can be obtained from the
right hand side of equations of system (6). Further, let $\omega,\theta$ and
$\vartheta$ are constant control parameters. Since all the model parameters
are constants and $Y_{1},Y_{2},Y_{3},Y_{4}$ and $Y_{5}$ are continuous.
Therefore $F_{S_{i}},F_{E_{i}},F_{I^{a}_{i}},F_{I^{s}_{i}}$ and $F_{R_{i}}$
are also continuous. Moreover, the partial derivatives $\frac{\partial
F_{Y_{k}}}{\partial Y_{k}}$ are also continuous for all $k=1,2,3,4,5.$
Therefore, a unique solution $(S_{i},E_{i},I^{a}_{i},I^{s}_{i},R_{i})$ is
satisfying the given initial conditions for the variables. Hence, there exists
a set of controls and the corresponding state variables. Thus, the first
condition is satisfied.

We have that all right-hand sides terms of system (6) are continuous and also
bounded above by a sum of bounded controls and state variables. It also can be
written as a linear function of $\omega,\theta$ and $\vartheta$ having the
time depend coefficients and state variables. From the system (6), we obtain

$\frac{dN_{i}(t)}{dt}\leq\Gamma_{i}-\mu_{i}N_{i}\implies\lim\sup_{t\rightarrow\infty}N_{i}(t)\leq\frac{\Gamma_{i}}{\mu_{i}}.$

Since, all solutions for the controlled system (6) are bounded. Therefore,
there exist positive constants $Z_{1},Z_{2},Z_{3},Z_{4},Z_{5}$ such that for
all $t\in[0,T]:$ $S_{i}(t)\leq Z_{1},E_{i}(t)\leq Z_{2},I^{a}_{i}(t)\leq
Z_{3},I^{s}_{i}(t)\leq Z_{4},$ and $R_{i}(t)\leq Z_{5}.$
Further, we consider

$\displaystyle F_{S_{i}}\leq$
$\displaystyle\Gamma_{i}+\omega_{i}(t)S_{i}(t)\sum_{j=1}^{n}\lambda_{ij}\left(\frac{I^{a}_{j}(t)+\eta_{j}I^{s}_{j}(t)}{N_{j}(t)}\right),$
$\displaystyle F_{E_{i}}\leq$ $\displaystyle\beta
S_{i}(t)(1+\eta_{j})\sum_{j=1}^{n}\lambda_{ij}-\beta
S_{i}(t)\omega_{i}(t)\sum_{j=1}^{n}\lambda_{ij}\left(\frac{I^{a}_{j}(t)+\eta_{j}I^{s}_{j}(t)}{N_{j}(t)}\right)-\theta_{i}(t)E_{i}(t),$
$\displaystyle F_{I^{a}_{i}}\leq$
$\displaystyle(1-\sigma)\alpha_{i}E_{i}(t)-\vartheta(t)I^{a}_{i}(t),$
$\displaystyle F_{I^{s}_{i}}\leq$
$\displaystyle\sigma\alpha_{i}E_{i}(t)-\vartheta(t)I^{s}_{i}(t),$
$\displaystyle F_{R_{i}}\leq$
$\displaystyle\gamma^{a}_{i}I^{a}_{i}(t)+\gamma^{s}_{i}I^{s}_{i}(t).$

Therefore, we can rewrite the control system (6) in the compact form as:

$F(S_{i},E_{i},I^{a}_{i},I^{s}_{i},R_{i})\leq\Lambda_{i}+A^{1}_{i}X^{1}_{i}(t)-B^{1}_{i}U_{i}(t),$

where
$F(S_{i},E_{i},I^{a}_{i},I^{s}_{i},R_{i})=\left(\begin{array}[]{c}F_{S_{i}}\\\
F_{E_{i}}\\\ F_{I^{a}_{i}}\\\ F_{I^{s}_{i}}\\\
F_{R_{i}}\end{array}\right),\Lambda_{i}=\left(\begin{array}[]{c}\Gamma_{i}\\\
0\\\ 0\\\ 0\\\ 0\end{array}\right),X_{i}(t)=\left(\begin{array}[]{c}S_{i}\\\
E_{i}\\\ I^{a}_{i}\\\ I^{s}_{i}\\\ R_{i}\end{array}\right),$
$U_{i}(t)=\left(\begin{array}[]{c}\omega_{i}\\\ \theta_{i}\\\
\vartheta_{i}\end{array}\right)$
$A^{1}_{i}=\left(\begin{array}[]{ccccc}0&0&0&0&0\\\
\beta(1+\eta)\sum_{j=1}^{n}\lambda_{ij}&0&0&0&0\\\
0&(1-\sigma)\alpha_{i}&0&0&0\\\ 0&\sigma\alpha_{i}&0&0&0\\\
0&0&\gamma^{a}_{i}&\gamma^{s}_{i}&0\end{array}\right)$ and
$B_{i}=\left(\begin{array}[]{ccc}-S_{i}(t)\sum_{j=1}^{n}\lambda_{ij}\left(\frac{I^{a}_{j}(t)+\eta
I^{s}_{j}(t)}{N_{j}(t)}\right)&0&0\\\
S_{i}(t)\sum_{j=1}^{n}\lambda_{ij}\left(\frac{I^{a}_{j}(t)+\eta
I^{s}_{j}(t)}{N_{j}(t)}\right)&E_{i}(t)&0\\\ 0&0&I^{a}_{i}(t)\\\
0&0&I^{s}_{i}(t)\\\ 0&0&0\\\ \end{array}\right).$ It gives a linear function
of control and state variable vectors. Since
$\left\lvert\left\lvert\Lambda_{i}\right\rvert\right\rvert\leq\Gamma_{i},$
therefore, we can obtain

$\displaystyle\left\lvert\left\lvert
F(S_{i},E_{i},I^{a}_{i},I^{s}_{i},R_{i})\right\rvert\right\rvert\leq$
$\displaystyle\left\lvert\left\lvert\Lambda_{i}\right\rvert\right\rvert+\left\lvert\left\lvert
A_{i}\right\rvert\right\rvert\left\lvert\left\lvert
X_{i}(t)\right\rvert\right\rvert+\left\lvert\left\lvert
B_{i}\right\rvert\right\rvert\left\lvert\left\lvert
U_{i}(t)\right\rvert\right\rvert$ (16)
$\displaystyle\Gamma_{i}+M^{\prime}_{i}\left(\left\lvert\left\lvert
X_{i}(t)\right\rvert\right\rvert+\left\lvert\left\lvert
U_{i}(t)\right\rvert\right\rvert\right)$

where $M^{\prime}_{i}=\max\left(\left\lvert\left\lvert
A_{i}\right\rvert\right\rvert,\left\lvert\left\lvert
B_{i}\right\rvert\right\rvert\right).$ Thus, we have that the right-hand side
is bounded by a sum of state and control vectors. Hence, third condition is
satisfied.
It is obvious that the objective function (7) is convex in $C_{ad}.$ Further,
we have that the state variables are bounded. So, let
$c_{1}=\sup_{t\in[0,T]}\left(E_{i}(t),I^{a}_{i}(t),I^{s}_{i}(t)\right)$ and
let there exist a constant $c_{2}$ such that
$c_{2}=\inf\left(A^{1}_{i},A^{2}_{i},A^{3}_{i}\right).$ Therefore, we obtain
$L\left(S_{i},E_{i},I^{a}_{i},I^{s}_{i},R_{i}\right)\geq-3c_{1}+\frac{c_{2}}{2}\left(\left\lvert\omega_{i}\right\rvert^{2}+\left\lvert\theta_{i}\right\rvert^{2}+\left\lvert\vartheta_{i}\right\rvert^{2}\right).$
Thus, the forth condition also hold. Hence, from results stated in [41]
(Theorem 2.1 of Chapter 3), we conclude that an optimal control exit for the
controlled system (6). ∎

Proof for Theorem 4.1

###### Proof.

From equation (8), we have the Hamiltonian at time $t$:

$H_{i}(t)=E_{i}(t)+I^{a}_{i}(t)+I^{s}(t)+\frac{A^{1}_{i}}{2}\left(\omega_{i}(t)\right)^{2}+\frac{A^{2}_{i}}{2}\left(\theta_{i}(t)\right)^{2}+\frac{A^{3}_{i}}{2}\left(\vartheta_{i}(t)\right)^{2}+\sum_{k=1}^{5}\psi^{k}_{i}(t)f_{k}\left(S_{i},E_{i},I^{a}_{i},I^{s}_{i},R_{i}\right)$

where

$\displaystyle f_{1}\left(S_{i},E_{i},I^{a}_{i},I^{s}_{i},R_{i}\right)=$
$\displaystyle\Gamma_{i}-\beta(1-\omega_{i}(t))S_{i}(t)\sum_{j=1}^{n}\lambda_{ij}\left(\frac{I^{a}_{j}(t)+\eta_{j}I^{s}_{j}(t)}{N_{j}}\right)-\mu_{i}S_{i},$
(17) $\displaystyle f_{2}\left(S_{i},E_{i},I^{a}_{i},I^{s}_{i},R_{i}\right)=$
$\displaystyle\beta
S_{i}(t)(1-\omega_{i}(t))\sum_{j=1}^{n}\lambda_{ij}\left(\frac{I^{a}_{j}(t)+\eta_{j}I^{s}_{j}(t)}{N_{j}}\right)-(\alpha_{i}+\mu_{i})E_{i}(t)-\theta_{i}(t)E_{i}(t),$
$\displaystyle f_{3}\left(S_{i},E_{i},I^{a}_{i},I^{s}_{i},R_{i}\right)=$
$\displaystyle(1-\sigma)\alpha_{i}E_{i}(t)-\left(\gamma^{a}_{i}+\delta^{a}_{i}+\mu_{i}\right)I^{a}_{i}(t)-\vartheta(t)I^{a}_{i}(t),$
$\displaystyle f_{4}\left(S_{i},E_{i},I^{a}_{i},I^{s}_{i},R_{i}\right)=$
$\displaystyle\sigma\alpha_{i}E_{i}(t)-\left(\gamma^{s}_{i}+\delta^{s}_{i}+\mu_{i}\right)I^{s}_{i}(t)-\vartheta(t)I^{s}_{i}(t),$
$\displaystyle f_{5}\left(S_{i},E_{i},I^{a}_{i},I^{s}_{i},R_{i}\right)=$
$\displaystyle\gamma^{a}_{i}I^{a}_{i}(t)+\gamma^{s}_{i}I^{s}_{i}(t)-\mu_{i}R_{i}(t),$

for time $t\in[0,T],$ adjoint equations and transversality conditions can be
obtained by using Pontryagin’s maximum principle, such that

$\displaystyle(\psi^{1}_{i})^{\prime}=-\frac{dH_{i}}{dS_{i}}=$
$\displaystyle\psi^{1}_{i}\mu_{i}+(\psi^{1}_{i}-\psi^{2}_{i})\beta\left(1-\omega_{i}(t)\right)\sum_{j=1}^{n}\lambda_{ij}\left(\frac{I^{a}_{j}(t)+\eta_{j}I^{s}_{j}(t)}{N_{j}}\right),$
(18) $\displaystyle(\psi^{2}_{i})^{\prime}=-\frac{dH_{i}}{dE_{i}}=$
$\displaystyle\left(\alpha_{i}+\mu_{i}+\theta_{i}(t)\right)\psi^{2}_{i}-(1-\sigma)\alpha_{i}\psi^{3}_{i}-\sigma\alpha_{i}\psi^{4}_{i},$
$\displaystyle(\psi^{3}_{i})^{\prime}=-\frac{dH_{i}}{dI^{a}_{i}}=$
$\displaystyle(\psi^{1}_{i}-\psi^{2}_{i})\beta(1-\omega_{i}(t))S_{i}(t)\sum_{j=1}^{n}\left(\frac{\lambda_{ij}}{N_{j}}\right)+\left(\gamma^{a}_{i}+\delta^{a}_{i}+\mu_{i}+\vartheta(t)\right)\psi^{3}_{i}+\gamma^{a}_{i}\psi^{5}_{i},$
$\displaystyle(\psi^{4}_{i})^{\prime}=-\frac{dH_{i}}{dI^{s}_{i}}=$
$\displaystyle(\psi^{1}_{i}-\psi^{2}_{i})\beta(1-\omega_{i}(t))S_{i}(t)\eta_{j}\sum_{j=1}^{n}\left(\frac{\lambda_{ij}}{N_{j}}\right)+\left(\gamma^{s}_{i}+\delta^{s}_{i}+\mu_{i}+\vartheta(t)\right)\psi^{4}_{i}+\gamma^{s}_{i}\psi^{5}_{i},$
$\displaystyle(\psi^{5}_{i})^{\prime}=-\frac{dH_{i}}{dR_{i}}=$
$\displaystyle\mu_{i}\psi^{5}_{i},$

with the tranversality conditions at time $T$
$\psi^{1}_{i}(T)=0,\psi^{2}_{i}(T)=1,\psi^{3}_{i}(T)=1,\psi^{4}_{i}(T)=1,\psi^{5}_{i}(T)=0.$
For $t\in[0,T]$, the optimal controls $\omega^{*}_{i},\theta^{*}_{i}$ and
$\vartheta^{*}_{i}$ can be obtained by solving the following optimality
conditions

$\frac{dH_{i}}{\omega_{i}}=0,\frac{dH_{i}}{\theta_{i}}=0,\,\text{and}\,\frac{dH_{i}}{\vartheta_{i}}=0.$

Therefore, we obtain

$\displaystyle
A^{1}_{i}\omega_{i}(t)+\left(\psi^{1}_{i}-\psi^{2}_{i}\right)\beta
S_{i}(t)\sum_{j=1}^{n}\lambda_{ij}\left(\frac{I^{a}_{j}(t)+\eta_{j}I^{s}_{j}(t)}{N_{j}}\right)$
$\displaystyle=0,$ (19) $\displaystyle
A^{2}_{i}\theta_{i}(t)-\psi^{2}_{i}E_{i}(t)$ $\displaystyle=0,$ $\displaystyle
A^{3}_{i}\vartheta_{i}(t)-\psi^{3}_{i}I^{a}_{i}(t)-\psi^{4}_{i}I^{s}_{i}(t)$
$\displaystyle=0.$

Form system of equations (19), we obtain

$\omega_{i}(t)=\frac{\left(\psi^{2}_{i}-\psi^{1}_{i}\right)\beta
S_{i}(t)}{A^{1}_{i}}\sum_{j=1}^{n}\lambda_{ij}\left(\frac{I^{a}_{j}(t)+\eta_{j}I^{s}_{j}(t)}{N_{j}}\right),$
$\theta_{i}(t)=\frac{\psi^{2}_{i}E_{i}(t)}{A^{2}_{i}},\,\vartheta_{i}(t)=\frac{\psi^{3}_{i}I^{a}_{i}(t)+\psi^{4}_{i}I^{s}_{i}(t)}{A^{3}_{i}}.$

Using the bounds of the controls which is given by $C_{ab},$ we obtain the
following optimal controls

$\omega^{*}_{i}=\max\left(0,\min\left(1,\frac{\left(\psi^{2}_{i}-\psi^{1}_{i}\right)\beta
S_{i}(t)}{A^{1}_{i}}\sum_{j=1}^{n}\lambda_{ij}\left(\frac{I^{a}_{j}(t)+\eta_{j}I^{s}_{j}(t)}{N_{j}}\right)\right)\right),$
$\theta^{*}_{i}=\max\left(0,\min\left(1,\frac{\psi^{2}_{i}E_{i}(t)}{A^{2}_{i}}\right)\right),\,\,\vartheta^{*}_{i}=\max\left(0,\min\left(1,\frac{\psi^{3}_{i}I^{a}_{i}(t)+\psi^{4}_{i}I^{s}_{i}(t)}{A^{3}_{i}}\right)\right)$

∎

## Appendix C

Figure 9: Age distribution population pyramids for India, Italy, Brazil and
USA.

Figure 10: The social contact matrices of various locations (home, school,
workplaces, and other locations) where the first column represents the
matrices for India, the second column represents the matrices for Italy, the
third column shows the matrices for Brazil and the last column for USA.