Global investigation for an "SIS" model for COVID-19 epidemic with asymptomatic infection

1.   Introduction

COVID-19 is one of the three Coronaviruses that has caused epidemic outbreaks over the last two decades. It can spread through close contact, coughing, sneezing, or talking. For COVID-19, there are two types of infected individuals: one is symptomatic infected individuals, defined as those who show symptoms (such as fever, cough, sore throat, etc.) after obtaining infected. The other is asymptomatic infected individuals, defined as those who do not show symptoms after obtaining infected ^[1,2]. Asymptomatic cases are not confirmed cases which are divided into two different states. The first one is the asymptomatic individuals who show no symptoms for the whole time of the infection. The second one is the asymptomatic individuals who exhibit symptoms after a period of the asymptomatic infection. In the generally tested subgroup, the proportion of asymptomatic individuals who are found to be positive for COVID-19 is alarmingly high. Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) shows that asymptomatic and mildly symptomatic infections may be the key to disease transmission ^[3,4]. Most asymptomatic infected individuals do not seek medical help because they have no obvious clinical signs which leads to the rapid spread of the disease. Therefore, the prevention and control of this specific type of patient on a global scale is a huge challenge and requires more attention from the world.

The global health crisis of the Coronavirus Covid-19 has brought to light the role of mathematical modeling in political and health decision-making ^[5,6,7,8]. We will study the SIS model inspired by the current Coronavirus outbreak. Mathematical models of infectious diseases have been used for over a decade to study and understand the mechanism of disease spread, predict the future of epidemics, and determine the best treatment strategy to protect human health. The famous SIR and SEIR mathematical epidemic models are the most used models providing good descriptions of infectious diseases especially when taking account of either symptomatic and asymptomatic infectious or time delay ^[9,10,11,12].

In this work, we have designed a compartmental ODE model to investigate the COVID-19 spreading dynamics incorporating non-linear saturated incidence rates for both asymptomatic and symptomatic infection. The incidence rate is an essential component for infectious disease transmission in a mathematical model for both symptomatic and asymptomatic infectious ^[13]. In the epidemiology perception, incidence rate measures the number of new cases of a disease within a specified time duration as a percentage of the number of persons at threat for the disease ^[14]. We have also included two control functions designing quarantine strategies which is used to reduce the contact among infected individuals. Here we have studied the qualitative nature of the model through investigation of local and global nature of equilibrium points. Finally we have characterised optimal control to reduce the outbreak size and the control implementation costs. Results of the optimal control model show that optimal control provides significant reduction in COVID-19 spread.

The paper is organised as follows: In Section 2, we have formulated the proposed model and have studied the boundedness and positivity of solutions. Steady state analysis with constant controls have been done in Section 3 with respect to the basic reproduction number ($\mathcal{R}$). In Section 4, we characterised the optimal control model and we give an efficiency analysis. Finally, some numerical simulations and discussions were given in Section 5. The conclusions of our study has been reported in Section 6.

2.   Model development

In this paper, we have studied a compartmental ODE model to investigate the COVID-19 spreading dynamics using non-linear saturated incidences. For this purpose we have formulated an ODE model by splitting the total host population $\Lambda$ into three disjoint epidemiological classes specifically susceptible individuals $S(t)$, asymptomatic infectious individuals $I_a(t)$ and symptomatic infectious individuals $I_s(t)$ i.e. $\Lambda = S(t)+I_a(t)+I_s(t)$.

To formulate the model we have considered the following factors:

● Susceptible individual once infected has a probability $0 < p < 1$ to be symptomatic and then a probability $0 < 1-p < 1$ to be asymptomatic.

● Susceptible individual is infected by asymptomatic infected individuals at a nonlinear increasing rate $f_a(I_a)S$ and by symptomatic infected individuals at a nonlinear increasing rate $f_s(I_s)S$. The importance of these increasing incidence rates is that the rate of effective contacts between infected individuals and susceptible individuals increases with the increase of number of infective individuals.

● asymptomatic and symptomatic individuals are recovered at constant rates $\sigma_a$ and $\sigma_s$, respectively.

● Recovered individuals can catch the diseases and then they are added to susceptible compartment.

● asymptomatic individuals become symptomatic individuals at a constant rate $\gamma$.

● Susceptible individuals are recruited at constant rate $m\Lambda$, natural mortality rate is $m$.

According to the epidemiological assumptions, we proposed the following mathematical model

with positive initial condition $(S^0, I_a^0, I_s^0) \in \mathbb{R}_+^3\;$.

3.   Mathematical analysis

Assumption 1. $f_a$ and $f_s$ are bounded, non-negative $C^1(\mathbb{R}_+)$, concave and increasing functions with $f_a(0) = f_s(0) = 0$.

Lemma 1. The general non-linear incidence rates $f_a$ and $f_s$ satisfy $f_a'(I) I \leq f_a(I) \leq f_a'(0) I$ and $f_s'(I) I \leq f_s(I) \leq f_s'(0) I, \; \forall I > 0$.

Proof. Let $I, I_1 \in \mathbb{R}_+$, and the function $g_1(I) = f_a(I)-I f_a'(I)$. Since $f_a'(I) \geq 0$ ($f_a$ is increasing function) and $f_a''(I) \leq 0$ ($f_a$ is concave) then $g_1'(I) = -If_a''(I)\geq 0$ and $g_1(I) \geq g_1(0) = 0$. Therefore, $f_a(I)\geq If_a'(I)$. Similarly, let $g_2(I) = f_a(I)-If_a'(0)$ then $g_2'(I) = f_a'(I)-f_a'(0)\leq 0$ once $f_a$ is a concave function. Thus $g_2(I) \leq g_2(0) = 0$ and $f_a(I)\leq If_a'(0)$. Similarly for the function $f_s$. It is necessary that the state variables $S (t), I_a(t)$ and $I_s(t)$ remain non-negative for all $t \geq 0$.

Proposition 1. $\Omega_0 = \left\{(S, I_a, I_s)\in \mathbb{R}_+^3\; /\; S+I_a+I_s = \Lambda \right\}$ is a positively invariant compact set for model (2.1).

Proof. Assume that the initial condition $(S^0, I_a^0, I_s^0) \in \mathbb{R}_+^3\;$. If $S = 0$ then $\dot S = m \Lambda + \sigma_a I_a + \sigma_s I_s > 0$ therefore $S(t) > 0$ for all $t > 0$. Similarly, assume that $I_a = 0$ then $\dot I_a = (1-p) Sf_s(I_s) \geq 0$ therefore $I_a(t)\geq 0$ for all $t > 0$. By the same way, assume that $I_s = 0$ therefore $\dot I_s = p Sf_a(I_a) + \gamma I_a \geq 0$ then $I_s(t)\geq 0$ for all $t > 0$. Therefore the model (3.2) admits a non-negative solution.

By summing the equations of (2.1), we get, for $T = S + I_a + I_s - \Lambda$, the following equation :

Hence

Hence, $\Omega_0$ is invariant for the model (2.1) due to all variables are non-negative. Using the conservation principles, we can compute $S$ as function of $I_a$ and $I_s$ :

Now, we can reduce the analysis of the original system (2.1) to the analysis of the following two-dimensional limiting system on the invariant set $\Omega$ :

with positive initial condition $(I_a^0, I_s^0) \in \mathbb{R}_+^2\;$.

$\Omega = \left\{(I_a, I_s)\in \mathbb{R}_+^2\; /\; I_a+I_s \leq \Lambda \right\}$ is a positively invariant compact set for model (3.2). It is obvious that $E_0 = (0, 0)$ is the only disease-free equilibrium of system (3.2).

The global behavior of our system inevitably depends on the basic reproduction number ($\mathcal{R}$), that is, the average number of secondary cases produced by an infectious individual who is introduced into an established population only of susceptible.

To drive the basic reproduction number ($\mathcal{R}$) for complex compartmental models, we use the next-generation operator approach proposed by Diekmann et al. ^[15,16]. Following the approach of van den Driessche and Watmough, we can rewrite (3.2)

where $\mathcal{F}$ denotes the rate of appearance of new infections, and $\mathcal{V}$ denotes the rate of transfer of individuals into or out of each population set. Furthermore,

Since det$(FV^{-1}) = 0$, then one of the eigenvalues is zero. Therefore, we can deduce the second eigenvalue since their sum is equal to the trace of the matrix $FV^{-1}$. Then, the basic reproduction number which is the spectral radius of $FV^{-1}$ (the maximum of the absolute values of its eigenvalues) is given by

Proposition 2. The model (3.2) admits a unique disease-free equilibrium $E_0 = \left(0, 0\right)$ and a unique endemic equilibrium $E^* = (I_a^*, I_s^*)$ such that $I_a^*, I_s^* > 0$.

Proof. Equilibria of (3.2) satisfy

which reduces to

Thus

We conclude that from (3.6)

We can write this equation on the form

where the function $g$ given by

● If $I_a = 0$, then $I_s = 0$. This equilibrium named the disease-free equilibrium will be noted here by $E_0 = (0, 0)$.

● If $I_a\not = 0$, then $g(I_a) = 0$. Let us calculate the derivative of the function $g$ given by

The functions $f_a$ and $f_s$ satisfy $f_a'(I)I \leq f_a(I)$ and $f_s'(I)I \leq f_s(I), \; \forall I \geq 0$ and all variables are non-negative, hence $g$ is an increasing function and satisfies $g'(I_a) < 0$.

An easy calculation gives

and

Since $\mathcal{R} > 1$, $g(\Lambda) < 0$ and $\lim_{I_a \to \, 0}g(I_a) < 0$, then $g(I_a) = 0$ has a unique positive solution $I_a^*$ in $(0, \Lambda)$ and therefore the equilibrium state $E^* = (I_a^*, I_s^*)$ is unique with $I_s^* = \dfrac{\gamma+p(m +\sigma_a)}{(1-p)(m +\sigma_s)} I_a^*$.

3.1.   Local analysis

In this subsection, the local stability behaviours of equilibria are discussed.

Theorem 1. $E_0$ is LAS when $\mathcal{R} < 1$ and it is unstable when $\mathcal{R} > 1$.

Proof. The Jacobian matrix given at the equilibrium point $E_0$ is

The trace is given by

Since $\mathcal{R} = \dfrac{p \Lambda f_s'(0)}{(m +\sigma_s)} + \dfrac{(1-p)\Lambda f_a'(0)}{(\gamma +m +\sigma_a)} + \dfrac{(1-p) \gamma \Lambda f_s'(0) }{(\gamma +m +\sigma_a)(m +\sigma_s)}$ then when $\mathcal{R} < 1$, it is obvious that $\dfrac{p \Lambda f_s'(0)}{(m +\sigma_s)} < 1$, $\dfrac{(1-p)\Lambda f_a'(0)}{(\gamma +m +\sigma_a)} < 1$ therefore Trace $(J_0) < 0$. The determinant is given by

If $\mathcal{R} < 1$, then we have negative eigenvalues. Therefore, $E_0$ is LAS. Whereas, if $\mathcal{R} > 1$, $E_0$ is therefore unstable.

Theorem 2. $E^*$ is LAS when $\mathcal{R} > 1$.

Proof. For the equilibrium point $E^*$, the Jacobian is given by $J^* = \left(\begin{array}{ccccccc} a_{11} & a_{12}\\ a_{21} & a_{22} \end{array} \right)$ where

The trace is given by

and the determinant is given by

By Lemma 1, the functions $f_a$ and $f_s$ satisfy $f_a'(I)I \leq f_a(I)$ and $f_s'(I)I \leq f_s(I), \; \forall I \geq 0$ and since $(\Lambda - I_a^* - I_s^*) > 0$, then, all terms of Det$(J^*)$ are either positive or non-negative and therefore Det$(J^*) > 0$ and the equilibrium point $E^*$ (which exists only if $\mathcal{R} > 1$) is LAS.

3.2.   Global analysis

Here, we discuss the global behaviour of the equilibrium points. For proving the global stability of the positive equilibrium points $E_0$ and $E^*$, we will exclude the existence of periodic solutions of system (3.2) by using generalized Bendixson-Dulac theorem. We first prove that the system (3.2) has no periodic orbit nor poly-cycle on $\Omega$

Theorem 3. System (3.2) has no periodic orbits nor poly-cycles on $\Omega$.

Proof. Consider a solution of system (3.2) on ${\Omega }$. Let $\xi_1 = \ln(I_a)$ and $\xi_2 = \ln(I_s)$ then the transformation of (3.2) gives the following new system:

We have

Thus using Dulac criterion ^[17,18], system (3.8) has no periodic trajectory. Therefore system (3.2) has no periodic orbit inside ${\mathcal \Omega }$.

Theorem 4. The solution of (3.2) converges asymptotically to :

● $E_0$ if $\mathcal{R} < 1$.

● $E^*$ if $\mathcal{R} > 1$.

Proof. We restrict the proof to the case where $\mathcal{R} > 1$. The other case can be done similarly. The system (3.2) admits two equilibrium points $E_0$ and $E^*$. $E_0$ is an unstable node, and only $E^*$ is a stable node. We aim to prove that $E^*$ is globally asymptotically stable. Let $I_a(0) > 0, I_s(0) > 0$ and $\omega$ the $\omega$-limit set of $(I_a(0), I_s(0))$. $\omega$ is an invariant compact set and $\omega\subset\bar{\Omega }$. $M$ can't be $E_0$ because $E_0$ is an unstable node and can't be a part of the $\omega$-limit set of $(I_a(0), I_s(0))$. System (3.2) has no periodic orbit inside $\Omega$. Using the Poincaré-Bendixon Theorem ^[17,18], $E^*$ is a globally asymptotically stable equilibrium point for system (3.2).

4.   Optimal control problem

The necessary conditions that an optimal pair must satisfy come from Pontryagin's Maximum Principle. Our goal is to minimize the number of infected individuals either asymptomatic or symptomatic while keeping those corresponding cost of the quarantine strategies low during the epidemic. Thus, we define a control function set as $V = \{v_1, v_2\}$ where $v_1(t)$ and $v_1(t)$ are the control variables for the quarantine strategies measuring the effort to reduce the contact between susceptible individuals and both kind of infected individuals (asymptomatic and symptomatic), respectively.

Therefore the model (2.1) is modified to the following new model:

Several optimal strategies for epidemic models were proposed in several previous works ^{[19,20,21,22]}. In our case, we discussed an optimal strategy that reduces the contact between infected and uninfected individuals to optimise the number of infective individuals.

Recall that the set $\Omega = \{(S, I_a, I_s)\in \mathbb{R}_+^3\; /\; S+I_a+I_s = \Lambda \}$ is a positively invariant compact set for (4.1).

Assume that $f_a$ and $f_s$ are globally Lipschitz with Lipschitz constants $L_a$ and $L_s$ where $\bar f_a = \sup_{I > 0} f_a(I)$ and $\bar f_s = \sup_{I > 0} f_s(I)$. Define the space

where $\bar v_1$ and $\bar v_1$ are two non-negative constants.

Our goal is to control the number infected individuals and minimize the cost of the effort to reduce the contact between susceptible and infected individuals $(v_1, v_2)$. The optimal control problem considered here is by minimizing an objective functional $J(v_1, v_2)$ as following

$\beta_1$ and $\beta_2$ are two constants.

Let $\varphi = \left(S, I_a, I_s\right)^t$, then we express (4.1) as

where $B = \left(\begin{array}{cccccc} -m & \sigma_a & \sigma_s\\ 0 & -(m+\sigma_a+\gamma) & 0 \\ 0 & \gamma & -(m+\sigma_s) \end{array} \right)$ and $\psi(\varphi) = \left(\begin{array}{ccc} m\Lambda-S((1-v_1)f_a(I_a)+(1-v_2)f_s(I_s))\\ (1-p) S((1-v_1)f_a(I_a)+(1-v_2)f_s(I_s)) \\ p S((1-v_1)f_a(I_a)+(1-v_2)f_s(I_s)) \end{array} \right)$.

Proposition 3. $Z(\varphi)$ is a Lipschitz continuous function.

Proof. First we shall prove that the function $\psi$ is uniformly Lipschitz and continuous since

where $M = 2 \max\left(2(1- \bar v_1)\Lambda L_a, 2(1- \bar v_1) \bar f_a, 2 (1- \bar v_2) \Lambda L_s, 2 (1- \bar v_2) \bar f_s\right)$.

where $\|.\|$ is the matrix norm. Then $|Z(\varphi_1)-Z(\varphi_2)| \leq H |\varphi_1-\varphi_2|$ where the constant $H = \max(M, \|B\|)$. Therefore $Z$ is Lipschitz continuous function. As the function $\psi$ is uniformly Lipschitz and continuous then there exists a unique solution for the system (4.3). First, introduce the Lagrangian for the optimal problem (4.1) and (4.2)

We derive necessary conditions on the optimal control by applying Pontryagin's Maximum Principle ^[23]. Let the Hamiltonian $H$ for the control problem (4.1) and (4.2):

where $\lambda_1, \lambda_2$ and $\lambda_3$ are the adjoint variables satisfying the following adjoint equations

Final conditions are given as follows: $\lambda_1(T) = 0, \lambda_2(T) = 0, \lambda_3(T) = 0$.

The Hamiltonian is minimized with respect to the control variables as following:

The derivative of the Hamiltonian $H$ with respect to $v_1$, and $v_2$ is given by

The optimal control is obtained by resolving the necessary conditions: $\dfrac{\partial H}{\partial v_1} = 0$ and $\dfrac{\partial H}{\partial v_2} = 0$ on some non-trivial intervals. In this case, the controls are expressed as

once

5.   Numerical simulations

For all numerical simulations, we used the Monod function as a non-linear incidence rates satisfying the assumption 1 $f_a(I) = \dfrac{\eta_a I}{\kappa_a+I}$ and $f_s(I) = \dfrac{\eta_s I}{\kappa_s+I}$. Here $\eta_a$, $\kappa_a$, $\eta_s$ and $\kappa_s$ are non-negative constants. The parameters used for the numerical simulations are given in Table 1.

We begin by some numerical results that confirm the stability of the equilibrium points of (3.2). In Figure 2, we give the results for the case where $\mathcal{R} < 1$. The numerical solution of the given model (3.2) approaches to the $DFE = (0, 0)$, which confirms that DFE is GAS when $\mathcal{R} < 1$.

In Figures 3 and 4, we give the results for the case where $\mathcal{R} > 1$. The numerical solution of given model (3.2) approaches asymptotically to the EE, which confirms that EE is GAS when $\mathcal{R} > 1$.

Sensitivity analysis can be helpful as to how the variability in the output of a mathematical model can be allocated to different sources of uncertainty in its input parameters ^[24,25]. Sensitivity analysis has several purposes; one is to determine the input parameters that most contribute to the system's dynamics. The other is to detect the impacts of each parameter on the other parameters and then determine the potential to simplify the model. In this study, the sensitivity study can be helpful as it will inform us how essential each parameter is to the transmission of the disease. We must discover the highest effect on the $\mathcal{R}$. Therefore, these input parameters will be critical targets for future intervention strategies. A fundamental approach expresses a relative change of a variable by a relative change of a parameter. Consequently, if the sensitivity index is positive, increasing the parameter value will cause an increase in the $\mathcal{R}$ value. Further, if the result is negative, then the parameter value and the $\mathcal{R}$ are inversely proportional.

Figure 5 shows that $\Lambda, p$ and $\gamma$ are the only three parameters with a positive sign. Therefore, as the value of $\Lambda, p$ and $\gamma$ increase, the value of $\mathcal{R}$ increases. Furthermore, since the remaining parameters increases, the value of $\mathcal{R}$ decreases. Figure 5 shows the behaviour of the $\mathcal{R}$ with respect to the model parameters.

The numerical results of the control problem were obtained using the an association between the Gauss-Seidel-like implicit finite-difference scheme and the first-order backward-difference. We used the same parameters for systems (4.1) and (4.2) as for the direct problem (3.2) with $\eta_a = 1$ and $\eta_s = 2$ (left) and with $\eta_a = 1$ and $\eta_s = 0.2$ (where EE is GAS) with variables $v_1$ and $v_2$ such that the initial condition $v_1(0) = v_2(0) = 0.5$. We plot in Figure 6 the behaviours $S(t), I_a(t)$ and $I_s(t)$ with respect to time.

As seen in Figure 2, susceptible compartment increases about $21\%$, however, the infected compartment decreases about $14\%$ for the asymptomatic infected compartment and $14\%$ for the asymptomatic infected compartment of their values at steady state in Figure 2. Note that there is no considerable influence of the values of $\alpha_1, \alpha_2, \beta_1$ and $\beta_2$ on the final values of susceptible and both infected compartments.

6.   Conclusions

In the present paper, a generalized "SIS" model with a non-linear incidence rate including asymptomatic and symptomatic infection was presented and analyzed. The basic reproduction number ($\mathcal{R}$) was calculated for the proposed system using the next generation matrix method. The system admits at most two equilibria: the disease-free equilibrium ($E_0$) and the endemic equilibrium ($E^*$). Local and global stability was carried. If $\mathcal{R} < 1$, the disease-free equilibrium $E_0$ is locally and globally asymptotically stable; if $\mathcal{R} > 1$, the endemic equilibrium $E^*$ is locally and globally asymptotically stable. Furthermore, we find several reasonable optimal control strategies to the prevention and control the disease. Finally, numerical simulations are presented to verify the above theoretical results are given. It is concluded that asymptomatic infections have a greater contribution to the diseases spread even in the case where the asymptomatic infections are less infectious.

One of the compartmental models' limitations is the assumption of a homogeneous population. However, variations in individuals in epidemiological characteristics are more realistic to consider. Further, the influences of time delay and intrinsic fluctuations are considered limitations to the compartmental models.

Competing interests

The author declares no conflict of interest.

Acknowledgements

The author is grateful to the unknown referees for many constructive suggestions, which helped to improve the presentation of the paper.

Appendix

Used numerical scheme (control problem)

Subdividing the interval $[0, T]$ as the following $[0, T] = \bigcup_{n = 0}^{N-1} [t_n, t_{n+1}]$ where $t_n = n dt$ and $dt = T/N$. Let $S^n$, $I_a^n$, $I_s^n$, $\lambda^n_1$, $\lambda^n_2$, $\lambda^n_3$, $v_1^n$ and $v_2^n$ approximate $S(t), I_a(t), I_s(t), \lambda_1(t), \lambda_2(t), \lambda_3(t)$, $v_1(t)$ and $v_2(t)$, respectively at the time $t_n$. $S^0, I_a^0, I_s^0$, $\lambda^0_1, \lambda^0_2, \lambda^0_3, v_1^0$ and $v_2^0$ be their values at initial time ($t = 0$). $S^N, I_a^N, I_s^N$, $\lambda^N_1, \lambda^N_2, \lambda^N_3, v_1^N$ and $v_2^N$ be the values at $t = T$. Then, we use the following scheme ^[6,8,21,22]

Therefore the following algorithm will be applied.