Age Dependent Epidemic Modeling of COVID-19 Outbreak in Kuwait, France, and Cameroon

Abstract

Revisiting the classical model by Ross and Kermack-McKendrick, the Susceptible–Infectious–Recovered (SIR) model used to formalize the COVID-19 epidemic, requires improvements which will be the subject of this article. The heterogeneity in the age of the populations concerned leads to considering models in age groups with specific susceptibilities, which makes the prediction problem more difficult. Basically, there are three age groups of interest which are, respectively, 0–19 years, 20–64 years, and >64 years, but in this article, we only consider two (20–64 years and >64 years) age groups because the group 0–19 years is widely seen as being less infected by the virus since this age group had a low infection rate throughout the pandemic era of this study, especially the countries under consideration. In this article, we proposed a new mathematical age-dependent (Susceptible–Infectious–Goneanewsusceptible–Recovered (SIGR)) model for the COVID-19 outbreak and performed some mathematical analyses by showing the positivity, boundedness, stability, existence, and uniqueness of the solution. We performed numerical simulations of the model with parameters from Kuwait, France, and Cameroon. We discuss the role of these different parameters used in the model; namely, vaccination on the epidemic dynamics. We open a new perspective of improving an age-dependent model and its application to observed data and parameters.

1. Introduction

1.1. Background

The influence of the median age of a country on the incidence of COVID-19 has been highlighted on numerous occasions, in particular through the appearance of many asymptomatic cases in the age groups below 50. For example, on the site of Johns Hopkins University dedicated to COVID-19 [1], data clearly shows this influence on the case fatality rate corresponding to the cumulative deaths recorded 5 months after the beginning of the outbreak (12 May 2020) vs. the median age of many countries in 2017 (Figure 1). This first observation has been confirmed by many studies in different countries [2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22] such as in France (Figure 2), and we will confirm in this paper that the age of the patients suffering from COVID-19 is a good predictor of severity.

This severity among different age groups increased COVID-19 incidence at the start of the pandemic in terms of hospitalization cases, intensive care unit (ICU) cases, prompting many countries to quickly implement intervention measures such as increasing the number of beds and ventilators in their various hospitals and setting up new isolation centers in order to mitigate its spread and save lives, particularly the elderly, who were disproportionately among those who were highly susceptible.

1.2. Related Work

The heterogeneity in the age of the populations studied here (Kuwait, France, and Cameroon) leads to considering a model in age groups with specific susceptibilities for each age class, which makes the prediction problem of the new infectious growth more difficult. Basically, there are three age groups of interest in the COVID-19 outbreak, which are, respectively, 0–19 years, 20–64 years, and >64 years, but here we only consider two (20–64 years and >64 years) age groups because the group 0–19 years is widely seen as being less infectious by SARS-CoV-2 [21,22] since this age group had a low infection rate throughout the period considered in this study, especially the countries data used for simulation.

The heterogeneity in individual age and the reaction to biological and environmental changes that have been observed in COVID-19 dynamics in terms of different reactions to vaccination by age group, severity of infection per age group, hospitalization, and intensive care unit (ICU) records show different patterns, which is why it is important to improve mathematical models for COVID-19 pandemic prediction to account for different proportions of ages in the population, which is a major factor in epidemic history [25]. In some age classes (namely young and adult classes), all infectious are neither symptomatic nor reported and some models incorporating asymptomatic [26], unreported [27], and super-spreaders [28] have been proposed for the COVID-19 outbreak. Many authors have also included age in their models as a gravity factor [29,30,31,32,33], but here, the originality is to propose an improvement of the classical Ross and Kermack-McKendrick model and use it for giving a theoretical and numerical framework for interpreting the relationship between demographic parameters such as age pyramid, fertility and mortality, and epidemiological parameters such as the basic reproduction number R₀ and vaccination rate.

In subsequent sections, we present (i) in Section 2, the materials and methods, (ii) in Section 3, the formulation of a new age-dependent dynamical model for the spread of COVID-19, (ii) in Section 4, we provide some mathematical analyses of the model, (iii) in Section 5, we carry out numerical simulations for some COVID-19 parameters from Kuwait, France, and Cameroon; lastly (iv), in Section 6 and Section 7, we provide a discussion and some conclusions.

2. Material and Methods

2.1. Material

The data material is from public epidemiologic and demographic databases [23,24] and the reference methods are both from classical demographic (such as Leslie) and epidemiologic (such as Ross and Kermack-McKendrick) models. In Figure 3 and Figure 4, we present visualizations for pandemic dynamics in different age groups and sexes for Kuwait and Cameroon, respectively, to support the motivation for this article that pandemic evolution and severity are related to age classes, and thus modeling is important as various researchers and health experts are investigating SAR-CoV-2 mutations. The countries under consideration have a higher proportion of young people, while the elderly have a lower proportion. Only 2% of the total population in Kuwait is over 65 and vulnerable to the pandemic, whereas 20.8% of the population in France is over 65 and vulnerable to the pandemic, and 2.7% of the population in Cameroon is over 65 and vulnerable to the pandemic.

2.2. Leslie Model

The first population growth model was proposed by Leslie in 1945 using the age pyramid vector $\mathrm{p}\left(\mathrm{t}\right)={\left({\mathrm{p}}_{\mathrm{j}}\left(\mathrm{t}\right)\right)}_{\mathrm{j}=1,2,\dots ,\mathrm{m}},$ where ${\mathrm{p}}_{\mathrm{j}}\left(\mathrm{t}\right)$ represents the size of the age class j at time t, with j ranging from the birth age 1 to the maximal death age m, whose discrete dynamics is governed by the matrix equation given as:

$$\begin{gathered} \mathrm{p}\left(\mathrm{t}\right)=\mathrm{Lp}\left(\mathrm{t}-1\right) \\ \mathrm{with} \mathrm{L}=\left({\mathrm{l}}_{\mathrm{jk}}\right)=\left[\begin{array}{ccccc}{\mathrm{f}}_1& {\mathrm{f}}_2& {\mathrm{f}}_3& \dots & \dots \\ {\mathrm{b}}_1& 0& 0& \dots & \dots \\ 0& {\mathrm{b}}_2& 0& \dots & \dots \\ ⋮& ⋮& \ddots & \dots & \dots \\ ⋮& ⋮& ⋮& \ddots & \dots \\ 0& 0& 0& \dots & {\mathrm{b}}_{\mathrm{m}-1}\end{array}\begin{array}{c}{\mathrm{f}}_{\mathrm{m}}\\ 0\\ 0 \\ ⋮\\ ⋮\\ 0\end{array}\right] \end{gathered}$$

where ∀ j = 1, …, m, b_j = 1 − μ_j ≤ 1 (where µ_j is the death rate at age j) is the survival probability between ages j and j + 1 and ${\mathrm{f}}_1$ is the birth rate at age j (i.e., the mean number of offsprings from an individual at age j). The dynamic stability of the Leslie system for a distance d is quantified by the tendency to return to its stationary pyramid p* after a perturbation p* + q, such as d(p*,L^m(p* + q)) < Ke^−mD tends exponentially to 0, when m goes to infinity, the parameter D being the stability module. For the $L^2$ distance to the stationary pyramid p*, D = |λ-λ’|, absolute value of the difference between the dominant and sub-dominant eigenvalues of L, i.e., λ and λ’ (λ = e^r, where r is the Malthusian growth rate, and p* is the left eigenvector of L corresponding to λ). For the distance (known as symmetrized divergence) of Kullback–Leibler, D = kH, where H is the entropy of p* and k is a constant [35,36,37,38,39,40,41,42,43,44,45,46].

2.3. Volterra Integral

The Volterra integral equation of the second kind is given as:

$$\mathrm{S}\left(\mathrm{t}\right)= \mathrm{S}\left(0\right)+{{\displaystyle \int }}_0^{\mathrm{t}}\mathcal{H}\left(\mathrm{t},\mathsf{ϵ}\right)\mathrm{S}\left(\mathsf{ϵ}\right)\mathrm{d}\mathsf{ϵ}$$

where $\mathcal{H}\left(\mathrm{t},\mathsf{ϵ}\right)$ is a kernel and $\mathrm{S}\left(\mathrm{t}\right)$ is the function to be solved.

2.4. The Ross-McKendrick SIR Model

Bernoulli [47] proposed a model for explaining the smallpox dynamics and since then, many discussions have occurred regarding the efficacy of firstly the inoculation and secondly the vaccination [48,49,50,51]. In Bernoulli’s model, the population is divided into susceptibles and immunized; the probability of belonging to them is denoted as u and v. The probability for a newborn individual susceptible and alive at age j and the probability of being immunized and alive at age j is denoted as ν_j. Ross [52] and McKendrick [53] (assisted later by W.O. Kermack) proposed an improvement of Bernoulli’s model called the SIR model with the following equations:

$$\frac{\mathrm{dS}}{\mathrm{dt}}=\mathsf{\delta }\mathrm{S}+\mathsf{\delta }\mathrm{I}+\left(\mathsf{\delta }+\mathsf{\gamma }\right)\mathrm{R}-\mathsf{\beta }\mathrm{SI}-\mu \mathrm{S},\frac{\mathrm{dI}}{\mathrm{dt}}=\mathsf{\beta }\mathrm{SI}-\left(\mu +\mathsf{\nu }\right)\mathrm{I}, \frac{\mathrm{dR}}{\mathrm{dt}}=\mathsf{\nu }\mathrm{I}-\left(\mu +\mathsf{\gamma }\right)\mathrm{R}$$

where S is the size of the susceptible, I the size of the infectious, and R the size of the recovered with the total population N defined as N = S + I + R. $\mathsf{\beta }$ is the contagion rate, $\mathsf{\delta }$ the birth rate (supposed to be equal to the death rate), $\mathsf{\gamma }$ is the loss of immune resistance, and $\mathsf{\nu }$ is the immunization rate. The basic reproduction number ${\mathrm{R}}_0=\frac{\mathsf{\beta }\mathrm{N}}{\mathsf{\nu }+\mathsf{\delta }}$ is the mean number of individuals secondarily infectious by one primary infectious. ${\mathrm{R}}_0$ predicts, if it is greater than 1, the occurrence of an epidemic wave. By defining age classes between 1 and m and by denoting ${\mathrm{S}}_{\mathrm{j}},{ \mathrm{I}}_{\mathrm{j}},{\mathrm{R}}_{\mathrm{j}}, \mathrm{j}=1, 2,\dots ,\mathrm{m},$ each subpopulation of S, I, and R at age j, we can define at any stationary state (${\mathrm{S}}^{*},{\mathrm{I}}^{*},{\mathrm{R}}^{*}$) the probabilities for a newborn individual of being alive and either susceptible, infectious or immune at age j by the following formula:

$${\mathrm{u}}^{*}\left(\mathrm{j}\right)=\frac{{\mathrm{S}}_{\mathrm{j}}{}^{*}}{{\mathrm{P}}^{*}},{ \mathrm{v}}^{*}\left(\mathrm{j}\right)=\frac{{\mathrm{I}}_{\mathrm{j}}{}^{*}}{{\mathrm{P}}^{*}},{ \mathrm{w}}^{*}\left(\mathrm{j}\right)=\frac{{\mathrm{R}}_{\mathrm{j}}{}^{*}}{{\mathrm{P}}^{*}}, \mathrm{where} \mathrm{P}*={\mathsf{\Sigma }}_{\mathrm{i}=1,\mathrm{m}} ({\mathrm{S}}_{\mathrm{i}}+{\mathrm{I}}_{\mathrm{i}}+{\mathrm{R}}_{\mathrm{i}})$$

That makes the link between the Bernoulli and the Ross and Kermack-McKendrick models, but the weakness of the latter still resides in many insufficiencies and approximations:

• When the population size of either susceptible or infectious populations tends to be very large, the quadratic term SI has to be replaced by a Michealian saturation term $\frac{\mathrm{SI}}{\left(\mathrm{k}+\mathrm{S}\right)\left( \mathrm{k} \prime +\mathrm{I}\right)}$
• the immunized infectious or healthy carriers are neglected in the absence of real information on their quantity and their influence on the spread of COVID-19;
• the total population size is supposed to be constant, the birth rate equaling the natural mortality. The Bernoulli model [47] takes implicitly into account the birth rate, and explicitly the natural mortality. The model by d’Alembert [54,55] improved Bernoulli’s model by distinguishing the specific mortality due to the infectious disease from the natural one, being more widely applicable than the model by Bernoulli which was restricted to immunizing infections. In the d’Alembert method, the only task was to calculate the survival function after eliminating the particular cause of death due to the infectious disease, but Bernoulli’s approach provided much more insight for a mechanistic interpretation of infectious disease real data;
• variables and parameters do not depend on space, i.e., neither migration nor population displacement;
• parameters do not depend on time which means no genetic adaptation of an infectious agent or a human population, even if it is very slow compared to the fast dynamics of epidemics.

Our proposed model is an improvement of the Ross and Kermack-McKendrick model by trying to compensate for a part of these defects. We will first introduce two age classes to account for adults and the elderly in the population and then take account of vaccination before applying our model to some countries chosen as examples.

3. New SIGR Model Formulation

We propose a Susceptible–Infectious–Goneanewsusceptible–Recovered (SIGR) model as an improvement of the Ross and Kermack-McKendrick models including age class and vaccination state for COVID-19 in a given population. Neglecting differences between kids and young adults, we only retain two age classes: adult and elderly, 1 and 2. We assume that all infectious are symptomatic and we consider birth and natural death rates ß and µ, as well as specific fatality rate $\mathsf{ϵ}$ due to the disease. Age groups (i = 1,2) concern individuals susceptible, infectious ${\mathrm{I}}_{\mathrm{i}},$ gone anew susceptible G_i, and a fully recovered and resistant R_i. We denote for each age group (i = 1,2) the transmission rates ${\mathsf{\phi }}_{\mathrm{i} }$, fertility, and loss of resistance rates β_ij (supposed to be equal inside an age class, for the sake of simplicity), natural death rates ${\mathsf{\mu }}_2{}^{\mathrm{S}},{ \mathsf{\mu }}_2{}^{\mathrm{I}},{ \mathsf{\mu }}_2{}^{\mathrm{G}}{ \mathrm{and} \mathsf{\mu }}_2{}^{\mathrm{R}}$, vaccination rates group ${\mathsf{\theta }}_{\mathrm{i}}^{\mathrm{G}}$ and ${\mathsf{\theta }}_{\mathrm{i}}^{\mathrm{R}},$ survival rates from age 1 to age 2 ${\mathrm{a}}_1^{\mathrm{S}},{ \mathrm{a}}_1^{\mathrm{I}},{ \mathrm{a}}_1^{\mathrm{G}}{ \mathrm{and} \mathrm{a}}_1^{\mathrm{R}},$ specific death rate due to the disease ${\mathsf{ϵ}}_{\mathrm{i} },$ relapsed rate ${ \mathsf{\eta }}_{\mathrm{i}}$ and recovery rate ${\mathsf{\gamma }}_{\mathrm{i}}$(cf.

Table 1. List of the parameters considered in the SIGR (Susceptible–Infectious–Goneanewsusceptible–Recovered) model.

Epidemiologic Parameters
Si	Susceptible individuals of age class i among age classes, adult 1 and elderly 2.
Ii	Infectious individuals of age class i
Gi	Gone anew susceptible individuals of age class i
Ri	Fully recovered and resistant individuals of age class i
Ro	Basic reproduction number
ϕi	Transmission rate of age class i
ηi	Relapsed rate of age class i
γi	Recovering rate of age class i
a1K	Survival rate from age 1 to age 2 for compartment K of age class 1
µ2K	Natural death rate of compartment K of age class 2
Demographic Parameters
βi1	Birth and loss of resistance of recovered from compartments of age class i
β22	Loss of resistance of recovered rate from compartments of age class 2
θiK	Vaccination rate from compartment K of age class i
${\mathsf{ϵ}}_{\mathrm{i}}$	Specific fatality rate $\mathsf{ϵ}$ due to the disease of age class i

).

The description above can be illustrated by the following set of non-linear differential Equation (1), while the graphical representation of the model is given in Figure 5:

$$\begin{gathered} {\mathrm{dS}}_1/\mathrm{dt}={\mathsf{\beta }}_{11} ({\mathrm{S}}_1+{\mathrm{G}}_1+{\mathrm{R}}_1)+{\mathsf{\beta }}_{21}({\mathrm{S}}_2+{\mathrm{G}}_2+{\mathrm{R}}_2) – ({\mathrm{a}}_1^{\mathrm{S}}+{\mathsf{\theta }}_1^{\mathrm{G}}+{ \mathsf{\theta }}_1^{\mathrm{R}}+{\mathsf{\phi }}_1({\mathrm{I}}_1 +{\mathrm{I}}_2)){\mathrm{S}}_1 \\ {\mathrm{dS}}_2/\mathrm{dt}={\mathrm{a}}_1^{\mathrm{S}}{\mathrm{S}}_1 – {\mathsf{\mu }}_2{}^{\mathrm{S}}{\mathrm{S}}_2 - ({\mathsf{\theta }}_2^{\mathrm{G}}+{ \mathsf{\theta }}_2^{\mathrm{R}}+{\mathsf{\phi }}_2({\mathrm{I}}_1+{\mathrm{I}}_2))\mathrm{S}{2}_{ }+{\mathsf{\beta }}_{22}({\mathrm{G}}_2+{\mathrm{R}}_2) \\ {\mathrm{dI}}_1/\mathrm{dt}={\mathsf{\phi }}_1\left({\mathrm{I}}_1+{\mathrm{I}}_2\right){\mathrm{S}}_1-({\mathrm{a}}_1^{\mathrm{I}}+{\mathsf{\eta }}_1){\mathrm{I}}_1 \\ {\mathrm{dI}}_2/\mathrm{dt}={\mathsf{\phi }}_2\left({\mathrm{I}}_1+{\mathrm{I}}_2\right){\mathrm{S}}_2+{\mathrm{a}}_1^{\mathrm{I}}{\mathrm{I}}_1-{\mathsf{\eta }}_2{ \mathrm{I}}_2 – {\mathsf{\mu }}_2{}^{\mathrm{I}}{\mathrm{I}}_2 \\ {\mathrm{dG}}_1/\mathrm{dt}={ \mathsf{\eta }}_1{\mathrm{I}}_1+{ \mathsf{\theta }}_1^{\mathrm{G}}{\mathrm{S}}_1 – ({\mathrm{a}}_1^{\mathrm{G}}+{\mathsf{\gamma }}_1+{\mathsf{ϵ}}_1){\mathrm{G}}_1 \\ {\mathrm{dG}}_2/\mathrm{dt}={ \mathsf{\eta }}_2{\mathrm{I}}_2+{ \mathrm{a}}_1^{\mathrm{G}}{\mathrm{G}}_1+{ \mathsf{\theta }}_2^{\mathrm{G}}{\mathrm{S}}_2 – ({\mathsf{\gamma }}_2+{ \mathsf{\mu }}_2{}^{\mathrm{G}}+{\mathsf{ϵ}}_2+{\mathsf{\beta }}_{22}){\mathrm{G}}_2 \\ {\mathrm{dR}}_1/\mathrm{dt}={\mathsf{\gamma }}_1{\mathrm{G}}_1+{\mathsf{\theta }}_1^{\mathrm{R}}{\mathrm{S}}_1-{ \mathrm{a}}_1^{\mathrm{R}}{\mathrm{R}}_1 \\ {\mathrm{dR}}_2/\mathrm{dt}={\mathsf{\gamma }}_2{\mathrm{G}}_2+{\mathsf{\theta }}_2^{\mathrm{R}}{\mathrm{S}}_2+{\mathrm{a}}_1^{\mathrm{R}}{\mathrm{R}}_1 – ({\mathsf{\mu }}_2{}^{\mathrm{R}}+{\mathsf{\beta }}_{22}){\mathrm{R}}_2 \end{gathered}$$

(1)

with ${\mathrm{S}}_1\left(\mathrm{t}\right)\ge 0,{\mathrm{S}}_2\left(\mathrm{t}\right)\ge 0,{ \mathrm{I}}_1\left(\mathrm{t}\right)\ge 0,{ \mathrm{I}}_2\left(\mathrm{t}\right)\ge 0,{ \mathrm{G}}_1\left(\mathrm{t}\right)\ge 0,{ \mathrm{G}}_2\left(\mathrm{t}\right)\ge 0,{ \mathrm{R}}_1\left(\mathrm{t}\right)\ge 0,{ \mathrm{R}}_2\left(\mathrm{t}\right)\ge 0.$

We will suppose in the following that the loss of resistance rates is neglectable.

4. Mathematical Analysis of the Model

4.1. Positivity of the Solution

Lemma 1.

Let the initial condition be given as follows:

$\left\{{\mathrm{S}}_1\left(0\right),{\mathrm{S}}_2\left(0\right),{ \mathrm{I}}_1\left(0\right),{ \mathrm{I}}_2\left(0\right),{ \mathrm{G}}_1\left(0\right),{ \mathrm{G}}_2\left(0\right),{ \mathrm{R}}_1\left(0\right),{ \mathrm{R}}_2\left(0\right)\ge 0\right\},$then solutions of system of

Equation (1) are positive for all$t>0.$

Proof 

From the first equation in the model Equation (1), we obtain:

$${\mathrm{dS}}_1/\mathrm{dt}={\mathsf{\beta }}_{11} ({\mathrm{S}}_1+{\mathrm{G}}_1+{\mathrm{R}}_1)+{\mathsf{\beta }}_{21} ({\mathrm{S}}_2+{\mathrm{G}}_2+{\mathrm{R}}_2) – ({\mathrm{a}}_1^{\mathrm{S}}+{\mathsf{\theta }}_1^{\mathrm{G}}+{ \mathsf{\theta }}_1^{\mathrm{R}}+{\mathsf{\phi }}_1({\mathrm{I}}_1+{\mathrm{I}}_2)){\mathrm{S}}_1,$$

Hence, dS₁/dt $\ge ({\mathsf{\beta }}_{11 }-{\mathrm{a}}_1^{\mathrm{S}}-{\mathsf{\theta }}_1^{\mathrm{G}}-{\mathsf{\theta }}_1^{\mathrm{R}}$) S₁

By using the separating variable method and then, integrating, we obtain:

$$\int \frac{{\mathrm{dS}}_1/\mathrm{dt}}{{\mathrm{S}}_1} \ge \int \left({\mathsf{\beta }}_{11 }-{\mathrm{a}}_1^{\mathrm{S}}-{\mathsf{\theta }}_1^{\mathrm{G}}-{\mathsf{\theta }}_1^{\mathrm{R}}\right)\mathrm{dt}$$

$$\ln {\mathrm{S}}_1\ge \left({\mathsf{\beta }}_{11 }-{\mathrm{a}}_1^{\mathrm{S}}-{\mathsf{\theta }}_1^{\mathrm{G}}-{\mathsf{\theta }}_1^{\mathrm{R}}\right)\mathrm{t}+\mathrm{k}$$

Finally, by writing S₁(0) = e^k, we have: ${\mathrm{S}}_1\left(\mathrm{t}\right)\ge {\mathrm{S}}_1\left(0\right){\mathrm{e}}^{\left({\mathsf{\beta }}_{11}-{\mathrm{a}}_1^{\mathrm{S}}-{\mathsf{\theta }}_1^{\mathrm{G}}-{\mathsf{\theta }}_1^{\mathrm{R}}\right)\mathrm{t}}\ge 0$.

By applying the same process to other equations in (1), we have:

$${\mathrm{S}}_2\left(\mathrm{t}\right)\ge {\mathrm{S}}_2\left(0\right){\mathrm{e}}^{\left({\mathrm{a}}_1^{\mathrm{S}}-{\mathsf{\theta }}_2^{\mathrm{G}}-{\mathsf{\theta }}_2^{\mathrm{R}}-{\mathsf{\mu }}_2{}^{\mathrm{S}}\right)\mathrm{t}}\ge 0$$

$${\mathrm{I}}_1\left(\mathrm{t}\right)\ge {\mathrm{I}}_1\left(0\right){\mathrm{e}}^{-\left({\mathrm{a}}_1^{\mathrm{I}}+{\mathsf{\eta }}_1\right)\mathrm{t}}\ge 0$$

$${\mathrm{I}}_2\left(\mathrm{t}\right)\ge {\mathrm{I}}_2\left(0\right){\mathrm{e}}^{-\left({\mathsf{\eta }}_2{ \mathrm{I}}_2+{ \mathsf{\mu }}_2{}^{\mathrm{I}}\right)\mathrm{t}}\ge 0$$

$${\mathrm{G}}_1\left(\mathrm{t}\right)\ge {\mathrm{G}}_1\left(0\right){\mathrm{e}}^{-\left({\mathrm{a}}_1^{\mathrm{G}}+{\mathsf{\gamma }}_1+{\mathsf{ϵ}}_1\right)\mathrm{t}}\ge 0$$

$${\mathrm{G}}_2\left(\mathrm{t}\right)\ge {\mathrm{G}}_2\left(0\right){\mathrm{e}}^{-\left({\mathsf{\gamma }}_2+{ \mathsf{\mu }}_2{}^{\mathrm{G}}+{\mathsf{ϵ}}_2\right)\mathrm{t}}\ge 0$$

$${\mathrm{R}}_1\left(\mathrm{t}\right)\ge {\mathrm{R}}_1\left(0\right){\mathrm{e}}^{-{\mathrm{a}}_1^{\mathrm{R}}\mathrm{t}}\ge 0$$

$${\mathrm{R}}_2\left(\mathrm{t}\right)\ge {\mathrm{R}}_2\left(0\right){\mathrm{e}}^{-{\mathsf{\mu }}_2{}^{\mathrm{R}}\mathrm{t}}\ge 0$$

Then, the solutions of Equation (1) $\left\{{\mathrm{S}}_1\left(\mathrm{t}\right),{\mathrm{S}}_2\left(\mathrm{t}\right),{ \mathrm{I}}_1\left(\mathrm{t}\right),{ \mathrm{I}}_2\left(\mathrm{t}\right),{ \mathrm{G}}_1\left(\mathrm{t}\right),{ \mathrm{G}}_2\left(\mathrm{t}\right),{ \mathrm{R}}_1\left(\mathrm{t}\right),{ \mathrm{R}}_2\left(\mathrm{t}\right)\right\}$ are positive for all $\mathrm{t}>0.$ □

4.2. Boundedness of the Solution

Let denote by S the total size of all individuals:

S(t) = S₁(t) + S₂(t) + I₁(t) + I₂(t) + G₁(t) + G₂(t) + R₁(t) + R₂(t)

Then, by adding all the model Equations (1), we have:

$$\mathrm{dS}/\mathrm{dt}={\mathsf{\beta }}_{11} ({\mathrm{S}}_1+{\mathrm{G}}_1+{\mathrm{R}}_1)+{\mathsf{\beta }}_{21} ({\mathrm{S}}_2+{\mathrm{G}}_2+{\mathrm{R}}_2) – ({{\mu }_2}^{\mathrm{S}}{\mathrm{S}}_2+{{\mu }_2}^{\mathrm{I}}{\mathrm{I}}_2+{{\mu }_2}^{\mathrm{G}}{{\mathrm{G}}_2}^{ }+ϵ{\mathrm{G}}_2+{{\mu }_2}^{\mathrm{R}}{{\mathrm{R}}_2}^{}$$

Let us denote µ = inf{µ₂^S, µ₂^I, µ₂^G+$𝜖$, µ₂^R}. By neglecting the fecundity rate of the young class and if β = β₂₁ ≤ µ, we have:

dS/dt ≤ β(S₂ + I₂ + G₂ + R₂) – µ(S₂ + I₂ + G₂ + R₂) ≤ 0,

Then, we can conclude that the total size S is bounded, which implies the boundedness of the partial sizes ${\mathrm{S}}_1\left(\mathrm{t}\right),{\mathrm{S}}_2\left(\mathrm{t}\right),{ \mathrm{I}}_1\left(\mathrm{t}\right),{ \mathrm{I}}_2\left(\mathrm{t}\right),{ \mathrm{G}}_1\left(\mathrm{t}\right),{ \mathrm{G}}_2\left(\mathrm{t}\right),{ \mathrm{R}}_1\left(\mathrm{t}\right),{ \mathrm{R}}_2\left(\mathrm{t}\right)$.

4.3. Disease-Free (Eradication) Equilibrium

Setting the right hand side of the model equations (1) to zero, i.e.,

$${\mathrm{dS}}_1/\mathrm{dt}={\mathrm{dS}}_2/\mathrm{dt}={\mathrm{dI}}_1/\mathrm{dt}={\mathrm{dI}}_2/\mathrm{dt}={\mathrm{G}}_1/\mathrm{dt}={\mathrm{dG}}_2/\mathrm{dt}={\mathrm{dR}}_1/\mathrm{dt}={\mathrm{dR}}_2/\mathrm{dt}=0$$

and supposing that all infectious class sizes are equal to zero which means that there is no disease eradication in the studied population, thus the disease-free equilibrium state is:

$$({\mathrm{S}}_1{}^{*},{ \mathrm{S}}_2{}^{*},{\mathrm{I}}_1{}^{*},{\mathrm{I}}_2{}^{*}{\mathrm{G}}_1{}^{*},{\mathrm{G}}_2{}^{*},{\mathrm{R}}_1{}^{*},{\mathrm{R}}_2{}^{*})=\left\{\frac{{\mathsf{\beta }}_{21 }{\mathrm{S}}_2}{{\mathrm{a}}_1^{\mathrm{S}}+{\mathsf{\theta }}_1^{\mathrm{G}}+{ \mathsf{\theta }}_1^{\mathrm{R}}-{\mathsf{\beta }}_{11 }},\frac{{\mathrm{a}}_1^{\mathrm{S}}{\mathrm{S}}_1}{{\mathsf{\mu }}_2{}^{\mathrm{S}}+{\mathsf{\theta }}_2^{\mathrm{G}}+{ \mathsf{\theta }}_2^{\mathrm{R}}},0,0,0,0,0,0\right\}$$

4.4. Proof of Stability of the Endemic State

Let us fix a set of values for the model parameters, where k is a scale parameter:

$${\mathsf{\phi }}_{1 }={ \mathsf{\phi }}_{2 }=\frac{4{\mathrm{k}}^2}{100},{ \mathsf{\mu }}_2{}^{\mathrm{S}}={\mathsf{ϵ}}_{1 }={\mathsf{\beta }}_{22 }=0,{ \mathsf{\mu }}_2{}^{\mathrm{I}}=\frac{499\mathrm{k}}{100},{ \mathsf{\mu }}_2{}^{\mathrm{G}}=\frac{49\mathrm{k}}{100},{ \mathsf{\mu }}_2{}^{\mathrm{R}}=\frac{2\mathrm{k}}{5},{ \mathsf{\beta }}_{11 }={ \mathsf{\beta }}_{21 }=\frac{\mathrm{k}}{5},$$

$${\mathsf{\theta }}_1^{\mathrm{G}}={\mathsf{\theta }}_1^{\mathrm{R}}={\mathsf{\gamma }}_1={\mathsf{\gamma }}_2=0.1\mathrm{k},{ \mathsf{\theta }}_2^{\mathrm{G}}={\mathsf{\theta }}_2^{\mathrm{R}}=\mathrm{k},{ \mathrm{a}}_1^{\mathrm{S}}={ \mathrm{a}}_1^{\mathrm{I}}={ \mathrm{a}}_1^{\mathrm{G}}={ \mathrm{a}}_1^{\mathrm{R}}=\frac{98\mathrm{k}}{96}, {\mathsf{ϵ}}_{2 }=\frac{4\mathrm{k}}{5},{ \mathsf{\eta }}_1={ \mathsf{\eta }}_2=0.2\mathrm{k}$$

The two stationary points are labeled with * (resp. **) for the eradication (resp. endemic) state:

$$({\mathrm{S}}_1{}^{*},{ \mathrm{S}}_2{}^{*},{\mathrm{I}}_1{}^{*},{\mathrm{I}}_2{}^{*}{\mathrm{G}}_1{}^{*},{\mathrm{G}}_2{}^{*},{\mathrm{R}}_1{}^{*},{\mathrm{R}}_2{}^{*})=\left(\frac{{\mathsf{\beta }}_{21 }{\mathrm{S}}_2}{{\mathrm{a}}_1^{\mathrm{S}}+{\mathsf{\theta }}_1^{\mathrm{G}}+{ \mathsf{\theta }}_1^{\mathrm{R}}-{\mathsf{\beta }}_{11 }},\frac{{\mathrm{a}}_1^{\mathrm{S}}{\mathrm{S}}_1}{{\mathsf{\mu }}_2{}^{\mathrm{S}}+{\mathsf{\theta }}_2^{\mathrm{G}}+{ \mathsf{\theta }}_2^{\mathrm{R}}},0,0,0,0,0,0\right)$$

and

$$({\mathrm{S}}_1{}^{**},{ \mathrm{S}}_2{}^{**},{\mathrm{I}}_1{}^{**},{\mathrm{I}}^{**},{\mathrm{G}}_1{}^{**},{\mathrm{G}}_2{}^{**},{\mathrm{R}}_1{}^{**},{\mathrm{R}}_2{}^{**})=\left(10,10,10,10,20,40,15,15\right), \mathrm{if} \mathrm{k}=1.$$

Let us show that the endemic state is locally stable. With the chosen parameter values and k = 1, from the model Equation (1), we have:

$$\begin{gathered} {\mathrm{dS}}_1/\mathrm{dt}=0.2{\mathrm{G}}_1-1.02{\mathrm{S}}_1+0.2{\mathrm{R}}_1+0.2{\mathrm{S}}_2+0.2{\mathrm{G}}_2+0.2{\mathrm{R}}_2-0.04{\mathrm{I}}_1{\mathrm{S}}_1-0.04{\mathrm{I}}_2{\mathrm{S}}_1 \\ {\mathrm{dS}}_2/\mathrm{dt}=1.02{\mathrm{S}}_1-2{\mathrm{S}}_2-0.04{\mathrm{I}}_1{\mathrm{S}}_2-0.04{\mathrm{I}}_2{\mathrm{S}}_2 \\ {\mathrm{dI}}_1/\mathrm{dt}=0.04{\mathrm{I}}_1{\mathrm{S}}_1+0.04{\mathrm{I}}_2{\mathrm{S}}_1-1.04{\mathrm{I}}_1 \\ {\mathrm{dI}}_2/\mathrm{dt}=0.04{\mathrm{I}}_1{\mathrm{S}}_2+0.04{\mathrm{I}}_2{\mathrm{S}}_2+1.02{\mathrm{I}}_1-5.19{\mathrm{I}}_2 \end{gathered}$$

(2)

$${\mathrm{dG}}_1/\mathrm{dt}=0.2{\mathrm{I}}_1+0.1{\mathrm{S}}_1-1.12{\mathrm{G}}_1$$

$${\mathrm{dG}}_2/\mathrm{dt}=0.2{\mathrm{I}}_2+1.02{\mathrm{G}}_1+{ \mathrm{S}}_2-1.39{\mathrm{G}}_2$$

$${\mathrm{dR}}_1/\mathrm{dt}=0.1{\mathrm{G}}_1+0.1{\mathrm{S}}_1-1.02{\mathrm{R}}_1$$

$${\mathrm{dR}}_2/\mathrm{dt}=0.1{\mathrm{G}}_2+{\mathrm{S}}_2+1.02{\mathrm{R}}_1 – 0.4{\mathrm{R}}_2$$

By calculating the Jacobian matrix M of the system of Equation (2) at the endemic state, where I is the identity matrix and finding the roots of its characteristic polynomial P_M, we have for the second stationary point the expression as follows:

$$\mathrm{M}-\mathsf{\lambda }\mathrm{I}=\left(\begin{array}{cccccccc}-1.02-\mathsf{\lambda }& 0.2& -0.4& -0.4& 0.2& 0.2& 0.2& 0.2\\ 1.02& -2-\mathsf{\lambda }& -0.4& -0.4& 0& 0& 0& 0\\ 0.8& 0& -1.04-\mathsf{\lambda }& 0& 0& 0& 0& 0\\ 0& 0.8& 1.02& -5.19-\mathsf{\lambda }& 0& 0& 0& 0\\ 0.1& 0& 0.2& 0& -1.12-\mathsf{\lambda }& 0& 0& 0\\ 0& 1& 0& 0.2& 1.02& -1.39-\mathsf{\lambda }& 0& 0\\ 0.1& 0& 0& 0& 0.1& 0& -1.02-\mathsf{\lambda }& 0\\ 0& 1& 0& 0& 0& 0.1& 1.02& -0.4-\mathsf{\lambda }\end{array}\right)$$

The roots of the characteristic polynomial of M, P_M, satisfy:

$${\mathrm{P}}_{\mathrm{M}}\left(\mathsf{\lambda })=\det (\mathrm{M}-\mathsf{\lambda }\mathrm{I}\right) = (\mathsf{\lambda }+0.706923)(\mathsf{\lambda }+5.13103)({\mathsf{\lambda }}^2+0.973528\mathsf{\lambda }+0.34052)$$

$$\left({\mathsf{\lambda }}^2+2.0371\mathsf{\lambda }+1.03861\right)\left({\mathsf{\lambda }}^2+3.69142\mathsf{\lambda } +3.53083\right)=0$$

The real parts of the eigenvalues of the matrix M are all negative, equal to:

−0706923, −5.13103, −0.486764, −1.01855, −1.84571

Hence, with the set of chosen parameter values, the stability of the endemic state is proved.

4.5. Existence and Unicity of the Solution

We want to establish existence and unicity of the solution for the model Equation (1). Let us denote the second member of (1) by $\mathcal{H}=\left({\mathcal{H}}_1,{\mathcal{H}}_2,{\mathcal{H}}_3,{\mathcal{H}}_4,{\mathcal{H}}_5,{\mathcal{H}}_6,{\mathcal{H}}_7,{\mathcal{H}}_8\right)$ and the state vector by $\mathrm{Z}=({\mathrm{S}}_1,{\mathrm{S}}_2,{ \mathrm{I}}_1,{ \mathrm{I}}_2,{ \mathrm{G}}_1,{ \mathrm{G}}_2,{ \mathrm{R}}_1,{ \mathrm{R}}_2) , \mathrm{with} \mathrm{Z}\left(0\right)={ \mathrm{Z}}_{\mathrm{o}}, \mathrm{Z}\left(\mathrm{t},\mathsf{ϵ}\right)=\mathcal{H}\left(\mathrm{t},\mathsf{ϵ}\right)\mathrm{Z}\left(\mathrm{t}\right).$ The Volterra integral equation formulation of the model Equation (1) is the following:

$$\begin{gathered} {\mathrm{S}}_1\left(\mathrm{t}\right)={\mathrm{S}}_1\left(0\right)+{{\displaystyle \int }}_0^{\mathrm{t}}{\mathcal{H}}_1\left(\mathrm{t},\mathsf{ϵ}\right){\mathrm{S}}_1\left(\mathsf{ϵ}\right)\mathrm{d}\mathsf{ϵ} \\ {\mathrm{S}}_2\left(\mathrm{t}\right)={\mathrm{S}}_2\left(0\right)+{{\displaystyle \int }}_0^{\mathrm{t}}{\mathcal{H}}_2\left(\mathrm{t},\mathsf{ϵ}\right){\mathrm{S}}_2\left(\mathsf{ϵ}\right)\mathrm{d}\mathsf{ϵ} \\ {\mathrm{I}}_1\left(\mathrm{t}\right)={\mathrm{I}}_1\left(0\right)+{{\displaystyle \int }}_0^{\mathrm{t}}{\mathcal{H}}_3\left(\mathrm{t},\mathsf{ϵ}\right){\mathrm{I}}_1\left(\mathsf{ϵ}\right)\mathrm{d}\mathsf{ϵ} \\ {\mathrm{I}}_2\left(\mathrm{t}\right)={\mathrm{I}}_2\left(0\right)+{{\displaystyle \int }}_0^{\mathrm{t}}{\mathcal{H}}_4\left(\mathrm{t},\mathsf{ϵ}\right){\mathrm{I}}_2\left(\mathsf{ϵ}\right)\mathrm{d}\mathsf{ϵ} \\ { \mathrm{G}}_1\left(\mathrm{t}\right)={\mathrm{G}}_1\left(0\right)+{{\displaystyle \int }}_0^{\mathrm{t}}{\mathcal{H}}_5\left(\mathrm{t},\mathsf{ϵ}\right){\mathrm{G}}_1\left(\mathsf{ϵ}\right)\mathrm{d}\mathsf{ϵ} \end{gathered}$$

(3)

$$\begin{gathered} { \mathrm{G}}_2\left(\mathrm{t}\right)={\mathrm{G}}_2\left(0\right)+{{\displaystyle \int }}_0^{\mathrm{t}}{\mathcal{H}}_6\left(\mathrm{t},\mathsf{ϵ}\right){\mathrm{G}}_2\left(\mathsf{ϵ}\right)\mathrm{d}\mathsf{ϵ} \\ {\mathrm{R}}_1\left(\mathrm{t}\right)={\mathrm{R}}_1\left(0\right)+{{\displaystyle \int }}_0^{\mathrm{t}}{\mathcal{H}}_7\left(\mathrm{t},\mathsf{ϵ}\right){\mathrm{R}}_1\left(\mathsf{ϵ}\right)\mathrm{d}\mathsf{ϵ} \\ {\mathrm{R}}_2\left(\mathrm{t}\right)={\mathrm{R}}_2\left(0\right)+{{\displaystyle \int }}_0^{\mathrm{t}}{\mathcal{H}}_8\left(\mathrm{t},\mathsf{ϵ}\right){\mathrm{R}}_2\left(\mathsf{ϵ}\right)\mathrm{d}\mathsf{ϵ} \end{gathered}$$

Let ${\mathcal{S}}_1$ be the solution for ${\mathrm{S}}_1 \mathrm{and}$ kernels ${\mathcal{H}}_{\mathrm{I}}, \mathrm{i}=1,2,\dots , 8$ satisfy Lipschitz conditions: ${\sup }_{0<\mathrm{t}\le 1}\Vert {\mathrm{S}}_1\Vert \le {\mathrm{c}}_1$, ${ \sup }_{0<\mathrm{t}\le 1}\Vert {\mathrm{S}}_2\Vert \le {\mathrm{c}}_2$, ${ \sup }_{0<\mathrm{t}\le 1}\Vert {\mathrm{I}}_1\Vert {\mathrm{c}}_3,{ \sup }_{0<\mathrm{t}\le 1}\Vert {\mathrm{I}}_2\Vert {\mathrm{c}}_4,{ \sup }_{0<\mathrm{t}\le 1}\Vert {\mathrm{G}}_1\Vert \le {\mathrm{c}}_5,$ ${\sup }_{0<\mathrm{t}\le 1}\Vert {\mathrm{G}}_2\Vert \le {\mathrm{c}}_6,{ \sup }_{0<\mathrm{t}\le 1}\Vert {\mathrm{R}}_1\Vert \le {\mathrm{c}}_7,{ \sup }_{0<\mathrm{t}\le 1}\Vert {\mathrm{R}}_2\Vert \le {\mathrm{c}}_{8 }{\mathrm{and} \mathrm{c}}_{\mathrm{i}}>0, \mathrm{for} \mathrm{i}=1, 2,\dots , 8.$ Then, following inequalities hold, using triangle inequality and properties of the ${\mathcal{H}}_{\mathrm{i}}^{’}\mathrm{s} \mathrm{norm}$:

$$\begin{gathered} \Vert {\mathcal{H}}_1\left({\mathrm{S}}_1\right)-{\mathcal{H}}_1\left({\mathcal{S}}_1\right)\Vert \le \left[{\mathsf{\beta }}_{11 }– \left({\mathrm{a}}_1^{\mathrm{S}}+{\mathsf{\theta }}_1^{\mathrm{G}}+{ \mathsf{\theta }}_1^{\mathrm{R}}+{ \mathsf{\phi }}_1\left({\mathrm{I}}_1+{\mathrm{I}}_2\right)\right)\right]{\mathrm{S}}_1-\left[{\mathsf{\beta }}_{11 }– ({\mathrm{a}}_1^{\mathrm{S}}+{\mathsf{\theta }}_1^{\mathrm{G}}+{ \mathsf{\theta }}_1^{\mathrm{R}}+{ \mathsf{\phi }}_1\left({\mathrm{I}}_1+{\mathrm{I}}_2\right))\right]{\mathcal{S}}_1\Vert \\ {}_{ } \\ \le \Vert {\mathsf{\beta }}_{11 }–({\mathrm{a}}_1^{\mathrm{S}}+{\mathsf{\theta }}_1^{\mathrm{G}}+{ \mathsf{\theta }}_1^{\mathrm{R}}+{ \mathsf{\phi }}_1\left({\mathrm{I}}_1+{\mathrm{I}}_2\right))\Vert \Vert {\mathrm{S}}_1-{\mathcal{S}}_1\Vert \\ \le \Vert {\mathsf{\beta }}_{11 }– ({\mathrm{a}}_1^{\mathrm{S}}+{\mathsf{\theta }}_1^{\mathrm{G}}+{ \mathsf{\theta }}_1^{\mathrm{R}}+{ \mathsf{\phi }}_1{\mathrm{c}}_3+{\mathsf{\phi }}_1{\mathrm{c}}_4)\Vert \Vert {\mathrm{S}}_1-{\mathcal{S}}_1\Vert \\ ={\partial }_1\Vert {\mathrm{S}}_1-{\mathcal{S}}_1\Vert \end{gathered}$$

where $\mathsf{\xi }={\mathsf{\beta }}_{11 }$ (G₁ + R₁) + β₂₁ (S₂ + G₂ + R₂), ${\partial }_1={\mathsf{\beta }}_{11 }–{ \mathrm{a}}_1^{\mathrm{S}}-{\mathsf{\theta }}_1^{\mathrm{G}}-{ \mathsf{\theta }}_1^{\mathrm{R}}-{ \mathsf{\phi }}_1{\mathrm{c}}_3-{\mathsf{\phi }}_1{\mathrm{c}}_4$.

Therefore, ${\mathcal{H}}_1$ satisfies the Lipschitz conditions. We can show in the same way that other functions ${\mathcal{H}}_{\mathrm{i}}, \mathrm{i}=2,\dots , 8$ in the model Equation (2) satisfy the Lipschitz conditions as follows:

$$\begin{gathered} \Vert {\mathcal{H}}_2\left({\mathrm{S}}_2\right)-{\mathcal{H}}_2\left({\mathcal{S}}_2\right)\le \Vert {\partial }_2\Vert {\mathrm{S}}_2-{\mathcal{S}}_2\Vert \\ \Vert {\mathcal{H}}_3\left({\mathrm{I}}_1\right)-{\mathcal{H}}_3\left({\mathcal{L}}_1\right)\le \Vert {\partial }_3\Vert {\mathrm{I}}_1-{\mathcal{L}}_1\Vert \\ \Vert {\mathcal{H}}_4\left({\mathrm{I}}_2\right)-{\mathcal{H}}_4\left({\mathcal{L}}_2\right)\le \Vert {\partial }_4\Vert {\mathrm{I}}_2-{\mathcal{L}}_2\Vert \\ \Vert {\mathcal{H}}_5\left({\mathrm{G}}_1\right)-{\mathcal{H}}_5\left({\mathcal{G}}_1\right)\le \Vert {\partial }_5\Vert {\mathrm{G}}_1-{\mathcal{G}}_1\Vert \\ \Vert {\mathcal{H}}_6\left({\mathrm{G}}_2\right)-{\mathcal{H}}_6\left({\mathcal{G}}_2\right)\le \Vert {\partial }_6\Vert {\mathrm{G}}_2-{\mathcal{G}}_2\Vert \\ \Vert {\mathcal{H}}_7\left({\mathrm{R}}_1\right)-{\mathcal{H}}_7\left({\Re }_1\right)\le \Vert {\partial }_7\Vert {\mathrm{R}}_1-{\Re }_1\Vert \\ \Vert {\mathcal{H}}_8\left({\mathrm{R}}_2\right)-{\mathcal{H}}_8\left({\Re }_2\right)\le \Vert {\partial }_8\Vert {\mathrm{R}}_2-{\Re }_2\Vert \end{gathered}$$

Let us now consider the following Neumann series:

$$\begin{gathered} {\mathrm{S}}_{1\mathrm{m}}\left(\mathrm{t}\right)={\mathrm{S}}_1\left(0\right)+{{\displaystyle \int }}_0^{\mathrm{t}}{\mathcal{H}}_1\left(\mathrm{t},\mathsf{ϵ}\right){\mathrm{S}}_{1\mathrm{m}-1}\left(\mathsf{ϵ}\right)\mathrm{d}\mathsf{ϵ} \\ {\mathrm{S}}_{2\mathrm{m}}\left(\mathrm{t}\right)={\mathrm{S}}_2\left(0\right)+{{\displaystyle \int }}_0^{\mathrm{t}}{\mathcal{H}}_2\left(\mathrm{t},\mathsf{ϵ}\right){\mathrm{S}}_{2\mathrm{m}-1}\left(\mathsf{ϵ}\right)\mathrm{d}\mathsf{ϵ} \\ {\mathrm{I}}_{1\mathrm{m}}\left(\mathrm{t}\right)={\mathrm{I}}_1\left(0\right)+{{\displaystyle \int }}_0^{\mathrm{t}}{\mathcal{H}}_3\left(\mathrm{t},\mathsf{ϵ}\right){\mathrm{I}}_{1\mathrm{m}-1}\left(\mathsf{ϵ}\right)\mathrm{d}\mathsf{ϵ} \\ {\mathrm{I}}_{2\mathrm{m}}\left(\mathrm{t}\right)={\mathrm{I}}_2\left(0\right)+{{\displaystyle \int }}_0^{\mathrm{t}}{\mathcal{H}}_4\left(\mathrm{t},\mathsf{ϵ}\right){\mathrm{I}}_{2\mathrm{m}-1}\left(\mathsf{ϵ}\right)\mathrm{d}\mathsf{ϵ} \\ {\mathrm{G}}_{1\mathrm{m}}\left(\mathrm{t}\right)={\mathrm{G}}_1\left(0\right)+{{\displaystyle \int }}_0^{\mathrm{t}}{\mathcal{H}}_5\left(\mathrm{t},\mathsf{ϵ}\right){\mathrm{G}}_{1\mathrm{m}-1}\left(\mathsf{ϵ}\right)\mathrm{d}\mathsf{ϵ} \\ {\mathrm{G}}_{2\mathrm{m}}\left(\mathrm{t}\right)={\mathrm{G}}_2\left(0\right)+{{\displaystyle \int }}_0^{\mathrm{t}}{\mathcal{H}}_6\left(\mathrm{t},\mathsf{ϵ}\right){\mathrm{G}}_{2\mathrm{m}-1}\left(\mathsf{ϵ}\right)\mathrm{d}\mathsf{ϵ} \\ {\mathrm{R}}_{1\mathrm{m}}\left(\mathrm{t}\right)={\mathrm{R}}_1\left(0\right)+{{\displaystyle \int }}_0^{\mathrm{t}}{\mathcal{H}}_7\left(\mathrm{t},\mathsf{ϵ}\right){\mathrm{R}}_{1\mathrm{m}-1}\left(\mathsf{ϵ}\right)\mathrm{d}\mathsf{ϵ} \\ { \mathrm{R}}_{2\mathrm{m}}\left(\mathrm{t}\right)={\mathrm{R}}_2\left(0\right)+{{\displaystyle \int }}_0^{\mathrm{t}}{\mathcal{H}}_8\left(\mathrm{t},\mathsf{ϵ}\right){\mathrm{R}}_{2\mathrm{m}-1}\left(\mathsf{ϵ}\right)\mathrm{d}\mathsf{ϵ} \end{gathered}$$

These Neumann series are convergent due to the Lipschtizian character of $\mathcal{H}’\mathrm{s}$, then:

$$\begin{gathered} \Vert {\mathrm{S}}_{1\mathrm{m}+1}-{\mathrm{S}}_{1\mathrm{m}}\Vert \le {{\displaystyle \int }}_0^{\mathrm{t}}\Vert {\mathcal{H}}_1\left(\mathrm{t},\mathsf{ϵ}\right){\mathrm{S}}_{1\mathrm{m}}\left(\mathsf{ϵ}\right)-{\mathcal{H}}_1\left(\mathrm{t},\mathsf{ϵ}\right),{\mathrm{S}}_{1\mathrm{m}-1}\left(\mathsf{ϵ}\right)\mathrm{d}\mathsf{ϵ}\Vert \\ \le {{\displaystyle \int }}_0^{\mathrm{t}}\Vert {\mathrm{S}}_{1\mathrm{m}}\left(\mathrm{t}\right)-{\mathrm{S}}_{1\mathrm{m}-1}\left(\mathrm{t}\right)\Vert \mathrm{d}\mathsf{ϵ} \le {\partial }_1\Vert {\mathrm{S}}_{1\mathrm{m}}\left(\mathrm{t}\right)-{\mathrm{S}}_{1\mathrm{m}-1}{\left(\mathrm{t}\right)\Vert }_{\infty } \end{gathered}$$

Other equations are given as follows:

$$\begin{gathered} \Vert {\mathrm{S}}_{2\mathrm{m}+1}-{\mathrm{S}}_{2\mathrm{m}}\Vert \le {\partial }_2\Vert {\mathrm{S}}_{2\mathrm{m}}\left(\mathrm{t}\right)-{\mathrm{S}}_{2\mathrm{m}-1}{\left(\mathrm{t}\right)\Vert }_{\infty } \\ \Vert {\mathrm{I}}_{1\mathrm{m}+1}-{\mathrm{I}}_{1\mathrm{m}}\Vert \le {\partial }_3\Vert {\mathrm{I}}_{1\mathrm{m}}\left(\mathrm{t}\right)-{\mathrm{I}}_{1\mathrm{m}-1}{\left(\mathrm{t}\right)\Vert }_{\infty } \\ \Vert {\mathrm{I}}_{2\mathrm{m}+1}-{\mathrm{I}}_{2\mathrm{m}}\Vert \le {\partial }_4\Vert {\mathrm{I}}_{2\mathrm{m}}\left(\mathrm{t}\right)-{\mathrm{I}}_{2\mathrm{m}-1}{\left(\mathrm{t}\right)\Vert }_{\infty } \\ \Vert {\mathrm{G}}_{1\mathrm{m}+1}-{\mathrm{G}}_{1\mathrm{m}}\Vert \le {\partial }_5\Vert {\mathrm{G}}_{1\mathrm{m}}\left(\mathrm{t}\right)-{\mathrm{G}}_{1\mathrm{m}-1}{\left(\mathrm{t}\right)\Vert }_{\infty } \\ \Vert {\mathrm{G}}_{2\mathrm{m}+1}-{\mathrm{G}}_{2\mathrm{m}}\Vert \le {\partial }_6\Vert {\mathrm{G}}_{2\mathrm{m}}\left(\mathrm{t}\right)-{\mathrm{G}}_{2\mathrm{m}-1}{\left(\mathrm{t}\right)\Vert }_{\infty } \\ \Vert {\mathrm{R}}_{1\mathrm{m}+1}-{\mathrm{R}}_{1\mathrm{m}}\Vert \le {\partial }_7\Vert {\mathrm{R}}_{1\mathrm{m}}\left(\mathrm{t}\right)-{\mathrm{R}}_{1\mathrm{m}-1}{\left(\mathrm{t}\right)\Vert }_{\infty } \\ \Vert {\mathrm{R}}_{2\mathrm{m}+1}-{\mathrm{R}}_{2\mathrm{m}}\Vert \le {\partial }_8\Vert {\mathrm{R}}_{2\mathrm{m}}\left(\mathrm{t}\right)-{\mathrm{R}}_{2\mathrm{m}-1}{\left(\mathrm{t}\right)\Vert }_{\infty } \end{gathered}$$

The above inequalities prove the existence of the function $\mathcal{H}.$ We now show the uniqueness of the solution by assuming that the kernels ${\mathcal{H}}_{\mathrm{i}}, \mathrm{i}=1,2,\dots , 8$ are separable, i.e., ${\mathcal{H}}_1\left(\mathrm{t},\mathsf{ϵ}\right)=\varnothing \left(\mathrm{t}\right)\mathsf{\sigma }\left(\mathsf{ϵ}\right)$. By denoting:

$$\mathsf{\zeta }\left(\mathrm{t}\right)={{\displaystyle \int }}_0^{\mathrm{t}}{\mathrm{S}}_1\left(\mathsf{ϵ}\right)\mathsf{\sigma }\left(\mathsf{ϵ}\right)\mathrm{d}\mathsf{ϵ}$$

then, because ${\mathrm{S}}_1\left(\mathrm{t}\right)={\mathrm{S}}_1\left(0\right)+{{\displaystyle \int }}_0^{\mathrm{t}}{\mathcal{H}}_1\left(\mathrm{t},\mathsf{ϵ}\right){\mathrm{S}}_1\left(\mathsf{ϵ}\right)\mathrm{d}\mathsf{ϵ}={\mathrm{S}}_1\left(0\right)+{{\displaystyle \int }}_0^{\mathrm{t}}\varnothing \left(\mathrm{t}\right)\mathsf{\sigma }\left(\mathsf{ϵ}\right){\mathrm{S}}_1\left(\mathsf{ϵ}\right)\mathrm{d}\mathsf{ϵ},$ we have:

$$\mathsf{\zeta } \prime \left(\mathrm{t}\right)=\varnothing \left(\mathrm{t}\right)\mathsf{\sigma }\left(\mathrm{t}\right)\mathsf{\zeta }\left(\mathrm{t}\right)+{ \mathrm{S}}_1\left(0\right)\mathsf{\sigma }\left(\mathrm{t}\right)$$

If $\mathsf{\sigma }\left(0\right)=0$, the solution $\mathsf{\zeta }\left(\mathrm{t}\right)$ follows $\mathsf{\zeta }\left(\mathrm{t}\right)={ \mathrm{e}}^{{{\displaystyle \int }}_0^{\mathrm{t}}\varnothing \left(\mathsf{ϵ}\right)\mathsf{\sigma }\left(\mathsf{ϵ}\right)\mathrm{d}\mathsf{ϵ}}$ and from the definition of $\mathsf{\zeta }\left(\mathrm{t}\right)$, the solution of ${\mathrm{S}}_1\left(\mathrm{t}\right)$ is given as:

$${\mathrm{S}}_1\left(\mathrm{t}\right)={\mathrm{S}}_1\left(0\right)+\varnothing \left(\mathrm{t}\right){{\displaystyle \int }}_0^{\mathrm{t}}{\mathrm{S}}_1\left(\mathsf{ϵ}\right)\mathsf{\sigma }\left(\mathsf{ϵ}\right)\mathrm{d}\mathsf{ϵ}={ \mathrm{S}}_1\left(0\right)+\varnothing \left(\mathrm{t}\right)\mathsf{\zeta }\left(\mathrm{t}\right)={\mathrm{S}}_1\left(0\right)+\varnothing \left(\mathrm{t}\right){\mathrm{e}}^{{{\displaystyle \int }}_0^{\mathrm{t}}\varnothing \left(\mathsf{ϵ}\right)\mathsf{\sigma }\left(\mathsf{ϵ}\right)\mathrm{d}\mathsf{ϵ}}$$

Hence, by the unicity of the solution $\mathsf{\zeta }\left(\mathrm{t}\right)$, there exists just only one continuous solution for ${\mathrm{S}}_1\left(\mathrm{t}\right)$. Following the same approach, we can obtain a unique solution for the remaining equations of the system (2).

4.6. Basic Reproduction Number

In this section, we apply the idea of a next-generation matrix by linearizing the model Equation (1) near the endemic state, for the infectious part of the system (1), and we obtain:

$$\begin{gathered} {\mathrm{dI}}_1/\mathrm{dt}={\mathsf{\phi }}_1\left({\mathrm{I}}_1+{\mathrm{I}}_2\right){\mathrm{S}}_1{}^{**}-({\mathrm{a}}_1^{\mathrm{I}}+{\mathsf{\eta }}_1){\mathrm{I}}_1 , \\ {\mathrm{dI}}_2/\mathrm{dt}={\mathsf{\phi }}_2\left({\mathrm{I}}_1+{\mathrm{I}}_2\right){\mathrm{S}}_2{}^{**}+{\mathrm{a}}_1^{\mathrm{I}}{\mathrm{I}}_1-{\mathsf{\eta }}_2{\mathrm{I}}_{2 }- {\mathsf{\mu }}_2{}^{\mathrm{I}}{\mathrm{I}}_2, \\ {\mathrm{dG}}_1/\mathrm{dt}={ \mathsf{\eta }}_1{\mathrm{I}}_1+{ \mathsf{\theta }}_1^{\mathrm{G}}{\mathrm{S}}_1{}^{**} – ({\mathrm{a}}_1^{\mathrm{G}}+{\mathsf{\gamma }}_1+{\mathsf{ϵ}}_1){\mathrm{G}}_1, \\ {\mathrm{dG}}_2/\mathrm{dt}={ \mathsf{\eta }}_2{\mathrm{I}}_2+{ \mathrm{a}}_1^{\mathrm{G}}{\mathrm{G}}_1+{ \mathsf{\theta }}_2^{\mathrm{G}}{\mathrm{S}}_2{}^{**} – ({\mathsf{\gamma }}_2+{ \mathsf{\mu }}_2{}^{\mathrm{G}}+{\mathsf{ϵ}}_2+{\mathsf{\beta }}_{22}){\mathrm{G}}_2 \end{gathered}$$

(4)

By summing Equation (4) and by denoting I as the size of all infectious, we have:

I(t) = I₁(t) + I₂(t) + G₁(t) + G₂(t), and

dI/dt = ${\mathsf{\phi }}_1\left({\mathrm{I}}_1+{\mathrm{I}}_2\right){\mathrm{S}}_1{}^{**}+{\mathsf{\phi }}_2\left({\mathrm{I}}_1+{\mathrm{I}}_2\right){\mathrm{S}}_2{}^{**}$− ${\mathsf{\mu }}_2{}^{\mathrm{I}}$ I₂$+{ \mathsf{\theta }}_1^{\mathrm{G}}{\mathrm{S}}_1{}^{**}-({\mathsf{\gamma }}_1+{\mathsf{ϵ}}_1){\mathrm{G}}_{1 }-({\mathsf{\gamma }}_2+{\mathsf{ϵ}}_2+{\mathsf{\beta }}_{22}){\mathrm{G}}_2$

The matrix J of the linearized system near a state (S₁, S₂) is:

$$\mathrm{J}=\left[\begin{array}{cccc}{\mathsf{\phi }}_1{\mathrm{S}}_1-{\mathrm{a}}_1^{\mathrm{I}}-{\mathsf{\eta }}_1& {\mathsf{\phi }}_1{\mathrm{S}}_1& 0& 0\\ {\mathsf{\phi }}_2{\mathrm{S}}_2+{\mathrm{a}}_1^{\mathrm{I}}& {\mathsf{\phi }}_2{\mathrm{S}}_2-{\mathsf{\eta }}_2-{ \mathsf{\mu }}_2{}^{\mathrm{I}}& 0& 0\\ {\mathsf{\eta }}_1& 0& -{\mathrm{a}}_1^{\mathrm{G}}-{\mathsf{\gamma }}_1-{\mathsf{ϵ}}_1& 0\\ 0& {\mathsf{\eta }}_2& {\mathrm{a}}_1^{\mathrm{G}}& -{\mathsf{\gamma }}_2-{ \mathsf{\mu }}_2{}^{\mathrm{G}}-{\mathsf{ϵ}}_2-{\mathsf{\beta }}_{22}\end{array}\right]$$

The corresponding characteristic polynomial P_J($\mathsf{\lambda }$) is equal to:

$$\begin{gathered} ({\mathsf{\phi }}_1{\mathrm{S}}_1-{\mathrm{a}}_1^{\mathrm{I}}-{\mathsf{\eta }}_1- \mathsf{\lambda })({\mathsf{\phi }}_2{\mathrm{S}}_2-{\mathsf{\eta }}_2-{ \mathsf{\mu }}_2{}^{\mathrm{I}}- \mathsf{\lambda }) \left(-{\mathrm{a}}_1^{\mathrm{G}}-{\mathsf{\gamma }}_1-{\mathsf{ϵ}}_1- \mathsf{\lambda }\right) \left(-{\mu }_2^{\mathrm{G}}-{\mathsf{\gamma }}_2-{\mathsf{ϵ}}_2-{\mathsf{\beta }}_{22}- \mathsf{\lambda }\right) \\ -{ \mathsf{\phi }}_1{\mathrm{S}}_1\left({\mathsf{\phi }}_2{\mathrm{S}}_2+{\mathrm{a}}_1^{\mathrm{I}}\right)\left(-{\mathrm{a}}_1^{\mathrm{G}}-{\mathsf{\gamma }}_1-{\mathsf{ϵ}}_1- \mathsf{\lambda }\right)\left(-{\mu }_2^{\mathrm{G}}-{\mathsf{\gamma }}_2-{\mathsf{ϵ}}_2-{\mathsf{\beta }}_{22}-\mathsf{\lambda }\right)= {\mathrm{P}}_{\mathrm{J}}\left(\mathsf{\lambda }\right) \end{gathered}$$

The positive eigenvalues of J are roots of the following polynomial:

$${\mathsf{\lambda }}^2-\left({\mathsf{\phi }}_1{\mathrm{S}}_1-{\mathrm{a}}_1^{\mathrm{I}}-{\mathsf{\eta }}_1+{\mathsf{\phi }}_2{\mathrm{S}}_2-{\mathsf{\eta }}_2-{ \mathsf{\mu }}_2{}^{\mathrm{I}}\right)\mathsf{\lambda }-{ \mathsf{\phi }}_1{\mathrm{S}}_1\left({\mathsf{\phi }}_2{\mathrm{S}}_2+{\mathrm{a}}_1^{\mathrm{I}}\right)+({\mathsf{\phi }}_1{\mathrm{S}}_1-{\mathrm{a}}_1^{\mathrm{I}}-{\mathsf{\eta }}_1)({\mathsf{\phi }}_2{\mathrm{S}}_2-{\mathsf{\eta }}_2-{ \mathsf{\mu }}_2{}^{\mathrm{I}})$$

These roots are equal to: B ± (B² − C)^1/2, where B and C are equal to:

$$\begin{gathered} \mathrm{B}=\frac{{\mathsf{\phi }}_1{\mathrm{S}}_1-{\mathrm{a}}_1^{\mathrm{I}}-{\mathsf{\eta }}_1+{\mathsf{\phi }}_2{\mathrm{S}}_2-{\mathsf{\eta }}_2-{ \mathsf{\mu }}_2{}^{\mathrm{I}}}{2} \\ \mathrm{C}=-{ \mathsf{\phi }}_1{\mathrm{S}}_1\left({\mathsf{\phi }}_2{\mathrm{S}}_2+{\mathrm{a}}_1^{\mathrm{I}}\right)+({\mathsf{\phi }}_1{\mathrm{S}}_1-{\mathrm{a}}_1^{\mathrm{I}}-{\mathsf{\eta }}_1)({\mathsf{\phi }}_2{\mathrm{S}}_2-{\mathsf{\eta }}_2-{ \mathsf{\mu }}_2{}^{\mathrm{I}}) \end{gathered}$$

Then, we have:

$$\mathrm{B}+({\mathrm{B}}^2-\mathrm{C}{)}^{1/2}=\frac{{\mathsf{\phi }}_1{\mathrm{S}}_1-{\mathrm{a}}_1^{\mathrm{I}}-{\mathsf{\eta }}_1+{\mathsf{\phi }}_2{\mathrm{S}}_2-{\mathsf{\eta }}_2-{ \mathsf{\mu }}_2{}^{\mathrm{I}}}{2}+[\frac{{({\mathsf{\phi }}_1{\mathrm{S}}_1-{\mathrm{a}}_1^{\mathrm{I}}-{\mathsf{\eta }}_1+{\mathsf{\phi }}_2{\mathrm{S}}_2-{\mathsf{\eta }}_2-{ \mathsf{\mu }}_2{}^{\mathrm{I}})}^2}{4}+{ \mathsf{\phi }}_1{\mathrm{S}}_1\left({\mathsf{\phi }}_2{\mathrm{S}}_2+{\mathrm{a}}_1^{\mathrm{I}}\right)-({\mathsf{\phi }}_1{\mathrm{S}}_1-{\mathrm{a}}_1^{\mathrm{I}}-{\mathsf{\eta }}_1)({\mathsf{\phi }}_2{\mathrm{S}}_2-{\mathsf{\eta }}_2-{ \mathsf{\mu }}_2{}^{\mathrm{I}}){]}^{1/2}$$

Hence, near unstable endemic state the positive dominant eigenvalue $\mathsf{\Lambda }$ is equal to:

$$\mathsf{\Lambda }=\frac{{\mathsf{\phi }}_1{\mathrm{S}}_1-{\mathrm{a}}_1^{\mathrm{I}}-{\mathsf{\eta }}_1+{\mathsf{\phi }}_2{\mathrm{S}}_2-{\mathsf{\eta }}_2-{ \mathsf{\mu }}_2{}^{\mathrm{I}}}{2}+{[\frac{{({\mathsf{\phi }}_1{\mathrm{S}}_1-{\mathrm{a}}_1^{\mathrm{I}}-{\mathsf{\eta }}_1-{\mathsf{\phi }}_2{\mathrm{S}}_2+{\mathsf{\eta }}_2+{ \mathsf{\mu }}_2{}^{\mathrm{I}})}^2}{4}+{ \mathsf{\phi }}_1{\mathrm{S}}_1\left({\mathsf{\phi }}_2{\mathrm{S}}_2+{\mathrm{a}}_1^{\mathrm{I}}\right)}^{}{]}^{1/2}$$

Therefore, the basic reproduction number R₀ equal to $\mathsf{\Lambda }$ near the endemic stationary state depends mainly on the infection rates ${\mathsf{\phi }}_1$ and ${ \mathsf{\phi }}_2$, when the sizes S₁ and S₂ are sufficiently important. If after a change of parameter values, the endemic state becomes unstable and R₀ is becoming more than 1, then an epidemic wave starts.

5. Numerical Simulation of the Model: Some Examples of the COVID-19 Outbreak in Kuwait, France and Cameroon

First, we provide some explanations for the data used for the simulations, including how we assumed some of the parameters, calculated others, and selected some from the literature cited in this article. In order to determine the susceptible classes sizes, we used Wikipedia data on the three studied countries to obtain their total population size. We then calculated the ratio between the young and elderly from the data presented in Figure 2, Figure 3 and Figure 4 and used this ratio to determine the sizes value for the two susceptible classes at the exponential phase considered (for Kuwait 28 December 2020, France 30 October 2021, and Cameroon 19 September 2021), because the Ross and Kermack-McKendrick model is only suited for the exponential growth phase of an epidemic wave. Because the progression rates to the classes reversed recovery and fully recovered in Figure 5 are not zero, we assumed that they should have some populations at the start of the wave, even if they were small. The transmission rate for Cameroon was chosen at the start of the second wave in January 2021 from [56,57], for France at the start of the fifth wave in December 2021 [58], and for Kuwait at the start of the fourth wave in December 2021 [59]. The values of specific death rates due to disease and natural death rates for Kuwait and France were taken from [58,59,60]. For Cameroon, the natural death rate was taken from [61], while the specific death rate due to disease was calculated from the cumulated deaths number due to the disease in two years divided by the cumulated infectious number in these two years [61]. The vaccination rate was chosen from [62], while the loss of resistance was chosen from [63]. Other parameters were assumed.

5.1. COVID-19 Outbreak in Kuwait

We have chosen the following set of parameter values corresponding to the COVID-19 outbreak at the start of the fourth wave in December 2021 in Kuwait:

$$\begin{gathered} {\mathsf{\phi }}_{1 }=1.7,{\mathsf{\phi }}_{2 }=0.9,{\mathsf{\mu }}_2{}^{\mathrm{S}}=0.003,{\mathsf{ϵ}}_{1 }=0.28,{\mathsf{\beta }}_{22 }=0,{ \mathsf{\mu }}_2{}^{\mathrm{I}}=0.0025,{ \mathsf{\mu }}_2{}^{\mathrm{G}}= \\ 0.002,{ \mathsf{\mu }}_2{}^{\mathrm{R}}=0.0021 ,{ \mathsf{\beta }}_{11 }=2.1 ,{\mathsf{\beta }}_{21 }=2.3 ,{ \mathsf{\theta }}_1^{\mathrm{G}} =0.765,{\mathsf{\theta }}_1^{\mathrm{R}}=0.678,{\mathsf{\gamma }}_1= \\ 0.62,{\mathsf{\gamma }}_2=0.74,{ \mathsf{\theta }}_2^{\mathrm{G}} =0.45 ,{\mathsf{\theta }}_2^{\mathrm{R}}=0.33,{ \mathrm{a}}_1^{\mathrm{S}}=0.7,{\mathrm{a}}_1^{\mathrm{I}}=0.54,{\mathrm{a}}_1^{\mathrm{G}}=0.65,{\mathrm{a}}_1^{\mathrm{R}}= \\ 0.8, {\mathsf{ϵ}}_{2 }=0.19 ,{\mathsf{\eta }}_1=0.3{ \mathrm{and} \mathsf{\eta }}_2=0.38. \end{gathered}$$

with initial values for S₁ = 4413099, S₂ = 51422, I₁ = 65, I₂ = 194, G₁ = 20, G₂ = 17, R₁ = 30 and R₂ = 73.

We present the visualization results for the simulated values in Figure 6.

5.2. COVID-19 Outbreak in France

We have chosen the following set of parameter values corresponding to the COVID-19 outbreak at the start of the fifth wave in December 2021 in France:

$$\begin{gathered} {\mathsf{\phi }}_{1 }= 1.2,{\mathsf{\phi }}_{2 }= 0.9,{\mathsf{\mu }}_2{}^{\mathrm{S}}= 0.009,{\mathsf{ϵ}}_{1 }= 0.28,{\mathsf{\beta }}_{22 }= 0,{ \mathsf{\mu }}_2{}^{\mathrm{I}}= 0.0025,{ \mathsf{\mu }}_2{}^{\mathrm{G}}= \\ 0.002,{ \mathsf{\mu }}_2{}^{\mathrm{R}}=0.0021 ,{ \mathsf{\beta }}_{11 }= 1.9,{\mathsf{\beta }}_{21 }= 2.2,{ \mathsf{\theta }}_1^{\mathrm{G}} = 0.735,{\mathsf{\theta }}_1^{\mathrm{R}}= 0.678,{\mathsf{\gamma }}_1= \\ 0.62,{\mathsf{\gamma }}_2=0.74,{ \mathsf{\theta }}_2^{\mathrm{G}} =0.45,{\mathsf{\theta }}_2^{\mathrm{R}}=0.33,{ \mathrm{a}}_1^{\mathrm{S}}= 0.7,{\mathrm{a}}_1^{\mathrm{I}}= 0.54,{\mathrm{a}}_1^{\mathrm{G}}= 0.65,{\mathrm{a}}_1^{\mathrm{R}}= \\ 0.8, {\mathsf{ϵ}}_{2 }= 0.19,{\mathsf{\eta }}_1= 0.3{ \mathrm{and} \mathsf{\eta }}_2= 0.38. \end{gathered}$$

with initial values for S₁ = 53372880, S₂ = 14017120, I₁ = 3912, I₂ = 1757, G₁ = 200, G₂ = 170, R₁ = 300 and R₂ = 730.

We present the visualization results for the simulated values in Figure 7.

5.3. COVID-19 Outbreak in Cameroon

We choose the following set of parameter values corresponding to the COVID-19 outbreak at the start of the second wave in January 2021 in Cameroon:

$$\begin{gathered} {\mathsf{\phi }}_{1 }=1.3,{\mathsf{\phi }}_{2 }=0.9,{\mathsf{\mu }}_2{}^{\mathrm{S}}=0.009,{\mathsf{ϵ}}_{1 }=0.28,{\mathsf{\beta }}_{22 }=0,{ \mathsf{\mu }}_2{}^{\mathrm{I}}=0.0025,{ \mathsf{\mu }}_2{}^{\mathrm{G}}= \\ 0.002,{ \mathsf{\mu }}_2{}^{\mathrm{R}}=0.0021 ,{ \mathsf{\beta }}_{11 }=4.5,{\mathsf{\beta }}_{21 }=4.1,{ \mathsf{\theta }}_1^{\mathrm{G}} =0.024,{\mathsf{\theta }}_1^{\mathrm{R}}=0.678,{\mathsf{\gamma }}_1= \\ 0.62,{\mathsf{\gamma }}_2=0.74,{ \mathsf{\theta }}_2^{\mathrm{G}} =0.45, {\mathsf{\theta }}_2^{\mathrm{R}}=0.33,{ \mathrm{a}}_1^{\mathrm{S}}=0.7,{\mathrm{a}}_1^{\mathrm{I}}=0.54,{\mathrm{a}}_1^{\mathrm{G}}=0.65,{\mathrm{a}}_1^{\mathrm{R}}= \\ 0.8, {\mathsf{ϵ}}_{2 }=0.19,{\mathsf{\eta }}_1=0.3{ \mathrm{and} \mathsf{\eta }}_2=0.38. \end{gathered}$$

with initial values for S₁ = 25828973, S₂ = 721027, I₁ = 151, I₂ = 21, G₁ = 75, G₂ = 33, R₁ = 56 and R₂ = 64.

We present the visualization results for the simulated values in Figure 8.

6. Discussion

The results for the three considered countries present some similarities but also some differences that we will discuss in the following.

First, in each case, the exponential growth of the infectious I₁ + I₂ and completely recovering R₁ + R₂ populations sizes correspond roughly to the data given in [57]. The model simulations give more, i.e., allows to see the part brought by each age class to the global growth.

Second, for Kuwait, Figure 6 shows a faster growth of infectious (I) in the young class (≤65 years) than in the older class (>65 years), which is also the case for France (Figure 7) and for Cameroon (Figure 8). On the other hand, this phenomenon is reversed in the three countries, with regard to the growth of populations immunized in a transient manner (G) and in a lasting manner (R). The phenomenon is more marked for Kuwait than for France, France itself having a more marked difference than for Cameroon. This is partly explained by the better vaccination rate of Kuwait than that of France and Cameroon in the young class, the effectiveness of the vaccination (which decreases in the older class) having been assumed to be equal for the three countries.

7. Conclusions

Taking age into account in modeling the COVID-19 pandemic makes it possible to simulate the differential dynamic behavior of the growth of infectious and immune populations, young and old, in order, for example, to adjust the vaccine policy according to the age. Future work should, for example, take into account more age groups, at least four: children (age ≤ 12 years), adolescents (12 years < age ≤ 18 years), young adults (18 years < age ≤ 65 years), and older adults (>65 years). The most important pitfall in simulations of such a model is the estimation (by observation, calculation, or assumption) of its parameters, already difficult with two age groups. A random choice of the values of the parameters in plausible intervals followed by a study of the sensitivity to the parameters of the model, could make it possible to partly overcome the constraint of parameter estimation in a more precise future epidemiological–demographic model, in particular with regard to age groups.

Another perspective would consist in introducing, thanks to Usher’s model [41,42], an accelerated aging due to infection, as well as an influence of exogenous determinants such as geo-climatic, socio-economic, and health-related factors [64,65,66,67,68,69,70,71,72,73,74,75,76], which weight differently on the different age groups, therefore, changing the growth dynamics specific to each of the sub-populations studied in this article. These perspectives will be examined in future works.